Identification of the key miRNA-mRNA regulatory network in lung adenocarcinoma
Original Article

Identification of the key miRNA-mRNA regulatory network in lung adenocarcinoma

Fangsu Xue1#, Wenjun You2#, Ying Liu3, Lou Zhong4

1Department of Respiration, Binhai County People’s Hospital, Yancheng, China; 2Department of Geriatrics, The Second People’s Hospital of Nantong, Affiliated Nantong Rehabilitation Hospital of Nantong University, Nantong, China; 3Department of Pathology, Affiliated Hospital of Nantong University, Nantong, China; 4Department of Thoracic Surgery, Affiliated Hospital of Nantong University, Nantong, China

Contributions: (I) Conception and design: F Xue, L Zhong; (II) Administrative support: L Zhong; (III) Provision of study materials or patients: F Xue, W You, Y Liu; (IV) Collection and assembly of data: F Xue, W You, Y Liu; (V) Data analysis and interpretation: L Zhong, F Xue; (VI) Manuscript writing: All authors; (VII) Final approval of the manuscript: All authors.

#These authors contributed equally to this work.

Correspondence to: Prof. Lou Zhong, MD. Department of Thoracic Surgery, Affiliated Hospital of Nantong University, Nantong 226001, China. Email: zhongl80@126.com.

Background: Lung adenocarcinoma (LUAD) is the most common type of lung cancer, has a high incidence rate and is a serious threat to human health. However, the pathogenesis of LUAD is still unclear. Further research on the pathogenesis of LUAD may provide targets for the early diagnosis and treatment of LUAD.

Methods: First, a transcriptome analysis was conducted to sequence the messenger RNA (mRNA) and micro RNA (miRNA) of the LUAD and adjacent control tissues. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were then performed for the functional annotation. Next, a differential miRNA-differential mRNA regulatory network was then constructed, and the function of the mRNAs in the network was analyzed and the key regulatory molecules (the hub molecules) were identified. Cytohubba was then used to analyze the top 20 hub molecules in the total miRNA-mRNA network, and the miRNAs regulating the top 20 hub genes (of which 2 were upregulated and 18 were downregulated). Finally, the key molecules were identified.

Results: By analyzing the function of the mRNA molecules in the regulatory network, we found that the immune response was inhibited, as were the movement and adhesion of immune-related cells; however, cell tumorigenesis, body death, and tumor cell proliferation were activated. The functions of the 20 hub molecules were mainly related to cytotoxicity, cell exosmosis, and cell adhesion mediated by immune cells. Further, we found that miR-5698, miR224-5p, and miR4709-3P regulate multiple key genes (e.g., PECAM1, CX3CR1, KLRD1, and CXCL12), and may be the key miRNAs regulating LUAD.

Conclusions: Immune response, cell tumorigenesis, and tumor cell proliferation play central roles in the overall regulatory network. miR-5698, miR224-5p, and miR4709-3p may be important biomarkers for the occurrence and development of LUAD and have great potential in the prognosis of LUAD patients and the development of new therapeutic targets.

Keywords: Lung adenocarcinoma (LUAD); micro RNA (miRNA); messenger RNA (mRNA)


Submitted Feb 15, 2023. Accepted for publication Apr 21, 2023. Published online Apr 27, 2023.

doi: 10.21037/jtd-23-400


Highlight box

Key findings

• This study found that immune response, cell tumorigenesis, and cell proliferation play central roles in the overall miRNA-mRNA regulatory network in lung adenocarcinoma (LUAD), and miR-5698, miR224-5p, and miR4709-3p may be important biomarkers for the occurrence and development of LUAD.

What is known and what is new?

• Targeting the key genes in the pathophysiological process of LUAD may be a potential way to treat LUAD.

• In this study, we constructed a microRNA-messenger RNA (miRNA-mRNA) regulatory network in LUAD and identified the key genes.

What is the implication, and what should change now?

• We constructed and analyzed a miRNA-mRNA regulatory network in LUAD, and identified new markers for evaluating the prognosis of LUAD. Our findings may inform the further study of the role of non-coding RNA in LUDA.


Introduction

Lung cancer has the highest incidence rate in the world and is also the most common cause of cancer-related death (1,2). At present, lung adenocarcinoma (LUAD) is the most common histological subtype of primary lung cancer (1). The relative incidence rate of LUAD continues to increase, and the prognosis of various histological subtypes is very poor (3-6).

Many treatment methods have been developed for lung cancer, including chemotherapy, radiotherapy, immunotherapy, and targeted treatment; however, it is still necessary to improve the overall survival rate through precision medicine (7,8). There is increasing evidence that malignant proliferation and distant metastasis are the main obstacles to the treatment of cancer (9-11). Extensive research has been conducted on the occurrence and development of LUAD; however, this research has failed to systematically study the molecular regulatory mechanism of the occurrence and development of LUAD overall, and consequently, the development of prognosis and treatment strategies for LUAD has been hindered. The targeting of the key genes in the pathophysiological process of LUAD may be a potential way to treat LUAD.

