Dual Regulatory Mechanisms of Expression and Mutation Involving Metabolism-Related Genes FDFT1 and UQCR5 during CLM.

Colorectal cancer (CRC) is the third most common cancer worldwide, and liver metastasis presents a major cause of CRC-associated death. Extensive genomic analysis has provided valuable insight into the pathogenesis and progression of CRC; however, a comprehensive proteogenomic characterization of CRC liver metastasis (CLM) has yet to be reported. Here, we analyzed the proteomes of 44 paired normal colorectal tissues and CRC tissues with or without liver metastasis, as well as analyzed genomics of CRC characterized previously by The Cancer Genome Atlas (TCGA) to conduct integrated proteogenomic analyses. We identified a total of 2,170 significantly deregulated proteins associated with CLM, 14.88% of which were involved in metabolic pathways. The mutated peptide number was found to have potential prognosis value, and somatic variants revealed two metabolism-related genes UQCR5 and FDFT1 that frequently mutated only in the liver metastatic cohort and displayed dysregulated protein abundance with biological function and clinical significance in CLM. Proteogenomic characterization and integrative and comparative genomic analysis provides functional context and prognostic value to annotate genomic abnormalities and affords a new paradigm for understanding human colon and rectal cancer liver metastasis.

Introduction

Colorectal cancer (CRC) is a significant contributor of cancer morbidity and mortality.1, 2, 3 Almost half of CRC patients die within 5 years of diagnosis due to the development of recurrent disease and metastasis.4, 5, 6, 7, 8 Therefore, it is important to illuminate the molecular basis of CRC liver metastasis (CLM) in hopes of developing new effective treatment modalities. The Cancer Genome Atlas (TCGA) has characterized the genomic features of many types of human cancers, including CRC,9, 10, 11 and the Clinical Proteomic Tumor Analysis Consortium has also performed CRC-integrated proteomic analyses.12, 13, 14, 15 However, the primary genetic basis of CLM has not been fully elucidated. Understanding the genetic and proteogenomic differences between primary colon cancer and associated metastases to the liver is essential for discovering metastasis-specific molecular biomarkers and for devising a better therapeutic approach for this disease.

Results and Discussion

To identify protein variants associated with CLM, we collected tumor samples and paired normal colorectal tissues (PNs) from 44 CRC tumor samples with (MT, n = 23) or without liver metastasis (NM, n = 21) in comparison with PNs (Table S1) to perform liquid chromatography-tandem mass spectrometry (LC-MS/MS)-based shotgun proteomics profiling (Figure S1). A total of 12,177,738 spectra were used in the Andromeda engine search, and 63,720 unique peptides were identified in an assembly of 5,758 protein groups with a protein-level false discovery rate (FDR) of 1.0%. Ingenuity pathway analysis with all 5,758 identified proteins showed that about 55% of the proteins were from the cytoplasm, 28% were from the nucleus, 7% were from the plasma membrane, and 2% were from the extracellular space, whereas 8% of proteins remained unclassified (Figure S2A). The random predicted cellular distribution of the proteins supports the quality of the sample preparation. A scatterplot of protein abundance between CRC tissues and PNs showed that there was a great variation between the NM or MT tumors and PN (Figures S2B and S2C). However, the protein expression between the MT and NM group was positively correlated (R2 = 0.86) (Figure S2D). These results suggest that metastatic and non-metastatic CRCs share similar protein profiles and that there are common molecular alterations at each stage of tumor development.16, 17, 18

Among the 5,758 proteins, a total of 1,679 proteins were significantly altered between CRC tissue and normal colorectal tissues (Figure 1A, left bars). There were 2,170 proteins with significant differences in MT tissues when compared with PNs or NM tissues (Figures 1A, middle and right bars, and 1B). To explore the functions of proteins that are dysregulated in CLM, we used DAVID (Database for Annotation, Visualization and Integrated Discovery) analysis software to classify the gene ontology of the 2,170 significantly altered proteins in MT tissues (Figures S3A–S3C). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed proteins revealed that metabolic pathways were altered in the MT samples, with 323 significantly differentially expressed proteins (14.88%) involved in metabolic pathways (Figure S3D), which suggests that metabolism-related pathways may play important roles in the liver metastasis of CRC.21, 22, 23, 24

Figure 1

Differently Expressed Protein and mRNA-Protein Correlation Analysis

(A) Numbers of differently expressed proteins (≥2-fold difference; p < 0.05 with a FDR q value < 0.05) between the 44 CRC and paired PN samples, and the 23 MT and 21 NM samples. (B) Significantly changed proteins among three groups from LC-MS/MS data. (C) KEGG pathway analysis of the 2,170 differentially expressed proteins. (D) Significantly changed genes among three groups from RNA-sequencing data. (E) 444 significantly changed genes showed significant mRNA-protein correlation, with a mean Spearman’s correlation coefficient of 0.71. Among these, 57 genes were enriched in metabolism pathways.

