The ES-2 human ovarian clear cell carcinoma cell line and HEK293T cell line were purchased from the American Culture Collection (ATCC) and cultured in Dulbecco’s modified Eagle’s medium (DMEM, Gibco, Grand Island, NY, USA) containing 10% Fetal Bovine Serum (FBS, HyClone, Logan, UT, USA) and 1% penicillin/streptomycin (P/S) (Gibco, Grand Island, NY, USA) at 37 °C with 5% CO${}_2$. Lentiviral packaging plasmids (Pspax2, PMD2G and ZsGREEN1 GFP, Thermo Fisher, Waltham, MA, USA) were transfected into HEK293T cells according to the manufacturer’s instructions to generate lentivirus, and a stable ES-2 cell line expressing GFP (ES-2-GFP) was obtained by monoclonal screening in 96-well plates after lentiviral transduction.

NCG (NOD/ShiltJGpt-Prkdcem26Cd52Il2rgem26Cd22/Gpt) strain mice, the severe immunodeficient strain mice lacking mature T cells, B cells, and NK cells used in the experiment (female, SPF grade, $>$8 weeks old) were purchased from JiCui Pharmaceutical Co., Ltd. All NCG mice were housed in a specific pathogen-free (SPF) animal room with free access to food and water. All animal experiments were performed and approved by the Animal Care and Use Committee of West China Hospital, Sichuan University. Exponentially growing ES-2-GFP cells were trypsinized to prepare a cell suspension in PBS at a concentration of 1 $\times$ 10${}^7$ cells/mL, and the multiorgan metastatic model was induced in NCG mice by i.v. injection of 1 $\times$ 10${}^6$ ES-2-GFP cells suspension per mouse.

The mice were killed and left open by cervical dislocation between 18 and 20 days postinoculation. The distribution of metastatic lesions in each organ of 12 mice was observed and then validated by fluorescence imaging. The static dissection of each organ was loaded into a frozen storage tube and frozen in liquid nitrogen.

ES-2 cells were washed twice with cold Phosphate buffered saline (PBS) and lysed in 1000 $\mu$L SDT lysis buffer containing 2% sodium dodecyl sulfate, 100 mM Tris/HCl pH 7.4, 100 mM dithiothreitol (Sigma‒Aldrich, St. Louis, MO, USA) and 1 $\times$ protease cocktail inhibitors (Roche Diagnostics, Sandhofer Strasse, Mannheim, Germany). After lysis, the mixture was transferred to a 1.5 mL centrifuge tube using cell scraping and sonicated on ice for one minute (5/5 s, 6 cycles, output 60 W, frequency 45 $\pm$, SCIENTZ-IID, Ningbo, China). The supernatants were collected after centrifugation at 14,000 g and 4 °C for 10 min. Approximately 100 mg of tissue was extracted from each metastasis sample, homogenized in 800 $\mu$L SDT lysis buffer, sonicated and centrifuged to collect supernatants under the same conditions. The protein concentration of each sample was determined by the tryptophan fluorescence emission method at 365 nm. Proteins in each sample (150 $\mu$g) were precipitated with cold acetone for 3 h at –20 °C. After centrifugation of the precipitated protein, the liquid layer was removed, and the cabinet was dried for 5 min. The protein precipitate was dissolved by the addition of 7 M guanidine hydrochloride and then incubated with 45 mM dithiothreitol for 30 min at 37 °C for reduction and alkylation and 106 mM iodoacetamide for 30 min at 25 °C in the dark. Then, the samples were transferred to a Vivacon 500 30 kDa Molecular weight cutoff (MWCO) centrifugation device (Sartorius, Gottingen, Germany), washed twice with an additional 100 $\mu$L 8 M urea solution and sequentially washed twice with 100 $\mu$L 50 mM ammonium bicarbonate. The samples were first digested with sequencing-grade trypsin (Promega, Madison, WI, USA) at a trypsin to protein ratio of 1:50 for 16 hours at 37 °C. After digestion, the peptides were dried in a vacuum and stored at –80 °C until use. All water used in the experiments was purified from a Milli-Q system (Millipore, Bedford, MA, USA), and other highest grade reagents and chemicals were purchased from Sigma‒Aldrich or Thermo Fisher Scientific.