MicroRNA (miRNA) is a type of highly conserved tissue-specific small non-coding RNA. MiRNA mainly affects the gene expression level by targeting messenger RNA (mRNA). The relative levels of miRNA and mRNA play an important role in tumor genesis and development (12-14). A recent study has shown that miRNA and mRNA participate in the regulation of tumorigenesis and development. For example, in gastric cancer, miR-21-5p derived from exocrine body was shown to promote peritoneal metastasis through mesothelium-to-mesenchymal transformation, and thus may be a new therapeutic target for peritoneal metastasis of gastric cancer (15). In ovarian cancer, miR-424-5p was shown to regulate ferroptosis by targeting Acyl-CoA synthetase 4 (ACSL4) in ovarian cancer cells, and thus may be a target for inhibiting the progression of ovarian cancer (16). In hepatocellular carcinoma, differentially expressed miRNAs, such as miR-221, may provide new clues for the diagnosis and treatment of hepatocellular carcinoma patients based on their interacting with differentially expressed genes (DEGs), such as estrogen receptor-alpha (ESR1) and C-X-C motif chemokine 1 (CXCL12) (17). In neuroblastoma, miR-15a-5p, miR-15b-5p, and miR-16-5p have been shown to inhibit tumor progression by directly targeting v-myc avian myelocytomatosis viral oncogene neuroblastoma derived homolog gene (MYCN). These miRNAs are considered potential targets for neuroblastoma treatment (18). In bladder cancer, it was found that Histone lysine demethylase 6a (KDM6A) mRNA delivery via mucosal adhesion nanoparticles capsule inhibits the metastasis of bladder cancer (19). In pancreatic cancer, a comprehensive data analysis of mRNA, miRNA, and the signaling pathways revealed that denticleless E3 ubiquitin protein ligase homolog (DTL), cadherin 11 (CDH11), collagen type V alpha 1 chain (COL5A1), integrin subunit alpha 2 (ITGA2), kinesin family member 14 (KIF14), structural maintenance of chromosomes 4 (SMC4), versican (VCAN), miR-210, miR-217, miR-216a, miR-216b, miR-375, and miR-634 could be used as new diagnostic and even therapeutic markers in future research (20). In addition, by interacting with the DEGs CD44, ATP citrate lyase (ACLY), actin related protein 3 (ACTR3), signal transducer and activator of transcription 1 (STAT1), laminin subunit gamma 2 (LAMC2), and tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein zeta (YWHAZ), differentially expressed miRNAs (such as miR-5580-3p) may serve as suitable candidate biomarkers for the diagnosis, prognosis, and treatment of oral squamous cell carcinoma (21). However, most of previous studies focused on protein-coding genes or single signaling pathway in LUAD (22,23). Therefore, an extended understanding of the specific mechanisms of miRNA-mRNA regulatory network in LUAD is warranted.