We performed unsupervised hierarchical clustering of the expression data from RNA-sequencing data to evaluate the genetic diversity between the primary and metastatic CRCs, based on the assumption that genetically similar CRCs and their matched metastatic foci would be closely related.25, 26, 27 As expected, the results showed that the MT and LM tissues had closely related expression profiles, indicating clonal and genetic similarity for these pairs. However, the NM samples were more distantly related, suggesting distinct genetic relationships between patients with and without metastasis (Figure S4A). A total of 2,299 genes from the RNA-sequencing were significantly changed in CLM (Figure S4B). KEGG pathway classification enrichment analysis indicated that 915 genes (39.8%) were enriched in metabolism pathways (Figure 1C), which was similar to the findings of the protein analysis (Figure S3D). These results provide additional confirmation that metabolism pathways may play a role in the carcinogenesis and development of CRC.

For further confirmation of these results, we analyzed RNA-sequencing data from TCGA portal. A total of 6,585 significantly changed genes in CRC tissues were identified (Figure 1D). When compared with the 2,170 significantly changed CLM-associated proteins identified by LC-MS/MS, 444 significantly changed genes showed significant positive mRNA-protein correlation (Figure 1E, left two panels). To determine whether the concordance between protein and mRNA variation is related to the biological function of the gene product, we performed KEGG enrichment analysis, which indicated that among the 444 significantly deregulated genes or proteins, 57 are enriched in metabolic pathways (p = 2.38E−7) (Figure 1E, right panel). These findings further verify the role of metabolic pathway genes in CLM.

To examine the clinical significance of our findings, we evaluated the clinical characters data of 374 CRC patients from TCGA database for the 444 significantly deregulated genes or proteins (Table S2). Our results showed that 17 dysregulated genes were significantly associated with overall survival (OS) and progression-free survival (PFS) of CRC patients. Among them, high expression of four upregulated genes, BGN, EHD2, UQCR5, or COL1A2, was associated with poor survival of CRC patients.29, 30, 31, 32 Furthermore, high expression of FDFT1 significantly prolongs the OS and PFS of CRC patients relative to low expression of FDFT1 (Table S3).

A fundamental goal of proteogenomics is to identify protein-coding alterations that are expressed at the protein level.33, 34, 35 However, standard database search approaches cannot identify variant peptides from MS/MS data.36, 37, 38 Therefore, we created a customized mutation database to search for single amino acid variants (SAAVs) in CRC. A SAAV library was prepared using 113,844 mutated sites in CRC tissues from cBioport, and 16,581 mutated proteins were identified, which constitutes 82.08% of 20,201 proteins in the human protein library (Figure 2A). We determined the total numbers of mutated and non-mutated peptides and tumor-specific mutant peptides (Figure 2B) and found that mutated peptide numbers in MT samples were significantly increased (Figure 2C), which indicates that the mutated peptide number has potential predictive value for CLM.

Figure 2

Numbers of SAAVs in Paired PN, NM, or MT Samples

(A) The proportion of mutated proteins and amino acids in CRC samples were calculated by comparing LC-MS/MS data for the standard protein library and SAAV library. (B) The mutated and non-mutated peptides numbers of 44 paired CRC tissues. (C) Numbers of SAAVs in 21 NM, 23 MT, and their PNs. (D) Numbers of NM-specific, MT-specific, and NM- and MT-shared SAAVs. The mutated peptides were identified by comparing LC-MS/MS data for the standard protein library and SAAV library.

Among those, 140 SAAVs in 131 proteins occurred only in NM patients (Figure 2D; Table S4), and 223 proteins in 18 MT patients had 256 SAAVs, of which 110 SAAVs in 100 proteins occurred in both NM and MT samples (Figure 2D; Table S5), and 203 SAAVs in 184 proteins only occurred in MT samples (Figure 2D; Table S6). Similarly to single nucleotide variants (SNVs), some mutations reported previously, such as those in TP53, antigen-presenting cells (APCs), vascular endothelial growth factor (VEGF), and SMAD4 were found.39, 40, 41, 42, 43 Among these, two mRNA-protein positively changed and metabolism-related proteins, FDFT1, a 47-kDa membrane-associated enzyme located at a branch point in the mevalonate pathway involved in the replication stage of the hepatitis virus C (HCV) life cycle and paclitaxel sensitivity in hypopharynx cancer cell, and UQCR5, a component of the ubiquinol-cytochrome c reductase complex, amplifying in primary breast cancer core biopsy samples46, 47, 48 and overexpressing in gastric cancer, were found to be an occurrence of somatic alterations, which was validated with Sanger sequencing, in the 18.2% (8/44) and 9.1% (4/44) of CRC, and are specific to CLM (Figure 3A). Mutation in FDFT1 was predicted to lose the function; however, mutation in UQCR5 was predicted to gain the function, and increased copy number enhanced UQCR5 gene expression and lead to increased protein expression when scale-invariant feature transform (SIFT) and PolyPhen2 were used for predicted function of FDFT1 and UQCR5 mutations (Figure 3B).