Sequential window acquisition of all theoretical mass spectra (SWATH-MS) analysis was performed on a quadrupole time-of-flight spectrometer (TripleTOF5600, AB/Sciex, Framingham, MA, USA) coupled to a micro-HPLC (NanoLC 415, Eksigent, Framingham, MA, USA). The LC flow phase consisted of mobile phase A (2% acetonitrile and 0.1% FA in ultrapure water) and mobile phase B (98% acetonitrile and 0.1% FA in ultrapure water) (all purchased from Thermo Fisher Scientific, Waltham, MA, USA). Samples were first captured using a ChromXP c18 precolumn (5 $\mu$m particle size, 120 Å, 75 $\mu$m id $\times$ 5 mm). Impurities and salt were removed for 6 min at a flow rate of 100% buffer a for 6 L/min. The peptides were then eluted onto a ChromXP c18 analysis column (3 µm particle size, 120 Å, 75 µm id × 15 cm). Mobile phases A and B were used to create the separation gradient of the eluted peptides in the column by increasing the solvent B concentration. A flow rate of 5 L/min was maintained throughout the flow phase, and the following gradient elution procedure was run: 5% B for 0.5 min; 5% to 25% B for 45 min; 25% to 35% B for 10 min; 35% to 55% B for 5 min; 55% to 80% B for 0.5 min; 80% B for 4.5 min; 80% to 5% B for 0.5 min; 5% B for 4.5 min. The key ion source parameters for data acquisition in the SWATH-MS mode are as follows: ion spray voltage float (ISVF) at 5500 V, curtain gas (CUR) at 30 psi, ion source gas 1 (GS1) at 17 psi, ion source gas 2 (GS2) at 13 psi, temperature (TEM) at 330 °C, and decluster potential (DP) at 100 V. The mass spectrometer was operated in cycle product ion mode. Full-MS scans were performed in the mass range of 350–1500 m/z in positive-ion mode. One hundred MS/MS acquisition windows with an overlap range of 1 m/z were constructed using a variable window mode. The coverage mass range was 100–1500 m/z in high sensitivity mode with a CES of 15 V. The maximum fill time for the MS scan was 250 ms and 30 ms for the MS/MS scan, with a final cycle time of 3.04 s. Peptide samples were prepared as 0.5 g/L, mixed with 12.5 fmol/L peptides from trypsin-digested bovine serum albumin (BSA), and analyzed for 5 L for each sample.

For HR-MRM experiments, LC‒MS/MS was performed on a quadrupole time-of-flight spectrometer (TripleTOF5600, AB/Sciex, Framingham, MA, USA) coupled to a micro-HPLC (NanoLC 415, Eksigent, Framingham, MA, USA). The LC flow phase consisted of mobile phase A (2% acetonitrile and 0.1% FA in ultrapure water) and mobile phase B (0.1% FA and 2% ultrapure water in acetonitrile) (all purchased from Thermo Fisher Scientific, Waltham, MA, USA). The sample was first loaded onto a ChromXP c18 precolumn (5 $\mu$m particle size, 120 Å, 75 $\mu$m id $\times$ 5 mm). Impurities and salts were removed at a flow rate of 6 $\mu$L/min in 100% buffer A for 6 minutes. The peptides were then eluted onto a ChromXP c18 analysis column (3 µm particle size, 120 Å, 75 µm id × 15 cm). Mobile phases A and B were used to create the separation gradient for eluting peptides in the column at an increasing concentration of solvent B. A flow rate of 5 $\mu$L/min was used throughout the flow phase, and the following gradient elution procedure was run: 5% to 8% B for 0.5 min; 8% to 30% B for 45 min; 30% to 50% B for 5 min; 50% to 80% B for 0.5 min; 80% B for 4.5 min; 80% to 5% B for 0.5 min; 5% B for 4.5 min. The key ion source parameters for data acquisition in the HR-MRM mode were as follows: ion spray voltage float (ISVF) at 5500 V, curtain gas (CUR) at 30 psi, ion source gas 1 (GS1) at 18 psi, ion source gas 2 (GS2) at 15 psi, temperature (TEM) at 350 °C, and cluster potential (DP) at 100 V. Individual candidate peptide TOF-MS scans were set at their m/z value followed by a MS/MS product ion scan from 100 to 1500 Da. The final cycle time was 2.16 s. Data acquisition was performed using Analyst software (version 1.6.3, Sciex-Foster, CA, USA).