In this study, we used transcriptome sequencing and a bioinformatics analysis to construct the miRNA-mRNA regulatory network in LUAD and identify the key genes. Our work extends understandings of the mechanism of miRNA and mRNA in LUAD. Their abnormal expression may be related to tumor proliferation, metastasis, and patient prognosis. Our findings may provide a new theoretical basis for the diagnosis, treatment, and prognosis of clinical malignant tumors. We present the following article in accordance with the STREGA reporting checklist (available at https://jtd.amegroups.com/article/view/10.21037/jtd-23-400/rc).


Methods

Samples

In this study, we collected the data of 3 primary LUAD patients undergoing surgical resections. The patients’ clinical information, such as age, gender, tumor node metastasis (TNM) classification, and stage, was basically consistent (Table 1). Cancer tissues and their adjacent tissues were collected for further sequencing analysis. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by the Institutional Ethics Committee of Binhai County People’s Hospital of Yancheng (No. 2020BYKYLL010) and informed consent was taken from all the patients.

Table 1

The patients’ clinical information

Patient Age (years) Gender Side T N M
No. 1 65 Male LL T1 N0 M0
No. 2 66 Male RL T1 N0 M0
No. 3 65 Male RL T1 N0 M0

LL, left lung; RL, right lung. T, tumor; N, node; M, metastasis.

Transcriptome sequencing

The mirVana™ miRNA isolation kit (Ambion, Austin, TX, USA) was used to extract the RNA. NanoDrop 2000 (Thermo Fisher Scientific Inc., Waltham, MA, USA) was used to quantify the total extracted RNA. In accordance with the manufacturer’s instructions, TruSeq Stranded Total RNA with Ribo Zero Gold (Cat. No. RS-122-2301, Illumina, CA, USA) was used to construct the mRNA library. The kit was used to digest the ribosomal RNA and break it into short fragments. A strand of complementary DNA (cDNA) was synthesized using the interrupted RNA as a template. To synthesize the other cDNA, 2'-deoxyuridine 5'-triphosphate (dUTP) was substituted with 2'-deoxythymidine 5'-triphosphate (dTTP), the different connectors were connected, 1 strand of cDNA containing dUTP was digested, and 1 strand of cDNA connecting the different connectors of the chain was retained. The purified cDNA strand was then repaired at the end, and a tail was added and connected to the sequencing connector, the fragment size was then selected, and finally, polymerase chain reaction (PCR) amplification was performed. The TruSeq Small RNA sample preparation kit (Cat. No. RS-200-0012, Illumina) was used to construct the Small RNA library. In short, the total RNA was connected to the connector at each end, and then the RNA connected to the connector was reverse transcribed into cDNA and amplified by PCR. The amplification products of 140–160 bp were isolated and purified into the miRNA library. The quality of the library was evaluated using the DNA high sensitivity chip of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Finally, the library was sequenced using the Illumina HiSeq X Ten platform (Illumina, CA, USA). The RNA sequencing and analysis were conducted by OE Biotech Co., Ltd. (Shanghai, China). The miRNA and mRNA expression profiles are available on NCBI SRA database, and the accession number is PRJNA932589.

Differential gene expression analysis

The raw reads obtained by the sequencing were filtered by the sortMeRNA (version v1.3.6, biocore) and Trimmomatic software (version 0.36, usadellab, German), and the FastQC (version 0.11.5, Babraham Bioinformatics, UK) was detected to obtain clean reads. HISAT2 (version 2.2.1.0, Johns Hopkins University, Center for Computational Biology, UK) was used to compare the clean reads, and RSeQC (version 2.6.4, Baylor College of Medicine, USA), suite was used to analyze the comparison results. The Stringtie software (version 1.3.3b, John Hopkins University, Baltimore, MD, USA) was used to assemble the reads matched to the gene to obtain a complete transcript. After the clean reads were compared to the transcription template by the Bowtie2 software (version 2.2.9, John Hopkins University, Baltimore, MD, USA), the transcripts were quantified by the eXpress software (version 1.5.1, Express Software Production, East Granby, CT, United States) to obtain the fragments per kilobase per million (FPKM) value and mRNA counts. The DESeq software (version 1.18.0, Fred Hutchinson Cancer Research Center, Seattle, WA, USA) was used to standardize the count of each sample mRNA and calculate the fold change (FC). A negative binomial distribution test was used to test the significance of the differences in the number of reads, and finally the differentially expressed mRNAs were screened using the following criteria: FC >2 or <–2, and a P value <0.05. The miRNAs were sequenced to obtain the raw reads. Cutadapt (version 1.14, National Bioinformatics Infrastructure, Uppsala, Sweden) was used to remove the connector sequence, and sequence lengths <15 and >41 bp were then removed. The fastx_Toolkit (version 0.0.13, Cancer Research UK Cambridge Institute, Cambridge, UK) and NGSQCToolkit (version 2.3.2, National Institute of Plant Genome Research, New Delhi, India) were used to obtain high-quality clean reads for the subsequent analysis. After using the Blastn software (version 2.2.26, National Center for Biotechnology, Bethesda, MD, USA) and RepeatMask (version 4.0.9, Institute for Systems Biology, USA) software to compare the clear reeds, the unannounced miRNA sequence was predicted and analyzed using the Mirdeep (version 2.0.0.8, Max Delbrueck Center for Molecular Medicine, Berlin, Germany) and RNAfold software (version 2.4.9, University of Vienna, Vienna, Austria), and miRNA expression was then calculated using the transcripts per million calculation metric based on the miRNA sequence, in which the total number of read matched pairs was used as the normalized expression standard. The Audic Claverie formula (Centre National de la Recherche Scientifique, Laboratory of Structural and Genetic Information, France) was used to calculate the p value, and miRNAs with a q value <0.05 and FC <–2 or >2 were then selected. Miranda (version 3.3a, http://www.microrna.org/microrna/, USA) as used to predict the target genes of the differentially expressed miRNAs. The following threshold parameters were used in the prediction of the miRNA target genes: S≥150, ΔG≤−30 kcal/mol and strict 5' seed pairing. The differential miRNA-differential mRNA regulatory network was thus constructed.

Gene function annotation

The mRNA transcript sequence was aligned to Swiss-Prot (http://www.gpmaw.com/html/swiss-prot.html, Swiss Institute of Bioinformatics & European Bioinformatics Institute, Swiss & UK), and the Gene Ontology (GO) annotation information of each mRNA transcript was obtained. The mRNA transcript sequence was compared to the KAAS database (http://www.genome.jp/tools/kaas/, Kyoto University Bioinformatics Center, Japan) to obtain the Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation information of each transcript sequence. A GO or KEGG group with a P value <0.05 was considered significantly enriched. Subsequently, Ingenuity Pathway Analysis (IPA, QIAGEN Redwood City, USA) software was used to analyze the molecular function of the DEGs, the methods of Z-score calculation can be found in causal analysis approaches in IPA (24), and the DEGs involved in the key molecular functions were screened for further analysis.

Network construction and analysis

STRING (version 11.5, http://string-db.org) was used to construct the protein-protein interaction (PPI) network. The key miRNA-mRNA subnetwork was constructed using Cytoscape (version 3.8.2, Institute for Systems Biology, Seattle, WA, USA). On the established PPI network, CytoHubba (version 0.1) and ClueGO (version 2.5.9) plugins of cystoscope were applied. The maximal clique centrality (MCC) algorithm was used to obtain the top genes and to build the core subnetworks of the different regulatory networks (25,26). ClueGO was used to analyze the main biological processes and KEGG pathways involved in the subnetwork.

Statistical analysis

When using RNA-sequencing data to compare and analyze whether the genes were differentially expressed in the two samples, the following two criteria were used: (I) the FC (i.e., the FC in the gene expression level in the two samples); and (II) the P value (mRNA) and q value (miRNA). The default filter criteria were a P value <0.05, and a default FC >2.


Results

Differential mRNA analysis of LUAD

RNA was extracted from the LUAD and paracancerous control tissues, and the whole transcriptome was sequenced. The principal component analysis (PCA) results showed that the gene expression of the LUAD samples differed significantly to that of the control group (Figure 1A). The volcanic map and thermal map results showed that 2,062 genes, of which 662 were upregulated and 1,400 were downregulated, were significantly changed in the LUAD tissue (Figure 1B,1C). The GO analysis of these 2,062 genes showed that the biological processes involved in these genes mainly included angiogenesis, the cell surface receiver signaling pathway, and the cell defense response (Figure 2A). The KEGG analysis showed that the pathways involved in these genes mainly included the cytokine-cytokine receptor interaction, the Ras-related protein 1 signaling pathway, and cell adhesion molecules (Figure 2B).

Figure 1 Changes in mRNA expression levels in lung adenocarcinoma. (A) PCA analysis of mRNAs; (B) volcano plot; (C) heat map of mRNAs. Differently expressed miRNAs (P<0.05; fold change <–2 or >2). The standardized expression level from blue to red represents the expression level from low to high; n=3. PC, principal component; LC, lung cancer; Nor, normal; FC, fold change; mRNA, messenger RNA; PCA, principal component analysis.
Figure 2 Differentially expressed mRNA function annotation. (A) GO analysis; the size and color of the circle indicate the number and −log10 (P value) of DEGs annotated in the specific GO biological processes, respectively; (B) KEGG analysis. BP, biological process; CC, cellular component; MF, molecular function; ABC, ATP binding cassette; HIF, hypoxia inducible factor; ECM, extracellular matrix; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differently expressed genes.

Differential miRNA analysis of LUAD

RNA was extracted from the LUAD and adjacent control tissues, and miRNA sequencing was performed. The PCA results showed that the gene expression of the LUAD samples differed significantly to that of the control group (Figure 3A). The volcanic and thermographic results showed that there were significant changes in the expression of the genes, of which 50 miRNAs were upregulated and 44 miRNAs were downregulated, in the LUAD tissues (Figures 3B,3C). The GO analysis results showed that the biological processes involved in these differentially expressed miRNA target genes mainly included microtubule anchoring, the regulation of microtubule based process, the positive regulation of T cell activation, and the positive regulation of extracellular matrix assembly (Figure 4A). The KEGG analysis showed that the pathways involved in these genes mainly included proteoglycans in cancer, the relaxin signaling pathway, and the phospholipase D signaling pathway (Figure 4B).

Figure 3 Changes in miRNA expression levels in lung adenocarcinoma. (A) PCA analysis of miRNA; (B) volcano plot; (C) heat map of miRNA. Differently expressed miRNAs (q<0.05; fold change <–2 or >2). The standardized expression level from blue to red represents the expression level from low to high; n=3. PC, principal component; LC, lung cancer; Nor, normal; FC, fold change; miRNA, microRNA; PCA, principal component analysis.
Figure 4 Functional annotation of the differentially expressed miRNA target genes. (A) GO analysis; the size and color of the circle indicate the number and −log10 (P value) of DEGs annotated in the specific GO biological processes, respectively; (B) KEGG analysis. The color of the bar indicated −log10 (P value) of DEGs annotated in the specific pathway. MHC, major histocompatibility complex; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix; GnRH, gonadotropin-releasing hormone; TRP, transient receptor potential; GO, Gene Ontology; DEGs, differently expressed genes.

Construction of the miRNA-mRNA regulatory network

Subsequently, the differential miRNA-differential mRNA regulatory network was constructed. The results showed that 40 upregulated miRNAs and 780 downregulated mRNAs constituted one miRNA-mRNA regulatory network, and 38 downregulated miRNAs and 412 upregulated mRNAs constituted the other miRNA-mRNA regulatory network. The most significant differences in the expression of the miRNAs in the regulatory network are shown in Table 2. The mRNAs in the regulatory network were further analyzed, and the results showed that they mainly participated in the following biological processes: the cell surface receiver signaling pathway, angiogenesis, the G-protein coupled receiver signaling pathway, and the regulation of immune response. The mRNAs mainly participated in the following signaling pathways: cell adhesion molecules, the cyclic adenosine monophosphate signaling pathway, the cyclic guanosine monophosphate-protein kinase G (cGMP-PKG) signaling pathway, and natural killer cell-mediated cytotoxicity (Figure 5).

Table 2

The most significant differentially expressed miRNA in the regulatory network

miRNA Log2 fold change P value
hsa-miR-659-5p 5.999457764 0.001341028
hsa-miR-5698 5.218436314 0.000864982
hsa-miR-4697-3p 4.964021707 0.0000239
hsa-miR-877-3p 4.773443991 0.002480896
hsa-miR-3617-5p 4.120845301 0.002295096
hsa-miR-6508-3p 4.063918642 0.002329101
hsa-miR-135b-3p 4.019552756 3.76E-10
hsa-miR-323b-3p 3.72081132 0.0000875
hsa-miR-551a 3.482910763 0.002455497
hsa-miR-135b-5p 3.442293207 3.72E-14
hsa-miR-148a-5p 3.294543791 0.000215177
hsa-miR-21-3p 3.255318989 0.00000561
hsa-miR-6772-3p 3.183203387 0.004553841
hsa-miR-182-3p 3.079397721 0.000000021
hsa-miR-4645-3p 3.022222859 0.002557073
hsa-miR-556-3p 3.020623129 0.002026974
hsa-miR-323a-3p 2.889515303 0.002841095
hsa-miR-5091 2.841382808 0.00055338
hsa-miR-224-5p 2.617202452 0.00000688
hsa-miR-9-5p 2.60547639 0.000715377
hsa-miR-584-3p −4.723274554 0.003618141
hsa-miR-486-5p −4.404575547 3.75E-18
hsa-miR-486-3p −4.324459857 8.09E-31
hsa-miR-4732-5p −4.267187316 0.00000133
hsa-miR-144-5p −4.030408818 7.43E-09
hsa-miR-942-5p −3.844620538 0.000000647
hsa-miR-451a −3.832194137 1.43E-09
hsa-miR-144-3p −3.636199209 0.00000301
hsa-miR-451b −3.363008176 0.000235816
hsa-miR-363-3p −3.292553099 0.000000102
hsa-miR-584-5p −3.15737839 7.15E-08
hsa-miR-20b-5p −3.056183711 0.000000396
hsa-miR-106a-5p −2.896005705 0.000000497
hsa-miR-3158-3p −2.71250167 0.00000378
hsa-miR-196b-5p −2.606892503 0.0000014
hsa-miR-4504 −2.499091633 0.00129907
hsa-miR-139-3p −2.359094898 0.000000307
hsa-miR-126-3p −2.293041988 0.00000275
hsa-miR-204-5p −2.121067315 0.000736143
hsa-miR-15a-5p −2.065622117 0.0000467

miRNA, microRNA.

Figure 5 mRNA function annotation in the miRNA-mRNA regulatory network. (A) GO analysis; the size and color of the circle indicate the number and −log10 (P value) of DEGs annotated in the specific GO biological processes, respectively; (B) KEGG analysis. The color of the bar indicated −log10 (P value) of DEGs annotated in the specific pathway. BP, biological process; CC, cellular component; MF, molecular Function; ECM, extracellular matrix; KEGG, Kyoto Encyclopedia of Genes and Genomes; ABC, ATP binding cassette; cAMP, Cyclic Adenosine 3',5'-monophosphate; cGMP, cyclic guanosine monophosphate; PKG, Protein kinase G; MAPK, Mitogen-activated protein kinases; miRNA-mRNA, microRNA-messenger RNA; GO, Gene Ontology; DEGs, differently expressed genes.

Functional analysis and key molecular analysis of the differentially expressed mRNAs in the regulatory network

The IPA software was used to analyze the function of the mRNA molecules in the regulatory network and the results showed that the immune response was inhibited, and the movement and adhesion of the immune-related cells were also inhibited, while the neoplasia of cells, organismal death, and tumor cell proliferation was activated (Figure 6). In total, 237 differentially expressed genes were found to be involved in proliferation regulation, of which 151 were downregulated, such as peptidyl arginine deiminase 4 (Padi4), and 86 were upregulated, such as cartilage oligomeric matrix protein (Comp). Through the construction of the protein-protein interaction network and the Cytohubba analysis, the following 10 key regulatory molecules were identified: epidermal growth factor receptor (EGFR), platelet and endothelial cell adhesion molecule 1 (PECAM1), CXCL12, interleukin 1 beta (IL1β), CDH1, selectin E (SELE), fms related receptor tyrosine kinase 1 (FLT1), peroxisome proliferator activated receptor gamma (PPARG), CD24, and epithelial cell adhesion molecule (EPCAM) (Figure 7). In total, 650 differentially expressed genes were found to be involved in cell tumorigenesis, of which 206 were upregulated, such as solute carrier family 25 member 48 (SLC25A48), and 444 were downregulated, such as INSC spindle orientation adaptor protein (INSC). In the interaction network, the following 10 key regulatory molecules were identified: protein tyrosine phosphatase receptor type C (PTPRC), CD34, PECAM1, IL1β, C-C motif chemokine ligand 5 (CCL5), CXCL12, SELE, CDH1, integrin subunit alpha X (ITGAX), and PPARG (Figure 8). Additionally, 129 differentially expressed genes were found to be involved in immune response, of which 107 were downregulated, such as Wnt3a, and 22 were upregulated, such as Wuc5b. The analysis also identified the following 10 key regulatory molecules: IL1β, toll like receptor 4 (TLR4), PECAM1, CCL5, CXCL12, SELE, C-X-C motif chemokine receptor 2 (CXCR2), EGFR, CXCR1, and selectin P ligand (SELPLG) (Figure 9).

Figure 6 IPA analysis of mRNA molecular function in the miRNA-mRNA regulatory network. The bars indicate z score, with >2 for activation and <−2 for inhibition. The black circle in line indicate statistical significance [−log10 (P value)] for each term. IPA, Ingenuity Pathway Analysis; miRNA-mRNA, microRNA-messenger RNA.
Figure 7 Differentially expressed genes related to proliferation regulation in the miRNA-mRNA regulatory network. (A) Differentially expressed genes; (B) Hub molecules. The scores are shown in red color. A darker color means a higher score. miRNA-mRNA, microRNA-messenger RNA.
Figure 8 Differentially expressed genes of cell tumorigenesis in the miRNA-mRNA regulatory network. (A) Differentially expressed genes; (B) Hub molecules. The scores are shown in red color. A darker color means a higher score. miRNA-mRNA, microRNA-messenger RNA.
Figure 9 Differentially expressed genes related to the inflammatory response in the miRNA-mRNA regulatory network. (A) Differentially expressed genes; (B) Hub molecules. The scores are shown in red color. A darker color means a higher score. miRNA-mRNA, microRNA-messenger RNA.

Cytohubba was used to analyze the top 20 hub molecules in the total miRNA-mRNA network (Figure 10). These molecules were highly coincident with the hub molecules in the above three processes, indicating that immune response, cell tumorigenesis, and tumor cell proliferation play central roles in the overall regulatory network. CluoGO was used to analyze these 20 core genes, which were found to be related to immune cell-mediated cytotoxicity, cell exosmosis, and cell adhesion.

Figure 10 Key differentially expressed genes in the miRNA-mRNA regulatory network. (A) Hub molecules. The scores are shown in red color. A darker color means a higher score; (B) function analysis. miRNA-mRNA, microRNA-messenger RNA.

Analysis of key miRNAs in regulatory networks

We then analyzed the miRNAs that regulate the above-mentioned 20 hub genes (of which 2 were upregulated and 18 were downregulated), and found that miR-5698, miR224-5p, and miR4709-3p regulate CD34, CX3CR1, and CXCL12 (Figure 11). These 3 miRNAs are presumed to be the key miRNAs that regulate LUAD.

Figure 11 Key differentially expressed miRNAs regulating the hub genes. The scores are shown in red color (miRNA) and blue color (mRNA). A darker color means a higher fold change. miRNA, microRNA.

Discussion

LUAD is a malignant tumor, has shown the fastest increase in incidence and mortality rates, and thus represents a great threat to human health and life (10). At present, its pathogenesis is complex and unclear. Some studies have shown that mRNA and miRNA play important roles in the development of LUAD (27,28). A bioinformatics analysis showed that 4 differentially expressed miRNAs (i.e., miR-196b-5p, miR-223-3p, miR-135a-5p, and miR-9) and 4 differentially expressed mRNAs (i.e., MEIS1, HOXA5, SCARB1, and RUNX2) participate in the regulation network of the molecular mechanism of bone metastasis in LUAD (29). CircRNA CRIM1 inhibits the invasion and metastasis of LUAD through the miR-182/miR-93 leukemia inhibitory factor receptor pathway (30). In this study, we extended understandings of the molecular mechanism of LUAD by constructing a miRNA-mRNA regulatory network to analyze the key genes involved in the development of LUAD.

According to the regulatory network, we found that miR-5698, miR-224-5p, and miR-4709-3p regulate multiple key hub genes identified by using cytoHubba’s MCC algorithm (31,32). These miRNAs are presumed to be the key miRNAs for regulating LUAD. These miRNAs have also been reported to be involved in the occurrence and development of other common cancers. By establishing a prediction model based on serum miRNA, researchers found that the serum level of miR-5698 was significantly correlated with the total survival time after the beginning of eribulin treatment, and thus could be used as a biomarker of patients’ responses to eribulin and to predict new distant metastasis in metastatic breast cancer (33). Research has also shown that miR-224-5p delivered by exosomes from cancer-associated fibroblasts enter clear cell renal cell carcinoma cells, and then promote the malignant behavior of clear cell renal cell carcinoma cells, which indicates that miR-224-5p has potential value as a therapeutic target for clear cell renal cell carcinoma (34). Further research has shown that miR-4709-3p inhibits the tumorigenesis of breast cancer cells by targeting testes-specific protease 50. A previous study showed that nuclear factor kappa beta and activin signal transduction are involved in miR-4709-3p-related tumor inhibition (35). In this study, miR-5698, miR-224-5p, and miR-4709-3p were considered important biomarkers of LUAD, but the specific pathogenesis of LUAD requires further study.

In this study, we also found that miR-5698 regulates CD34, ITGAL, PTPRC, PRF1, CD247, PECAM1, CX3CR1, KLRD1, CXCL12, and TBX21, miR-224-5p regulates FCGR3B, TBX21, ITGAX, PTPRC, PECAM1, CX3CR1, KLRD1, and CXCL12, and miR-4709-3p regulates CD247, NCR1, IL1B, PECAM1, CX3CR1, KLRD1, CXCL12, FASLG, CD34, and ITGAX. The common genes regulated by miR-5698, miR-224-5p, and miR-4709-3p are PECAM1, CX3CR1, KLRD1, and CXCL12. These genes may also be important therapeutic targets for LUAD. There is evidence that PECAM1, CX3CR1, KLRD1, and CXCL12 are involved in the progression of multiple malignant tumors. For example, in non-small cell lung cancer, a survival analysis showed that the high expression of PECAM1 is correlated with an improved survival rate, which suggests that PECAM1 may be a potential prognostic marker of lung cancer (36). In colorectal cancer, CX3CR1 was shown to contribute to the recruitment and regulation of immune infiltrating cells and the polarization of macrophages in the progression of colorectal cancer, which suggests that CX3CR1 can be used as a prognostic biomarker for colorectal cancer (37). In a univariate Cox analysis of 43 mRNAs, KLRD1 was identified as a prognostication-related gene (P<0.05) in chromophobe renal cell carcinoma (38). An analysis of the comprehensive characterization of the microenvironment of prostate cancer showed that the CXCR4/CXCL12 interaction could serve as a potential new target to interfere with the angiogenesis of prostate cancer (39). In this study, we found that the four common genes functions mainly involve immune response, cell tumorigenesis, and tumor cell proliferation, which is consistent with the findings of the above-mentioned studies. This study had some limitations. First, Timer and Cibersort methods could be used to detect the differences in the proportion of immune cells between samples to analyze the immune infiltration pattern of LUAD. Secondly, we performed transcriptome sequencing but did not carry out a biofunctional verification of the sequencing results. We aim to fill this knowledge gap in future studies.

In conclusion, in this study, we tried to identify the expression, potential function, and clinical value of miRNAs and mRNAs in LUAD. Our findings may contribute to the personalized prediction of the prognosis of LUAD patients, and as potential biomarkers, miRNAs and mRNAs may reflect the response of LUAD patients to the specific targeted treatment of LUAD markers. This study highlighted the important role of miRNAs and mRNAs in the occurrence and development of LUAD, and identified potential guiding biomarkers for the selection of treatment methods.


Conclusions

In this study, we constructed and analyzed a miRNA-mRNA regulatory network in LUAD through bioinformatics and screened out the key genes. We conducted a functional analysis and found that immune response, cell tumorigenesis and tumor cell proliferation play central roles in the overall regulatory network. In addition, we speculated that miR-5698, miR224-5p, and miR4709-3p may be important biomarkers for the occurrence and development of LUAD and have great potential in the prognosis and treatment of LUAD patients. Our study identified new markers for evaluating the prognosis of LUAD patients and provides important evidence for further research on the role of non-coding RNA in LUAD.


Acknowledgments

Funding: This work was supported by the Natural Science Research Project of Nantong Science and Technology Bureau (No. MS12020029).


Footnote

Reporting Checklist: The authors have completed the STREGA reporting checklist. Available at https://jtd.amegroups.com/article/view/10.21037/jtd-23-400/rc

Data Sharing Statement: Available at https://jtd.amegroups.com/article/view/10.21037/jtd-23-400/dss

Peer Review File: Available at https://jtd.amegroups.com/article/view/10.21037/jtd-23-400/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jtd.amegroups.com/article/view/10.21037/jtd-23-400/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by the Institutional Ethics Committee of Binhai County People’s Hospital of Yancheng (No. 2020BYKYLL010) and informed consent was taken from all the patients.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.


References

  1. Bray F, Ferlay J, Soerjomataram I, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018;68:394-424. [Crossref] [PubMed]
  2. Wu Y, Jiang M. The revolution of lung cancer treatment: from vaccines, to immune checkpoint inhibitors, to chimeric antigen receptor T therapy. Biotarget 2017;1:7. [Crossref]
  3. Meza R, Meernik C, Jeon J, et al. Lung cancer incidence trends by gender, race and histology in the United States, 1973-2010. PLoS One 2015;10:e0121323. [Crossref] [PubMed]
  4. Lortet-Tieulent J, Soerjomataram I, Ferlay J, et al. International trends in lung cancer incidence by histological subtype: adenocarcinoma stabilizing in men but still increasing in women. Lung Cancer 2014;84:13-22. [Crossref] [PubMed]
  5. Travis WD, Brambilla E, Nicholson AG, et al. The 2015 World Health Organization Classification of Lung Tumors: Impact of Genetic, Clinical and Radiologic Advances Since the 2004 Classification. J Thorac Oncol 2015;10:1243-60. [Crossref] [PubMed]
  6. Noma D, Inamura K, Matsuura Y, et al. ALK-rearranged lung adenocarcinoma showing intra-bronchial protrusion: a case of actually peripheral origin with a rare spreading pattern. Biotarget 2017;1:15. [Crossref]
  7. Zhang L, Li M, Wang W, et al. Celecoxib alleviates denervation-induced muscle atrophy by suppressing inflammation and oxidative stress and improving microcirculation. Biochem Pharmacol 2022;203:115186. [Crossref] [PubMed]
  8. Wu L, Cheng D, Yang X, et al. M2-TAMs promote immunoresistance in lung adenocarcinoma by enhancing METTL3-mediated m6A methylation. Ann Transl Med 2022;10:1380. [Crossref] [PubMed]
  9. Chaffer CL, Weinberg RA. A perspective on cancer cell metastasis. Science 2011;331:1559-64. [Crossref] [PubMed]
  10. Zhao Y, Brasier AR. Uronic acid pathway metabolites regulate mesenchymal transition and invasiveness in lung adenocarcinoma. Biotarget 2019;3:19. [Crossref] [PubMed]
  11. Ohtsuka M, Tanemura M, Akamatsu H. Long noncoding RNAs regulate malignant phenotypes in colorectal cancer. Biotarget 2018;2:4. [Crossref]
  12. Hill M, Tran N. miRNA interplay: mechanisms and consequences in cancer. Dis Model Mech 2021;14:dmm047662. [Crossref] [PubMed]
  13. Mishra S, Yadav T, Rani V. Exploring miRNA based approaches in cancer diagnostics and therapeutics. Crit Rev Oncol Hematol 2016;98:12-23. [Crossref] [PubMed]
  14. Scheunemann D, Pradhan AK, Das SK, et al. Wnt7a and miR-370-3p: new contributors to bladder cancer invasion. Biotarget 2018;2:14. [Crossref] [PubMed]
  15. Li Q, Li B, Li Q, et al. Exosomal miR-21-5p derived from gastric cancer promotes peritoneal metastasis via mesothelial-to-mesenchymal transition. Cell Death Dis 2018;9:854. [Crossref] [PubMed]
  16. Ma LL, Liang L, Zhou D, et al. Tumor suppressor miR-424-5p abrogates ferroptosis in ovarian cancer through targeting ACSL4. Neoplasma 2021;68:165-73. [Crossref] [PubMed]
  17. Mou T, Zhu D, Wei X, et al. Identification and interaction analysis of key genes and microRNAs in hepatocellular carcinoma by bioinformatics analysis. World J Surg Oncol 2017;15:63. [Crossref] [PubMed]
  18. Chava S, Reynolds CP, Pathania AS, et al. miR-15a-5p, miR-15b-5p, and miR-16-5p inhibit tumor progression by directly targeting MYCN in neuroblastoma. Mol Oncol 2020;14:180-96. [Crossref] [PubMed]
  19. Kong N, Zhang R, Wu G, et al. Intravesical delivery of KDM6A-mRNA via mucoadhesive nanoparticles inhibits the metastasis of bladder cancer. Proc Natl Acad Sci U S A 2022;119:e2112696119. [Crossref] [PubMed]
  20. Sohrabi E, Rezaie E, Heiat M, et al. An Integrated Data Analysis of mRNA, miRNA and Signaling Pathways in Pancreatic Cancer. Biochem Genet 2021;59:1326-58. [Crossref] [PubMed]
  21. Amiri-Dashatan N, Koushki M, Jalilian A, et al. Integrated Bioinformatics Analysis of mRNAs and miRNAs Identified Potential Biomarkers of Oral Squamous Cell Carcinoma. Asian Pac J Cancer Prev 2020;21:1841-8. [Crossref] [PubMed]
  22. Wu C, Dong B, Huang L, et al. SPTBN2, a New Biomarker of Lung Adenocarcinoma. Front Oncol 2021;11:754290. [Crossref] [PubMed]
  23. Xu JY, Zhang C, Wang X, et al. Integrative Proteomic Characterization of Human Lung Adenocarcinoma. Cell 2020;182:245-261.e17. [Crossref] [PubMed]
  24. Krämer A, Green J, Pollard J Jr, et al. Causal analysis approaches in Ingenuity Pathway Analysis. Bioinformatics 2014;30:523-30. [Crossref] [PubMed]
  25. Teng L, Shen L, Zhao W, et al. SLAMF8 Participates in Acute Renal Transplant Rejection via TLR4 Pathway on Pro-Inflammatory Macrophages. Front Immunol 2022;13:846695. [Crossref] [PubMed]
  26. Chin CH, Chen SH, Wu HH, et al. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol 2014;8:S11. [Crossref] [PubMed]
  27. Cao X, Xue F, Chen H, et al. MiR-202-3p inhibits the proliferation and metastasis of lung adenocarcinoma cells by targeting RRM2. Ann Transl Med 2022;10:1374. [Crossref] [PubMed]
  28. Zhang L, Liang J, Han Z, et al. Micro-ribonucleic acids (miRNAs) and a proteomic profile in lung adenocarcinoma cases with brain metastasis. Ann Transl Med 2022;10:1389. [Crossref] [PubMed]
  29. Han L, Yao Z, Xie L, et al. Transcriptome Sequencing reveals the expressed profiles of mRNA and ncRNAs and regulate network via ceRNA mediated molecular mechanism of lung adenocarcinoma bone metastasis in Xuanwei. Transl Cancer Res 2021;10:73-87. [Crossref] [PubMed]
  30. Wang L, Liang Y, Mao Q, et al. Circular RNA circCRIM1 inhibits invasion and metastasis in lung adenocarcinoma through the microRNA (miR)-182/miR-93-leukemia inhibitory factor receptor pathway. Cancer Sci 2019;110:2960-72. [Crossref] [PubMed]
  31. Sun Y, Cai D, Hu W, et al. Identifying hub genes and miRNAs in Crohn's disease by bioinformatics analysis. Front Genet 2022;13:950136. [Crossref] [PubMed]
  32. Yang Q, Li K, Li X, et al. Identification of Key Genes and Pathways in Myeloma side population cells by Bioinformatics Analysis. Int J Med Sci 2020;17:2063-76. [Crossref] [PubMed]
  33. Satomi-Tsushita N, Shimomura A, Matsuzaki J, et al. Serum microRNA-based prediction of responsiveness to eribulin in metastatic breast cancer. PLoS One 2019;14:e0222024. [Crossref] [PubMed]
  34. Sun H, Sun J, Li M, et al. Transcriptome Analysis of Immune Receptor Activation and Energy Metabolism Reduction as the Underlying Mechanisms in Interleukin-6-Induced Skeletal Muscle Atrophy. Front Immunol 2021;12:730070. [Crossref] [PubMed]
  35. Chen Z, Chen X, Ji Y, et al. A narrative review of the role of m6A in oxidative stress and inflammation. Biotarget 2021;5:1. [Crossref]
  36. Cao S, Wang Y, Li J, et al. Prognostic Implication of the Expression Level of PECAM-1 in Non-small Cell Lung Cancer. Front Oncol 2021;11:587744. [Crossref] [PubMed]
  37. Yue Y, Zhang Q, Sun Z. CX3CR1 Acts as a Protective Biomarker in the Tumor Microenvironment of Colorectal Cancer. Front Immunol 2021;12:758040. [Crossref] [PubMed]
  38. Borkamo ED, Dahl O, Bruland O, et al. Kinetics study on markers of the immune system by gene expression profiling of an in vivo heated tumor. Int J Hyperthermia 2009;25:41-6. [Crossref] [PubMed]
  39. Heidegger I, Fotakis G, Offermann A, et al. Comprehensive characterization of the prostate tumor microenvironment identifies CXCR4/CXCL12 crosstalk as a novel antiangiogenic therapeutic target in prostate cancer. Mol Cancer 2022;21:132. [Crossref] [PubMed]

(English Language Editor: L. Huleatt)

Cite this article as: Xue F, You W, Liu Y, Zhong L. Identification of the key miRNA-mRNA regulatory network in lung adenocarcinoma. J Thorac Dis 2023;15(4):2037-2050. doi: 10.21037/jtd-23-400

Download Citation