Figure 3

The Frequency and Distribution of Mutational Proteins Specifically Altered in MT for 44 CRC

(A) Heatmap comparing the frequency and distribution of 12 mutational proteins in 44 CRC (including 23 MT and 21 NM). Red, amplification or high expression; green, deletion or low expression. (B) Scaled probabilities for entire protein score of predicted mutation using SIFT and PolyPhen2 for FDFT1 and UQCR5 mutations.

Moreover, FDFT1 knockdown (Figures S5A and S5B) or UQCR5 overexpression (Figures S5C and S5D) led to a significant increase of migrated distance (Figure 4A) and an obvious increase in matrigel invasion (Figure 4B). FDFT1 was found to downregulate in CRC tissues when compared with normal colorectal tissues (Figure 4C). Determination of UQCR5 gene expression identified overexpression (Figure 4D) and negatively correlated with FDFT1 expression in normal colorectal tissues and CRC tissues (Figure 4E).

Figure 4

The Biological Function and Clinical Significance of FDFT1 and UQCR5 in CRC

Wound-healing assay (A) and migration abilities (B) of the parental and shFDFT1 or UQCR5 overexpressed SW480 cells. LC-MS/MS to quantify FDFT1 (C) and UQCR5 (D) levels in 44 paired CRC (including 23 MT and 21 NM) and normal colorectal tissues. (E) The expression correlation between FDFT1 and UQCR5 expression levels in 44 adjacently normal colorectal tissues or CRC tissues.

Patients with lower FDFT1 expression were found to have shorter median OS (Figures S6A and S7A) and PFS (Figures S6B and S7B), and high UQCR5 expression was an important risk factor for OS (Figure S6A and S7A) and PFS (Figures S6B and S7B). Moreover, concomitant low expression of FDFT1 and high expression of UQCR5 correlated with a shorter median OS (Figures S6A and S7A) and PFS (Figures S6B and S7B) in CRC patients. These studies provide direct evidence for contribution of metabolism-related genes.

To the best of our knowledge, this is the first comprehensive study to use proteogenomic profiling of primary CRCs from patients with or without liver metastasis to define the dominant events of metastatic lesions in terms of their expression and mutation. Our comprehensive integrative analysis of 44 colorectal tumors and normal pairs provides a number of insights into the biology of CLM and identifies potential therapeutic targets. Moreover, our characterization of the annotated metastatic CRC proteome clarifies the power of integrating genomics (SNVs) and proteomics (SAAVs). This approach provides new insights into the roles of these protein alterations in CLM, which can be broadly extended to understand the roles of protein mutation in other cancers.

Materials and Methods

Details of sample preparation and analysis are described in the Supplemental Materials and Methods. The study was examined and approved by the Ethics Committee of the Shanghai Tenth People’s Hospital, Tongji University School of Medicine. This study was registered with ClinicalTrials.gov: NCT02917707). Forty-four paired fresh primary tumors from the colorectal and PNs were used for proteogenomic analysis. The reliability of the exome analysis and somatic variant identification strategies was assessed using PCR and Sanger sequencing. Nine specimens from three CRC patients with metastasis and three specimens from three CRC patients without liver metastasis were obtained for RNA-sequencing analysis. Human CRC cell line SW480 was used for scratch-wound and in vitro invasion assays. All calculations were performed with SPSS 20.0 software.

Whole-exome sequencing data from this study are available for download through NCBI Sequence Read Archive: PRJNA358865. All RNA-sequencing data have been deposited in the Gene Expression Omnibus: GSE92914. All of the MS proteomics data have been deposited to iProX (https://www.iprox.org/index): IPX00083203 and IPX00083210.

Author Contributions

Y.-S.M., J.-B.L., X.-M.W., and D.F. designed and supervised the research and experiments. Y.-S.M., Z.-J.W., B.C., H.-W.Z., T.H., H.-D.L., H.X., Y.-Z.Z., Y.-Z.Y., R.-T.X., L.L.T., and D.F. analyzed the genetic data. Y.-S.M., Z.-J.W., B.C., H.-W.Z., T.H., S.-B.X., L.L., C.-L.L., J.-B.L., X.-M.W., and D.F. performed the proteogenomic experiments and analyzed the data. Y.-S.M., Z.-J.W., B.C., H.-W.Z., T.H., S.-B.X., L.L., R.-T.X., J.-B.L., and X.-M.W. performed tissue acquisition and clinical data collection. Y.-S.M., Z.-J.W., B.C., H.-W.Z., S.-B.X., L.L., R.-T.X., J.B.L., and D.F. performed functional validation experiments. All authors contributed to the final version of the manuscript.

Conflicts of Interest

The authors declare no competing interests.