Introduction
For patients with colorectal cancer (CRC), the overall survival benefits of systemic treatments, including tumor resection, 5-fluorouracil (5-FU)-based chemotherapy and targeted therapies, have been firmly established (1,2). Nevertheless, recurrence and metastasis caused by the drug resistance of tumor cells continue to obstruct significant therapeutic efficacy.
The development of CRC is a multi-gene, multi-step and multi-factor interactive process. The improvement of therapeutic efficacy significantly relies on improved studies of gene functions and mechanisms (3). It has been widely recognized that abnormal bioenergetics is one of the most common phenotypes of a number of tumor cells, including CRC (4). Changes in ATP synthetase are found in tumor cells and are considered to be associated with the energy metabolism and drug resistance within the tumor (5,6).
ATP synthase-coupling factor 6, mitochondrial (ATP5J) is a protein connecting F0 and F1, which are two components of ATP synthetase (7). ATP5J has been hypothesized to recover the activity of ATPase inhibited by oligomycin, as well as the exchange between ATP and inorganic phosphorus (8). Our previous study reported that ATP5J may be a tumor biomarker. Higher expression of ATP5J was found in CRC tissue compared with that in normal tissue, and in metastatic lymph nodes compared with that in primary tumors (9). In vitro experiments demonstrated that upregulation of ATP5J expression in DLD1 cells enhances cell migration and 5-FU resistance, while downregulation reverses this reaction. Investigation of ATP5J has been mainly focused on cardiovascular disease. Only a few studies have discussed the role of ATP5J in cancer, including renal carcinoma and hepatocellular carcinoma (10,11). The molecular mechanism of ATP5J in CRC cells remains unknown.
Based on previous studies, the following study investigated differentially expressed mRNAs (DEmRNAs) and lncRNAs (DElncRNAs) using microarray and bioinformatics methods. The obtained results may provide valuable information on the function and mechanism of ATP5J in CRC, and a better understanding may be beneficial for the improvement of CRC management.
Materials and methods
Cell lines
In our previous study (9), the human colon cancer DLD1 cell line was used for cell transfection and stable colony selection. The expression of ATP5J was highly downregulated in clone 4, while it was highly upregulated in clone 2. Consequently, 4 DLD1 cells that were stably transfected with pcDNA3.1(+)/ATP5J plasmid (DLD1/A2), ATP5J short hairpin (sh)RNA plasmid (DLD1/SA4) and the corresponding control vectors (DLD1/C6 and DLD1/CA6), respectively, were selected and used in the present study.
RNA extraction and microarray
Total RNA was extracted from the 4 cell types by TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA) according to the manufacturer's protocols. RNA quality was detected by an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). The microarray analysis was performed by CapitalBio Corporation (Beijing, China). Total RNA was extracted from the cells and the cDNA that was then produced by reverse transcription served as the template to synthesize the fluorescently labeled (Cy5 or Cy3) cDNA using Klenow fragment polymerase. The labeled cDNA was hybridized on the lncRNA + mRNA Human Gene Expression Microarray (4×180K; Agilent Technologies, Inc.). The microarrays were washed three times with wash solution (0.2% SDS, 2X SSC at 42°C for 120 sec) and scanned using an Agilent G2565CA Microarray Scanner (Agilent Technologies, Inc.), according to the manu facturer' s protocols. The raw data were analyzed and then normalized using percentile normalization. Only the mRNAs or lncRNAs with a fold-change ≥2.0 and a P-value cutoff of <0.05 were considered as DEmRNAs or DElncRNAs.
Pathway enrichment analysis of DEmRNAs
The Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) provides a comprehensive set of functional annotation tools for investigators to extract biological meaning from genes or proteins. The Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.kegg.jp/) is a database resource for functional interpretation and annotation of enriched pathways using large-scale gene datasets. The Gene Ontology (GO; http://www.geneontology.org/) database is commonly used to unify the representation of gene and gene product attributes. GO and KEGG pathway enrichment analyses were conducted for differentially expressed genes using DAVID. A P-value of <0.05 was selected as the cut-off criterion.
Protein-protein interaction (PPI) network construction
The Search Tool for the Retrieval of Interacting Genes (STRING; http://www.string-db.org/) database is an online database containing direct and indirect associations of proteins. DEmRNAs were mapped to STRING to assess the interactions. Next, PPI networks were established using Cytoscape software (version 3.6.0; http://www.cytoscape.org/). The Molecular Complex Detection (MCODE) plug-in was utilized to filter the modules of the PPI network in Cytoscape. An MCODE score of >5 was the selection criteria. Hub genes were exported. Function and pathway enrichment analyses were performed for DEmRNAs in the modules. A P-value of <0.05 was considered statistically significant.
Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)
RNA was extracted as aforementioned. RNA was quantified using a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Inc.). cDNA was generated using an RNeasy Mini kit (Takara Bio, Inc., Otsu, Japan). qPCR analysis was performed with SYBR-Green Master mix (Takara Bio, Inc.) for lncRNA and mRNA detection. The qPCR was performed at 95°C for 2 min, 40 cycles of 95°C for 5 sec and 60°C for 30 sec, one cycle of 95°C for 5 sec, 60°C for 1 min and 95°C for 15 sec, and finally, 50°C for 30 sec. Relative expression was analyzed using the 2−ΔΔCq method (12). Human GAPDH was used as an endogenous reference gene. The sequences of the primers are shown in Table I.
Construction of DEmRNA-DElncRNA co-expression network
A DElncRNA-DEmRNA co-expression network was constructed to explore the association between the RNAs. Pearson's correlation coefficient (PCC) was calculated between DElncRNAs and DEmRNAs according to corresponding expression levels. The criteria of |PCC| ≥0.90 was used to determine DElncRNA-DEmRNA pairs for network construction. The DElncRNA-DEmRNA network was constructed using Cytoscape 3.5.1 (http://cytoscape.org/).
Statistical analysis
The results of the qPCR were presented using Graph Pad Prism (version 6.0; GraphPad Software, Inc., La Jolla, CA, USA). Data are presented as the mean ± standard error of the mean. The independent samples t-test was performed for data comparison. P<0.05 was considered to indicate a statistically significant difference.
Results
Identification of DEmRNAs
According to microarray analysis, 32,206 mRNAs were detected in two pairs of cells, including 2,190 DEmRNAs (503 upregulated and 1,687 downregulated), which were significantly differentially expressed with a fold-change of ≥2.0. In addition, the DEmRNAs were divided into two sets: Set A contained all DEmRNAs from cell pair DLD1/A2 vs. DLD1/C6, and set B contained all DEmRNAs from DLD1/SA4 vs. DLD1/CA6. Heat map representation showed the top 100 differentially DEmRNAs (Fig. 1).
RT-qPCR validation of DEmRNAs
To validate the results of microarray analysis, RT-qPCR was performed in ATP5J and four randomly selected mRNAs. As shown in Fig. 2, cytotoxic and regulatory T cell molecule and cystatin SN were downregulated in DLD1/A2 vs. DLD1/C6 and upregulated in DLD1/SA4 vs. DLD1/CA6, while opposite results were obtained for ATP5J, cluster of differentiation 44 and eukaryotic translation elongation factor 1α2. The PCR results were basically concurrent with the microarray data.
Functional and pathway enrichment analysis
In order to investigate potential gene and gene products, GO functional and KEGG pathway enrichment analyses were performed for sets A and B, respectively. The DEmRNAs from set A were mainly enriched in biological processes associated with 'positive regulation of gene expression', 'cell-cell signaling' and 'regulation of cell growth'. The 10 most significantly enriched GO terms for DEmRNAs from sets A and B are shown in Fig. 3A and C, while the 10 most significantly enriched KEGG pathways in sets A and B are presented in Fig. 3B and D.
PPI network analysis
By submitting all DEmRNAs from sets A and B into STRING, the PPI network was obtained (data not shown). Next, the MCODE plug-in was used to analyze the results of STRING in Cytoscape software. The 3 most significant modules were acquired; the functions were mainly associated with 'cytokine-cytokine receptor interaction', 'ubiquitin-mediated proteolysis' and 'cell adhesion molecules (CAMs)' (Fig. 4). Moreover, a total of 20 genes were considered as hub genes (Table II), among which AKT serine/threonine kinase 2 (AKT2) obtained the highest value for neighborhood connectivity, outdegree and closeness centrality. This meant that AKT2 was the most important gene within the PPI network and represented the closest connection with other genes in the module (13).
Co-expression network of DEmRNA-DElncRNA
Based on microarray analysis, 39,304 lncRNAs were detected. A total of 3,002 DElncRNAs (2,319 in DLD1/A2 vs. DLD1/C6; 683 in DLD1/SA4 vs. DLD1/CA6) were significantly differentially expressed with a fold-change of ≥2.0. For further analysis, DElncRNAs upregulated in DLD1/A2 vs. DLD1/C6 and downregulated in DLD1/SA4 vs. DLD1/CA6 were selected, which meant that the expression of DElncRNAs was positively associated with the expression of ATP5J in the two cell pairs. The DElncRNAs that were negatively associated with ATP5J in the two cell pairs were also included. Certain DEmRNAs were selected according to the aforementioned criteria. Therefore, 51 DEmRNAs and 30 DElncRNAs were identified and listed in Tables III and IV.
To further investigate the potential association between DEmRNAs and DElncRNAs, PCC was calculated according to the expression levels of 51 DEmRNAs and 30 DElncRNAs. Using |PCC| ≥0.90, 343 DEmRNA/DElncRNA pairs were used to construct the co-expression network, including 239 pairs exhibiting positive correlation. The network is shown in Fig. 5. The pair of premature ovarian failure 1B (POF1B) and ENSG00000236654 presented the most significant positive correlation coefficiency, while the most marked negative correlation coefficiency was found in the pair of glycogen phosphorylase L and ENSG00000170011. One DElncRNA may target multiple DEmRNAs; for example, XIST was associated with POF1B, chorionic somatomammotropin hormone 1 (CSH1) and calmin (CLMN), leading to multiple signaling pathways or mechanisms.
RT-qPCR validation of hub genes and DElncRNA
RT-qPCR was used to validate AKT2 and XIST. As shown in Fig. 6, AKT2 and XIST were downregulated in DLD1/A2 vs. DLD1/C6 (P<0.01) and upregulated in DLD1/SA4 vs. DLD1/CA6 (P<0.05). In general, the validation results were consistent with the microarray data.
Discussion
ATP5J has been associated with several cancer types, including diffuse large B-cell lymphoma renal carcinoma and hepatocellular carcinoma (10,11,14). Furthermore, our previous study confirmed that ATP5J is associated with CRC cell migration and 5-fluorouracil resistance. The present study further investigated the molecular mechanism and potential key genes interacting with ATP5J in CRC via microarray and bioinformatics analyses.
The mRNAs from set A, which were differentially expressed in the upregulation of ATP5J, were mainly enriched in the GO terms of 'extracellular matrix organization', 'cell adhesion', 'positive regulation of gene expression' and 'regulation of cell growth'. These biological processes are essential for tumor cell survival and growth. Emerging data have confirmed the crucial role of the extracellular matrix in tumor metastasis (15,16). Changes in the content of the extracellular matrix may influence tumor cell properties, including proliferation and motility (17). Moreover, the primary enriched GO terms of set B contained enzyme regulation ('regulation of transcription from RNA polymerase II promoter') and biological responses ('immune response' and 'inflammatory response'). It is widely accepted that a causal association exists between inflammation, immunity and cancer (18). Human immunity responses, such as immune surveillance, are part of an innate and efficient antineoplastic system. Therefore, it could be concluded that the upregulation of ATP5J may affect cell migration and 5-FU sensitivity by influencing the biological processes of the intracellular and extracellular environment, thus promoting cell growth; its downregulation may operate by regulating immune and inflammatory responses. The KEGG pathways of set A were enriched primarily in tumor-related processes, including the phosphoinositide 3-kinase-protein kinase B signaling pathway, the tumor necrosis factor (TNF) signaling pathway and the Janus kinase-signal transducer and activator of transcription signaling pathway (19,20).
According to PPI network analysis, several biomarkers were identified. Some of these genes are associated with CRC biological behavior, including AKT2, TNF receptor superfamily member 19, claudin 7, cholecystokinin B receptor, glutathione peroxidase 3 and homer scaffolding protein 2. AKT2, one of the most significantly meaningful genes, has been shown to be essential in several cellular pathways involved in cell proliferation, metastasis and drug resistance (21). Overexpression of AKT2 has been detected in breast, ovarian and colon cancer (22,23). Previous studies have investigated the role of AKT2 in CRC via the construction of cell and animal models. The results showed that AKT2 deficiency may decrease the metastatic ability of CRC cells in mice (24). In addition, an shRNA-mediated AKT2-knockdown system has been shown to decrease cell survival and proliferation in primary tumors (25). Furthermore, low AKT2 expression in HCT116 cells has been shown to be correlated with increased chemosensitivity to paclitaxel (26). Taken together, the regulation of AKT2 may be a critical process for ATP5J exerting an effect on CRC cell biological functions.
Numerous studies have confirmed the association between lncRNAs and tumors. Previous studies have indicated that lncRNAs could have an important regulatory role in gene expression and are associated with cancer development (27). In the present study, microarray analysis was performed for lncRNAs. From the lncRNA expression profiles, DElncRNAs meeting the criteria were selected for further analysis. According to the results, TP73-AS1, HAS2-AS1, SFTA1P and XIST were highlighted and have been previously associated with the development and progression of certain cancer types, including hepatocellular carcinoma and oral squamous cell carcinoma (28–30). Among them, XIST was the only lncRNA that has been previously investigated in CRC (31). XIST, located in the X chromosome, was one of the first lncRNAs to be found that was involved in the progression of cancer, including non-small cell lung cancer and human glioblastoma (32,33). Certain studies have found the mechanism of XIST in several cancer types; for example, the lncRNA XIST has been shown to promote the malignancy of esophageal squamous cell carcinoma through the modulation of miR-101/enhancer of zeste homolog 2 (34). In addition, by regulating miR-186-5p, XIST controls cell proliferation and invasion in non-small cell lung cancer (35). Moreover, significantly high XIST expression has been found in CRC tumor tissues compared with that in paired adjacent normal tissues (36). High-level expression of XIST could predict poor disease-free survival times in CRC patients. Moreover, XIST is considered to be an oncogene that could promote cell proliferation through the miR-132-3p/mitogen-activated protein kinase 1 axis (31).
To further elucidate the potential functions of lncRNAs and find out their target genes in the cell pairs of the present study, a co-expression network of DEmRNA/DElncRNA was constructed. XIST was associated with 3 mRNAs (POF1B, CSH1 and CLMN). The gene encoding POF1B was critical for ovarian function (37). However, its impact was also confirmed in the regulation of cell adhesion in human intestinal cell lines (38). CLMN has been shown to regulate the cell cycle in a report involving neuroblastoma cells (39). Therefore, we hypothesized that the co-repression of XIST and POF1B and CLMN was involved in the molecular mechanisms of ATP5J in CRC cells, which will be confirmed by further cell line experiments in future studies.
In conclusion, the present study determined the expression profiles of mRNAs and lncRNAs in CRC cells overexpressing/underexpressing ATP5J. According to bioinformatics analyses, AKT2 and XIST expression was identified as a potential biomarker participating in the effect of ATP5J in CRC, which in turn was associated with cell migration and 5-FU resistance. Further research is necessary to confirm this hypothesis.