Mass spectrometry acquisition data were processed using Protein-Pilot software (AB/Sciex, v.5.0.1, Framingham, MA, USA). The nonredundant human UniProtKB/Swiss-Prot protein database (https://www.uniprot.org/uniprotkb?facets=reviewed%3Atrue&query=proteome%3AUP000005640), BSA peptide sequence and mouse UniProtKB/Swiss-Prot protein database (https://www.uniprot.org/uniprotkb?facets=reviewed%3Atrue&query=proteome%3AUP000000589) were combined to form a heterozygous library (https://www.uniprot.org/uniprotkb/P02769/entry) for IDA data alignment to identify proteins. Search parameters were set as follows: sample type: identification; Cys alkylation: iodoacetamide; digestion: trypsin; instrument: TripleTOF 5600; special factors: none; species: Homo sapiens; search effort: thorough ID; results quality: 0.05 (detected protein threshold) performed Paragon searches. The false discovery rate (FDR) for proteins and peptides was set to 0.01 (minimum length set to seven amino acids) for analysis. Data from all IDA method identification runs performed using ProteinPilot software were merged into one batch and used for SWATH assay library construction. After database retrieval, the mouse polypeptides were removed by using the self-compiled R language software package (R 3.6.0 programming language, R Foundation for Statistical Computing, Vienna, Austria) to establish a human protein quantification database. Protein quantification information was extracted using PeakView software 2.2 and the swath plugin, where the BSA peptides were used to calibrate the retention time calibration. To obtain better quantitative results using the SWATH acquisition method, the following parameters were set: peptide filter: 15 peptides per protein; six transitions per peptide; 95% peptide confidence; 1% FDR threshold; exclusion of modified peptides; XIC option: 8 min XIC extraction window; 50 ppm XIC width. Missing values in the results are automatically populated by the PeakView software (AB/Sciex, version 2.2, Framingham, MA, USA) built-in SWAP plugin. Subsequent statistical analysis of the data was performed using R statistical software, version 3.6.0. The median method was chosen to normalize the data, and the p value of protein expression was calculated with the Kruskal‒Wallis test. A Benjamini‒Hochberg adjustment p value 0.05 was used to screen for differentially expressed proteins (DEPs).

Bioinformatics analysis was performed using the Gene Ontology (GO) database (http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genome (KEGG) (http://www.genome.jp/kegg/pathway.html) based on the DEPs. In addition, Ingenuity Pathway Analysis (IPA) (http://www.ingenuity.com/products/ipa) was used to determine the canonical pathways among DEPs. Fisher’s exact test was used to calculate the p values.

Considering the inevitable mixing of human tumor cells with murine organ tissue during tumor sampling, we designed primers for differential bases in human and mouse mRNAs, which is used to solve the problem that the human metastasis samples cannot be completely separated from the tissues of NCG mice due to surgical operation during the RT‒PCR experiments. Primer synthesis was commissioned by Beijing Qingke Biological Company. The PCR primer sequences are shown in Supplementary Table 1.

Total RNA from ES-2 cells and metastatic foci of various organs was extracted using TRIzol reagent (Gibco, Life Technologies, Carlsbad, CA, USA) and reverse-transcribed into cDNA using a RevertAid™ First Strand cDNA Synthesis Kit with a DNase I Reverse Transcription Kit (ThermoFisher, USA). PCR amplification was performed using the TransStart Taq DNA Polymerase PCR Kit (TransGen Biotech, Inc., Beijing, China). Fluorescence quantitative PCR experiments were performed using a CFX96 dual-channel real-time PCR instrument (Bio-Rad, Hercules, CA, USA) with 2$\times$ TSINGKE® Master qPCR Mix (SYBR Green I) (TSINGKE, Beijing, China). The cycling conditions were as follows: 50 °C for 2 min, followed by 35 cycles of 95 °C for 2 min, 95 °C for 15 s, and 60 °C for 30 s. Data were analyzed by Bio-Rad CFX Manager software. The relative expression levels of the target genes versus the reference gene (GAPDH) were calculated using the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{CT}}$ method. All experiments were performed in triplicate, providing three technical repeats.

Quantitative data are represented as the mean $\pm$ standard deviation. Data analysis was performed using the one-way ANOVA method in GraphPad Prism 8.0.2 software (GraphPad Software Inc., San Diego, CA, USA). A value of p 0.05 was considered statistically significant.

To facilitate observation and verification of metastases, ES-2 cells were labeled using lentiviral transfection with green fluorescent protein (GFP), and distinctly fluorescent green was observed under fluorescence microscopy. At 18 to 22 weeks after cell injection, mice were sacrificed, and multiorgan metastases of ES-2 cells, such as in the liver, lung, ovary and peritoneum, were clearly observed in NCG mice. In addition, fluorescence imaging proved that the number of metastases formed by the ES-2 cell lines labeled with GFP varied significantly in different organs (Fig. 2A–C).

Fig. 2.

Tissue harvesting and fluorescence imaging of ES-2 cell line. (A) Fluorescent imaging of the metastasis formed by ES-2 cells. (B) In vivo optical imaging of liver, lung, and ovary metastases formed by ES-2 cells. (C) In vivo optical imaging of connective tissue metastases formed by ES-2 cells. The white arrowhead indicates metastases. The connective tissues in different organs are collectively referred to as Other.

From the data acquired by the SWATH strategy, a total of 6023 proteins were identified from 32,602 different unique peptides by matching to the heterozygous library containing human and mouse protein sequences and bovine serum albumin peptide sequences. After excluding 17,809 mouse source peptides, the number of remaining human peptides was 14,793, from which a total of 4503 human proteins were identified. Ultimately, a human protein quantitative database was established. p values were calculated using the Kruskal‒Wallis test, which were corrected using the Benjamini‒Hochberg strategy with an FDR 0.1 (p. adjusted), and a total of 2868 proteins with a p value 0.05 after correction were screened.

To understand the gene regulation and protein profile expression changes from circulating ES-2 cells in blood vessels to different organs in colonization and metastasis, we first performed a hierarchical clustering analysis (HCA) on protein data from ES-2 cell lines and the metastases of different organs. From the heatmap (Fig. 3A), the tumors metastasized to organs showed significant differences in protein expression profiles compared with the ES-2 cell lines. Moreover, various profiles were also clearly observed among the metastases from different organs. Subsequently, a principal component analysis (PCA) of the protein expression data of different metastasis groups was carried out, and the results showed obvious distinctions in the distribution of these metastatic samples, each with significant organ-dependent separation. Interestingly, the protein spectrum of liver metastases showed the most pronounced deviation compared with metastases from other organs (Fig. 3B). Next, we performed statistical analysis of differential proteins between the ES-2 cell line and each organ metastatic tumor group separately. According to the screening criteria of FDR 0.02 and fold change (FC) $>$1.5 (volcano plots are shown in Supplementary Fig. 1), compared with the original ES-2 cell line, we obtained 1938 differential proteins in the liver metastasis group, 1499 differential proteins in the lung metastasis group, 1850 differential proteins in the ovary metastasis group and 1641 differential proteins in the peritoneum metastasis group (Fig. 3C). To study the commonness and organ specificity of gene expression of ES-2 cells in the formation of different organ metastases, we built a Venn diagram to show the collection of proteins expressed in different organ metastases, shown as Fig. 3D. Notably, there are 1020 proteins in the intersection of multiple organ metastases. We chose the 227 proteins uniformly upregulated in various organ metastasis for GO enrichment analysis and KEGG pathway analysis and results showed these proteins are mainly located in the cytoplasm and nucleus, as well as in exocytic exosomes, mainly involved in the intracellular protein transport and hydrolysis process, and their molecular functions focus on the cell adhesion with ATP binding and cadherin participation (Supplementary Fig. 2A). At the same time, these differential proteins were mainly enriched in metabolism as well as oxidative phosphorylation pathways, including the active glycolysis of tumor cells in order to adapt to the metabolic stress in the microenvironment, as well as the nucleotide metabolic pathways (Supplementary Fig. 2B). These results implied the proteins in the intersection of multiple organ metastases might be associated with general gene expression regulation in solid metastatic tumor formation from tumor cells. On the other side, Venn diagram also reflected the organ specificity of gene expression of ES-2 cells in different organ metastases: the mounts of proteins specifically regulated in metastases of different organs are separately 158 for liver, 80 for lung, 149 for ovary and 55 for other connective tissues, and these proteins may be involved in tumor colonization in corresponding target organs.

Fig. 3.

Multivariate statistical analysis of the total protein expressed by different organs metastasis and ES-2 cell lines. (A) Hierarchical cluster analysis of the 2868 DEPs. Heatmap showed that clustering of differential proteins between ES-2 cells and different organ metastases. (B) Principal Component Analysis of the remanent DEPs from metastases in different organs after removing the proteins expressed by the cells. PCA analysis showed that there was a significant outlier trend of proteins expressed by liver metastases. The p-value was calculated using the Kruskal-Wallis test. A Benjamini-Hochberg adjustment p-value 0.05 (FDR 0.1) was used to screen for DEP. (C) Distribution of different protein quantities expressed by different organ metastases. (D) According to the screening criteria of FDR 0.02 and fold change (FC) $>$1.5, Veen plot shows the intersectional relationship of proteins expressed by different organ metastases.

Since in all metastases, liver metastases showed the most significant differences in protein expression compared to original ES-2 cells and showed obvious deviation from other organ metastases in PCA scoring, we focused on the 158 proteins specifically expressed in liver metastases (Listed in Supplementary Table 2) as representatives in subsequent bioinformatical studies. We first performed GO analysis and KEGG pathway analysis to understand the protein expression changes observed in liver metastases. GO is used to describe gene functions and is divided into three categories: biological processes, cellular components, and molecular functions. GO enrichment analysis showed that 158 proteins specifically expressed in the liver were mainly located in the cytosol and nucleus and were mainly involved in DNA transcription, and molecular function focused on ATP binding (Fig. 4A). Further KEGG pathway analysis of the above 158 proteins showed that these proteins were mainly enriched in metabolic pathways (Fig. 4B).

Fig. 4.

Summary of the 158 proteins specifically expressed in liver metastases. (A) GO annotation of the 158 proteins specifically expressed in liver metastases. Only the top 10 words are shown according to the order of the p-adjusted. BP, biological process; CC, cellular component; MF, molecular function. (B) KEGG pathway enrichment analysis of 158 proteins specifically expressed by liver metastases. (C) Diagram of classical overlapping pathway analysis of 158 protein genes specifically expressed in liver metastases (Classical overlapping pathways share one or more genes by interconnection, where the deeper the red, the more significant the contribution).

Next, we introduced these protein genes into the IPA for classical overlapping canonical pathway analysis and screened genes related to the liver and ovarian cancer background as well as tumor colonization and metastasis. Pathway analysis showed that the genes specifically expressed by liver metastases were mainly concentrated in the oxidative phosphorylation and glutathione biosynthetic pathways (4C). Based on IPA pathway analysis results and SWATH quantitative data, we screened 17 proteins that were upregulated in liver metastases with functional importance and sufficient MS response intensity for further verification by targeted quantification using HR-MRM acquisition. Six proteins (Listed in Table 1), including ferritin light chain (FTL), lactate dehydrogenase A (LDHA), eukaryotic translation initiation factor 2 subunit 1 (EIF2S1), carbamoyl-phosphate synthase (CPS1), long-chain-fatty-acid–CoA ligase 1 (ACSL1) and glutamate–cysteine ligase catalytic subunit (GCLC), which could be effectively quantified by HR-MRM acquisition and showed exact upregulation in liver metastases, were screened out (Fig. 5A), and this result further narrowed the range of target proteins.

Fig. 5.

RT-qPCR results. (A) Exact upregulation of six proteins quantified by HR-MRM acquisition. (B) PCR preliminary experiment of specific primer for target protein gene (ACSL1, LDHA and FTL), negative control for normal liver tissue of NCG mice. (C) Real-time PCR analysis of genes upregulated in liver metastases (ACSL1, LDHA, and FTL) expression in various organ metastases of NCG mice (n = 3). The positive control (POS) was the ES-2 cell line, and the negative control (CON) was the liver tissue of normal NCG mice, GAPDH was used as a reference gene. *p 0.05, **p 0.01, ***p 0.001.

To distinguish between human and mRNAs, we designed specific primers aimed at differential bases for distinguishing human- and murine-derived cDNA and performed PCR preexperiments. PCR experiments demonstrated that this series of primers are targeted to human-derived cDNA (Fig. 5B). To examine the reliability of the proteomic data in the study, six proteins that were completely upregulated in liver metastases resulting from the HR-MRM validation were selected for RT‒PCR validation. The results showed that the three differential protein gene expression levels derived from the six candidate target proteins were consistent with the proteomic data, including ACSL1, LDHA, and FTL, which were significantly higher in liver metastases than in lung, ovary, and other metastases (Fig. 5C).