Skip to main content

Comprehensive identification of RNA transcripts and construction of RNA network in chronic obstructive pulmonary disease

Abstract

Background

Chronic obstructive pulmonary disease (COPD) is one of the world’s leading causes of death and a major chronic disease, highly prevalent in the aging population exposed to tobacco smoke and airborne pollutants, which calls for early and useful biomolecular predictors. Roles of noncoding RNAs in COPD have been proposed, however, not many studies have systematically investigated the crosstalk among various transcripts in this context. The construction of RNA functional networks such as lncRNA-mRNA, and circRNA-miRNA-mRNA interaction networks could therefore facilitate our understanding of RNA interactions in COPD. Here, we identified the expression of RNA transcripts in RNA sequencing from COPD patients, and the potential RNA networks were further constructed.

Methods

All fresh peripheral blood samples of three patients with COPD and three non-COPD patients were collected and examined for mRNA, miRNA, lncRNA, and circRNA expression followed by qRT-PCR validation. We also examined mRNA expression to enrich relevant biological pathways. lncRNA-mRNA coexpression network and circRNA-miRNA-mRNA network in COPD were constructed.

Results

In this study, we have comprehensively identified and analyzed the differentially expressed mRNAs, lncRNAs, miRNAs, and circRNAs in peripheral blood of COPD patients with high-throughput RNA sequencing. 282 mRNAs, 146 lncRNAs, 85 miRNAs, and 81 circRNAs were differentially expressed. GSEA analysis showed that these differentially expressed RNAs correlate with several critical biological processes such as “ncRNA metabolic process”, “ncRNA processing”, “ribosome biogenesis”, “rRNAs metabolic process”, “tRNA metabolic process” and “tRNA processing”, which might be participating in the progression of COPD. RT-qPCR with more clinical COPD samples was used for the validation of some differentially expressed RNAs, and the results were in high accordance with the RNA sequencing. Given the putative regulatory function of lncRNAs and circRNAs, we have constructed the co-expression network between lncRNA and mRNA. To demonstrate the potential interactions between circRNAs and miRNAs, we have also constructed a competing endogenous RNA (ceRNA) network of differential expression circRNA-miRNA-mRNA in COPD.

Conclusions

In this study, we have identified and analyzed the differentially expressed mRNAs, lncRNAs, miRNAs, and circRNAs, providing a systematic view of the differentially expressed RNA in the context of COPD. We have also constructed the lncRNA-mRNA co-expression network, and for the first time constructed the circRNA-miRNA-mRNA in COPD. This study reveals the RNA involvement and potential regulatory roles in COPD, and further uncovers the interactions among those RNAs, which will assist the pathological investigations of COPD and shed light on therapeutic targets exploration for COPD.

Introduction

Chronic obstructive pulmonary disease (COPD) is one of the world’s leading causes of death characterized by severe respiratory symptoms and airflow limitation [1]. With the increasing smoking rates in developing countries and the aging population in high-income countries, World Health Organization predicted that the prevalence rate of COPD will keep rising in the next 40 years and that the number of patients who die of COPD and related diseases will exceed 5.4 million per annum by 2060 [2]. At present, treatment of COPD mainly focuses on controlling symptoms and limiting disease progression, and the development of targeted therapeutic drugs is slow [3]. Genetic basis of COPD remains elucidated, with specific and effective biomolecules for COPD early clinical prevention and treatments are still lacking [4].

Dysregulation of gene expression is a universal theme of human diseases, and the resulting abnormal mRNAs and non-coding RNAs play important roles in the occurrence and progression of various diseases, including COPD [5]. Next-generation sequencing facilitates our understanding of differentially expressed RNA and its significance in COPD. For example, 82 gene locus were found to be associated with COPD through a genome-wide association study (GWAS), suggesting a genetic basis for the clinical heterogeneity seen in COPD [6]. Peripheral blood could also be used to capture relevant lung pathobiology through RNA-seq profiling [7]. A series of COPD-related genes have also been identified in gene expression profiling, with more comprehensive assessments of rare genetic variations are increasingly being used in COPD genetic research [8].

MiRNAs are critical regulators involved in multiple physiological and pathological processes, including development, organogenesis, apoptosis, and cell proliferation and differentiation [9]. Dysregulation of miRNA often leads to impaired cellular function and disease [10]. MicroRNA-218 acts by inhibiting TNFR1-mediated activation of NF-κB, which is associated with the overproduction of MUC5AC and inflammation in smoking-induced bronchiolitis of COPD [11]. In terms of treatment, through bronchial biopsies of 69 moderate/severe COPD patients and related in vitro experiments, miR-320d inhibits the expression of inflammatory cytokines and regulates the pro-inflammatory response of the airway epithelium in patients with COPD as a novel mediator of Inhaled Corticosteroids (ICS) [12].

LncRNAs are widely involved in various biological processes, such as regulation of transcription, RNA splicing, etc. [13]. Aberrant expressions of lncRNA are closely related to various human diseases. For instance, SIRT1/p53 and FOXO3a signaling pathways mediated by senescence-associated lncRNA1 (SAL-RNA1) regulate the expression of SIRT1 and FOXO3a, promoting the senescence of type II alveolar epithelial cells (AECII) in the pathogenesis of COPD [14]. In a microarray study of genome-wide lncRNA expression profiling in lung tissue of patients with COPD, 120 lncRNAs were up-regulated and 43 lncRNAs were down-regulated in the lung tissue of COPD patients compared with normal smokers without COPD [15].

CircRNAs are a special class of non-coding RNAs with a covalently-closed loop generally through 5’ to 3’ ends back splicing in animals, involved in various diseases including COPD. circRNAs play a significant role in regulating gene expression by binding to specific microRNAs or RNA binding proteins. For instance, hsa_Circ_0016070 promotes the proliferation of pulmonary artery smooth muscle cells (PASMCs) and participates in vascular remodeling in pulmonary hypertension (PAH) via the miR-942/CCND1 axis [16]. Epithelial–mesenchymal transition (EMT) may promote small airway remodeling and fibrosis in patients with COPD and increase the likelihood of lung cancer [17], circ0061052 regulates the expression of FoxC1/Snail by competitively binding miR-515-5p in human bronchial epithelial cells and involves in the CS-induced EMT and airway remodeling in COPD [18].

In this study, we have comprehensively identified differentially expressed mRNAs, miRNAs, lncRNAs, and circRNAs in the peripheral blood of COPD patients through high-throughput RNA sequencing. GSEA analysis showed that these differentially expressed RNAs correlate with several RNA biological processes such as ncRNA metabolic process, ncRNA processing, ribosome biogenesis, rRNA metabolic process, tRNA metabolic process, and tRNA processing. To demonstrate the potential interactions between lncRNAs and their interacting mRNAs, we have constructed the lncRNA-mRNA co-expression network in COPD. For the validated circRNAs, we have analyzed the features of the circRNAs and predicted their potential functional mechanisms in COPD. We further for the first time constructed the circRNA-miRNA-mRNA network in COPD, which provides a valuable prognostic or predictive resource for future COPD research and clinical treatments.

Methods

Clinical samples

All fresh peripheral blood samples of three patients with COPD and three non-COPD patients from the Second Affiliated Hospital of Anhui Medical University, which was approved by the Ethics Committee of the Second Affiliated Hospital of Anhui Medical University [YX2021-148(F1)]. Written informed consent was obtained from each patient for this study. Information of the patients is included in Tables 1 and 2.

Table 1 Information of the patients for RNA-seq
Table 2 Information of the patients for RT-qPCR validation

Total RNA extraction

The peripheral blood samples and TRizol reagent (Invitrogen, USA) are evenly mixed in a ratio of 1:3, a total of 12 ml. Total RNA was extracted by using TRizol reagent according to the manufacturer’s instructions.

Library construction and high-throughput sequencing

Total RNAs from peripheral blood of three patients with COPD and their corresponding control patients without COPD were extracted for high-throughput sequencing. Whole transcriptome libraries were constructed by the TruSeq Ribo Profile Library Prep Kit (IL, USA), according to the manufacturer’s instructions. In brief, 10 mg total RNA was depleted rRNA with an Illumina Ribo-Zero Gold kit and purified for end repair and 50-adaptor ligation. Then, reverse transcription was performed with random primers containing 30 adaptor sequences and randomized hexamers. Finally, the cDNAs were purified and amplified with a thermo cycler. The PCR products of 300–500 bp were purified, quantified, and stored at -80℃ before sequencing. The libraries were subjected to 150 nt paired-end sequencing with Illumina novaseq 6000.

Transcriptome sequencing data analysis

To analyze mRNAs and lncRNAs, the high-throughput sequencing tools, Hisat2, and feature_counts, were used to map clean reads to Homo sapiens reference genome (hg19) and calculate the gene expression level which was normalized to FPKM. To determine the differentially expressed mRNA and lncRNA, the “DEseq2” package in R software was used with the corresponding cutoff (q < 0.05, |log2(fold change)|> 1 for mRNA and lncRNA). To predict circular RNA (circRNA), we identified the possible candidates with find_circ, the junction reads were normalized to TPM. Standard of P < 0.05, and|log2(fold change)|> 1 was used to identify differentially expressed circRNAs. The prediction of microRNA response element (MRE) and RNA binding protein (RBP) among differentially expressed circRNAs were analyzed by the CircRNA database (CSCD) (http://gb.whu.edu.cn/CSCD/).

Small RNA data analysis

For small RNA (sRNA) sequencing, six sRNA libraries were generated with TruSeq small RNA (IL, USA) according to the manufacturer’s instructions. Then the prepared libraries were sequenced with an Illumina novaseq 6000. After filtering out the reads shorter than 15 nt, the remaining reads were mapped to the human genome (hg19) and the miRNA database in miRBase with bowtie (− v 1). The differentially expressed miRNAs were determined by DEseq2 with the cutoff of P < 0.1, |log2(fold change)|1.

Construction of lncRNA-mRNA co-expression network

To construct co-expression network between significantly differentially expressed lncRNA and mRNA. Firstly, the Pearson’s correlation coefficients between lncRNA and mRNA were calculated by cor function in R software. A criterion of the coefficient parameter |R|> 0.99 was used for the remaining RNAs to further construct the network. For visualization of the network, the Cytoscape software (https://cytoscape.org, v3.4.0) was used.

Construction of circRNA-miRNA-mRNA network

To construct a circRNA-miRNA-mRNA network between differentially expressed circRNA, miRNA, and mRNA. Firstly, 6 up-regulated circRNAs were verified by RT-qPCR. Next, the miRanda (http://www.microrna.org/microrna/home.do) and TargetScan (http://www.targetscan.org/vert_72/) were used to identify miRNAs targeting on these circRNAs. The miRDB (http://www.mirdb.org/) and miRTarBase (http://miRTarBase.mbc.nctu.edu.tw/) were applied to predict the target mRNAs of miRNAs. Then, the Pearson correlation coefficient (R) in these interactions were calculated by the cor function in R software. These correlation coefficient were used to estimate the co-expression relationship between circRNAs, mRNAs, and miRNAs. The strong interactions between miRNAs, circRNAs, and mRNAs (with a threshold of |R|> 0.99) were selected to construct a circRNA-miRNA-mRNA network using Cytoscape (https://cytoscape.org, v3.4.0).

Gene set enrichment analysis (GSEA)

The GSEA analysis was performed to annotate the function of protein-coding genes (PCGs) in the peripheral blood of COPD patients. These PCGs were pre-ranked and the GSEA analysis was implemented in GSEA software (version 4.1.0, https://www.gsea-msigdb.org).

Reverse transcription and real-time quantitative PCR

cDNA was synthesized using the GoScript Reverse Transcription System (Promega, USA) according to the manufacturer’s protocol. Quantitative real-time PCR was performed with GoTaq SYBR Green qPCR Master Mix (Promega, USA) on a PikoReal 96 real-time PCR system followed by 40 amplification cycles (Thermo Fisher Scientific, USA) according to standard procedures. Actually, all amplification curves already reached the stationary stage before 35 amplification cycles, and the Ct value was obtained at the exponential stage. Relative RNA level was normalized to GAPDH mRNA. All primers are shown in Additional file 3: Table S1.

Statistical analysis

The nortest (v 1.0.4) R package was used to assess whether the data follows normal distribution in all experiments. Student’s t-tests were used to calculate P-values by t.test function in R software. The values reported in the graphs represent average of actual number of independent experiments, with error bars showing SD. The cor function in R software was used to calculate correlation coefficients of all groups. The cluster analysis in heatmaps was conducted by cluster (v 2.1.2) R package. The heatmaps were generated by pheatmap R package (v 1.0.12). The bar plots and scatter plots were generated by ggplot2 (v 3.3.3). After analysis of variance with F-tests with R software, the statistical significance and P-values were evaluated with Student’s t-tests.

Results

Identification of differentially expressed mRNAs and lncRNAs in COPD

To identify COPD-related RNAs, we used fresh peripheral blood samples from COPD patients for RNA sequencing. Fresh peripheral blood from three patients without COPD was used as controls (Fig. 1A). Patient characteristics are shown in Table 1. All the COPD patients were clinically stable without any acute infection and were not on corticosteroid therapy. RNA-seq data of sufficient quality was obtained for mRNA, lncRNA, miRNA, and circRNA expression profiles. We first set out to more specifically investigate the differentially expressed mRNAs and lncRNAs with filtration criteria (fold changes ≥ 2.0 and p-values ≤ 0.05) (Fig. 1B, C). We identified 282 mRNAs and 146 lncRNAs that were differentially expressed in COPD (Additional file 4: Table S2, Additional file 5: Table S3). Among them, 115 mRNAs were upregulated, 167 mRNAs were downregulated; 64 lncRNAs were upregulated and 82 lncRNAs were downregulated in COPD. Accumulative evidence suggests that lncRNAs can interact with co-expressed mRNAs through the formation of complementary hybrids, working from both nearby and distant sites. We, therefore, constructed the lncRNA-mRNA co-expression to reveal the co-regulatory relationship network of COPD based on their potential multi-reciprocal interactions (Fig. 1D). We performed the correlation analysis for mRNAs and lncRNAs by calculating the Pearson correlation coefficient in all samples and selected lncRNA-mRNA pairs with |R|> 0.99 for co-expression network construction. Of these interactions, most coding genes are associated with distinct lncRNAs. A few coding genes were regulated by multiple lncRNAs, such as SCARA5, SNX7, KLHL4, VMO1, SKIV2L, C1orf167, GZMB, NUTM2F, PCDHGA9, AC068775.1, AC104389.6, CLIC6, and COLCA2. These lncRNA-mRNA co-expression networks provided the potential functional interactions of lncRNA-mRNA for future studies.

Fig. 1
figure 1

Analytical procedure and the analysis of mRNA and lncRNA in COPD patients and healthy individuals. The analytical procedure of COPD patient and health individual. A Hierarchical cluster heatmaps display differentially expressed transcripts among mRNA (B) and lncRNA (C). D The co-expression network between differentially expressed lncRNA and mRNA (correlation coefficient absolute value > 0.99). The boxes represent mRNA, the circles represent lncRNA. Red, upregulated; blue, downregulated

GSEA analysis of the differentially expressed mRNAs

To investigate the potential RNA-related pathways and biological processes that these differentially expressed mRNAs involved in COPD, we performed GSEA analysis with either upregulated or downregulated mRNAs (Fig. 2). Of 3972 gene sets, 3756 gene sets (94.56%) were downregulated in the COPD group compared to the normal group. 444 gene sets are significantly enriched at FDR < 25%, and 396 gene sets are significantly enriched at nominal p-value ≤ 1% (Additional file 8: Tables S6, Additional file 9: Table S7). According to KEGG pathway enrichment analysis, downregulated gene sets in COPD were most enriched in six pathways: “ncRNA metabolic process”, “ncRNA processing”, “ribosome biogenesis”, “rRNAs metabolic process”, “tRNA metabolic process” and “tRNA processing”, all of which were downregulated in COPD group. 465 genes were enriched in “ncRNA metabolic process”, among which the top five genes correlated with COPD were TDRD1, LYAR, HENMT1, TRIM71, and SLFN13; 378 genes were enriched in “ncRNA processing”, among which the top five genes correlated with the COPD were LYAR, TSEN54, HSD17B10, ELAC1, and NSA2; 300 genes were enriched in “ribosome biogenesis”, among which the top five genes correlated with the COPD were LYAR, MRPL11, RSL24D1, NSA2, and TSR2; 232 genes were enriched in “rRNAs metabolic process”, among which the top five genes correlated with the COPD were LYAR, SLFN13, NSA2, TSR2, and NHP2; 180 genes were enriched in “tRNA metabolic process”, among which the top five genes correlated with the COPD were SLFN13, TSEN54, HSD17B10, ELAC1 and TARS1; 127 genes were enriched in “tRNA processing”, among which the top five genes correlated with the COPD were TSEN54, HSD17B10, ELAC1, TRMT61B, and OSGEPL1.

Fig. 2
figure 2

GSEA analysis from COPD RNA-seq. A Heatmap of gene expression level in six GO terms related to RNA processes. B Enrichment plot of above six GO ontologies among RNA processes

Validation of differentially expressed mRNA and lncRNAs in clinical samples

To validate the differentially expressed mRNAs and lncRNAs from the COPD RNA-seq, we randomly chose the differentially expressed mRNAs and lncRNAs respectively for RT-qPCR validation with COPD patient sample pairs (Fig. 3). The housekeeping gene GAPDH was used as the endogenous control. Patient characteristics are shown in Table 2. In the mRNA group, CLGN, GPR42, and RSAD2 were significantly upregulated in the blood of COPD patients compared to the normal controls; HLA-DPA1 and TMEM17 were significantly downregulated compared to the normal control (Fig. 3A). In the lncRNA group, AC099332.1 and LINC02573 was significantly upregulated in COPD patients, while AL445493.3, HCG27, and AP000465.1 were significantly downregulated in COPD patients, compared to the normal controls (Fig. 3B). The present results showed that the RT-qPCR verification was in accordance with RNA-seq results with high confidence.

Fig. 3
figure 3

Validation of mRNAs and lncRNAs identified from RNA-seq data by RT-qPCR. The validation results of (A) differentially expressed mRNAs, and (B) differentially expressed lncRNAs. Data are median SD. *p < 0.05, **p < 0.01 are based on the Student’s t-test

Analysis of the differentially expressed miRNAs and circRNAs in COPD

CircRNAs and miRNAs, as well as their functional interactions, play vital roles in the regulation of many biological processes including COPD. We then set out to analyze the differentially expressed circRNAs and miRNAs with filtration criteria (circRNA: fold changes ≥ 2.0 and p values < 0.05; miRNA: fold changes ≥ 2.0 and p values < 0.1) (Fig. 4A, B) [19,20,21]. 85 miRNAs were differentially expressed, with 60 miRNAs (70.59%) were downregulated, while 25 miRNAs (29.41%) were upregulated (Additional file 6: Table S4). Totally 81 circRNAs were differentially expressed, and a heatmap was performed to show either upregulated or downregulated circRNAs. 47 circRNAs (58.02%) were downregulated, while 34 circRNAs (41.98%) were upregulated (Additional file 7: Table S5). To validate the differentially expressed circRNAs from COPD RNA-seq, we chose 6 differentially expressed circRNAs for RT-qPCR validation with COPD patient sample pairs (Fig. 4C). The housekeeping gene GAPDH was used as the endogenous control. Patient characteristics are shown in Table 2. CircFCHO2, circMBOAT2, circPTPN22, circTBC1D22A, cirCOPDADM, and circCKAP5 were significantly upregulated in the blood of COPD patients compared to the normal control, which was all in accordance with RNA-seq results with high confidence.

Fig. 4
figure 4

The analysis of circRNA and miRNA in COPD. A Heatmap display of differentially expressed circRNAs. B Volcano plot display differentially expressed transcripts of miRNAs. C Validation of 6 differentially expressed circRNAs by RT-qPCR

Construction of circRNA-miRNA-mRNA network of COPD

CircRNAs participate in a variety of biological processes through multiple mechanisms due to their unique structures and properties. We then investigated the above 6 validated circRNAs based on their structures, and performed functional prediction (Fig. 5A). Three criteria for circRNAs were mainly focused: MRE (MicroRNA Element) for potential miRNA binding capabilities; RBP (RNA Binding Protein) for circRNA-RBP interaction; ORF (Open Reading Frame) for circRNA as a potentially translatable template [22,23,24,25,26]. CircRNAs with miRNA binding sites could bind directly to the corresponding miRNAs to inhibit miRNA activity and thus regulate the expression of target genes [22, 27]. All the 6 circRNAs have miRNA binding sites, indicating they might function through miRNA sponging, which is also the most reported circRNA function. CircRNAs can act as scaffolds or decoys of RBPs, and thus exert their functions through interactions with the RBPs. All the circRNAs examined but circPTPN22 have predicted RBP binding sites, suggesting they might function as RNA-binding protein scaffold. Some circRNAs could act as translation templates, and their polypeptide products have been shown to play roles in physiology and disease. All the circRNAs have ORF, which indicated that they possessed capability for translation and may serve as functional polypeptides. Of note, circRNAs might function through one or multiple mechanisms depending on specific interacting factors and/or various contexts, which need further investigations and verifications for the predicted mechanisms.

Fig. 5
figure 5

Prediction of the functional mechanism of differentially expressed circRNAs and their functional network in COPD. A Functional mechanism prediction of 6 differentially expressed circRNAs. MRE, miRNA response elements. RBP, RNA binding protein. ORF, open reading frame. B Construction of circRNA-miRNA-mRNA network for 6 differentially expressedcircRNAs according to the interactions among circRNAs, mRNAs, and miRNAs (correlation coefficient absolute value > 0.99)

To investigate the potential function of circRNAs acting as miRNA sponges in COPD, we hence constructed a circRNA-miRNA-mRNA network based on our validated 6 circRNAs and their corresponding regulating miRNAs and mRNAs in COPD. 18 miRNAs and 16 mRNAs were found to be correlated with this ceRNA network (Fig. 5B; Additional file 10: Table S8), which indicates that these identified six circRNAs might function through interaction directly or indirectly with these miRNAs and mRNAs.

Discussion

COPD is a major heterogeneous disease, which becomes one of the leading causes of death worldwide[28]. Bronchodilators and anti-inflammatory agents provide the major treatment for COPD patients, however, with poor treatment outcome. Better understandings of the pathophysiology of disease development are recognized to be of increasing clinical importance, and identification of critical therapeutic targets specifically for COPD is of utmost importance [29]. RNA therapeutics provide high specificity and represent safer and reversible alternatives to DNA-based gene therapies, and RNA drugs potentially offer unique opportunities to expand the range of therapeutic targets with several of them have been approved or being currently in clinical trials [30]. This study uses RNA-seq to preliminarily explore the differentially expressed mRNAs, lncRNAs, miRNAs and circRNAs during the process of COPD, and tentatively reveal the interactions and potential functional mechanisms of these RNAs. We hypothesized that certain dysregulated RNAs could exert their functions either individually or cooperatively, and some of them could be exploited as potential biomarkers.

High-throughput sequencing technologies have facilitated the detection of aberrantly expressed both protein-coding and noncoding genes in humans at the transcriptome level [31]. The number of studies focusing on the disease-related RNAs is fast increasing as differential expression of specific genes would positively or negatively correlate with disease pathology. Increasing evidence long suggests that lncRNAs, miRNAs, and circRNAs along with mRNAs could massively involved in the progression of many diseases, with some of them being identified as potentially suitable biomarkers [32, 33]. However, there is only a limited number of studies for the universal investigations on the differentially expressed RNAs of COPD, especially circRNAs [34].

In this study, we utilized RNA deep-sequencing from peripheral blood of COPD patients and identified 282 mRNAs, 146 lncRNAs, 85 miRNAs, and 80 circRNAs that were differentially expressed. To the best of our knowledge, it is the first comprehensive study that identifies and analyzes the differentially expressed mRNA, miRNA, lncRNA, and circRNA in COPD. We also validated some of the differentially expressedRNAs which were in accordance with the RNA-seq with high confidence. We also performed GSEA analysis with either upregulated or downregulated differentially expressed RNAs. Several critical pathways that associate with COPD were found in the analysis, especially noncoding RNA metabolic processing, RNA processing, ribosome biogenesis, rRNA metabolic process, tRNA metabolic process, and tRNA processing. It is acknowledged that these biological processes might be correlated with the progression of COPD, yet with limited evidence. For example, m6A RNA methylation regulators contribute to COPD progression, and the expressions of m6A RNA methylation regulators (IGF2BP3, FTO, METTL3, and YTHDC2) [35], which have significant associations with some key genes enriched in the signaling pathway and biological processes that promote the development progression of COPD, are highly correlated with the occurrence of COPD [36, 37]. Of note, six mRNAs (LYAR, SLFN13, TSEN54, HSD17B10, ELAC1, NSA2) listed in the six pathways of GESA analysis were significantly downregulated in at least three pathways, in which LYAR, TSEN54, and ELAC1 were enriched in at least four pathways. LYAR is a cell growth-regulating nucleolar protein, which is reported to potentiate rRNA synthesis and is important for influenza A virus replication [38]. LYAR is also found down-regulated proteins in well-differentiated normal human primary bronchial/tracheal epithelial cells compared with undifferentiated cells [39]. TSEN54 is a putative causal gene for lung function and is associated with the genetics of smoking behavior, lung function, and COPD (even in nonsmokers) [40]. Further studies are worthy of investigating whether the differentially expressed RNAs function in these biological processes [41]. In addition, we also found pathways known to be closely related to COPD pathology in the GSEA analysis, such as oxidative phosphorylation, protein folding, cytolysis, amino acid activation, rRNA transcription, ribosomal large subunit biogenesis [42,43,44,45,46,47] (Fig. S2).

lncRNA-mRNA interaction is one of the most frequently studied functional mechanisms of lncRNA [48]. We also constructed a lncRNA-mRNA co-expression network in COPD, based on the ceRNA mechanism which is similar to circRNA-miRNA interaction. It is long acknowledged that the RNA-seq-based lncRNA-mRNA network has become a useful tool to predict functional lncRNAs and their potential functional mechanisms, which is widely used in many disease models and clinical samples [49]. Our lncRNA-mRNA network in COPD would hopefully provide useful information for future COPD investigations.

Lines of evidence reveal that circRNAs participate in various biological processes, including cell proliferation, protein metabolism, autophagy, tumor immunology, signal transduction, genome stability, etc. [29, 50]. In this study, we presented 6 circRNAs that were also validated in the following experiments, and we demonstrated their potential mechanisms according to their nucleotide sequences. CircRNAs are generally found to function with multiple mechanisms, among them, the three main functions are miRNA sponge, protein scaffold, transcriptional regulation, and translation [51]. Of note, circRNA acting as miRNA sponge was first found in 2013 and has become one of the most investigated mechanisms in almost all disease models [50]. It is believed that circRNA, which harbors miRNA binding sites, would sequester miRNA(s) like a sponge and thus alleviate miRNA’s suppression of its target protein(s) [16]. Some circRNAs can act as RBP scaffolds or decoys to exert their function through RNA-RBP interaction. Also, some endogenous circRNAs can be used as scaffolds to regulate protein–protein interactions. Some other circRNAs can initiate peptide translation through internal ribosomal entry sites (IRES) or N6-methyladenosine (m6A) [52, 53]. Based on this, we first constructed the circRNA-miRNA-mRNA network in COPD, and 6 circRNAs, 18 miRNAs, 16 mRNAs were potentially correlated and were involved in the occurrence and progression of COPD, which indicates potential biological or clinical mechanisms for future studies. However, there is a certain possibility that these circRNAs might functionally interact with proteins or translate to polypeptides in other contexts due to the cell- and tissue-specificity of circRNAs [51]. Together with the lncRNA-mRNA co-expression network, we believe that a comprehensive understanding of the complex networks of interactions that these differentially expressed RNAs coordinate would provide a unique opportunity for better therapeutic interventions.

The findings of the study reveal potentially important RNAs, pathways and RNA networks that are associated with COPD, and provide a prospective view towards the underlying functions and functional mechanisms of these differentially expressed RNAs. Compared to the previous work in COPD, we have first conducted comprehensive bioinformatics analysis towards RNAs, especially circRNA, an emerging RNA molecule with unique criteria and clinical potential. To date, few studies have investigated the differential expression of multiple RNAs in COPD context, and this study provides valuable RNA resources for future in-depth COPD investigations. We have also explored the functional mechanisms based on these analyzed RNAs, and most importantly, established the interaction networks, especially the circRNA-miRNA-mRNA network in COPD, which provides sufficient novelty to this study.

However, some further questions remain that require future investigations. First, more COPD samples are expected for our RNA sequencing and bioinformatics analysis based on sample size estimation, which will optimize the study design for differential expression detection with high confidence. Second, it should be noted that we used peripheral blood from COPD patients for RNA-seq rather than bronchus due to the lack of these samples in clinical practices. The direct associations of the RNAs analyzed and COPD should be critically investigated in the future study.Third, the RNA regulatory networks were only based on bioinformatics predictions and lack practical experiments for verification which requires future in-depth studies with in vivo and in vitro experiments. The use of COPD animal models is also a useful tool in COPD investigations, which would facilitate the research on potential functions and evolutionary conservation of the RNAs. Our team is currently committed to studying these potential functions and mechanisms.

Availability of data and materials

The datasets generated and analyzed in the current study are available in the GEO (Gene Expression Omnibus). The accession number is GEO: GSE198740.

Abbreviations

ceRNA:

Competing endogenous RNAs

circRNA:

Circular RNA

COPD:

Chronic obstructive pulmonary disease

FEV1:

Forced expiratory volume in 1 s

FVC:

Forced vital capacity

GSEA:

Gene set enrichment analysis

lncRNA:

Long non-coding RNA

miRNA:

MicroRNA

mRNA:

Message RNA

ncRNA:

Non-coding RNA

qRT-PCR:

Quantitative real-time PCR

RNA-seq:

RNA-sequencing RNA

References:

  1. Agustí A, Vogelmeier C, Faner R. COPD 2020: changes and challenges. Am J Physiol Lung Cell Mol Physiol. 2020;319:L879-l883.

    PubMed  Article  CAS  Google Scholar 

  2. Disease Gifcol: Global strategy for the diagnosis, management, and prevention of chronic obstructive pulmonary disease (2021 REPORT). 2020.

  3. Vogelmeier CF, Román-Rodríguez M, Singh D, Han MK, Rodríguez-Roisin R, Ferguson GT. Goals of COPD treatment: focus on symptoms and exacerbations. Respir Med. 2020;166: 105938.

    PubMed  Article  Google Scholar 

  4. Nambiar S, Bong How S, Gummer J, Trengove R, Moodley Y. Metabolomics in chronic lung diseases. Respirology. 2020;25:139–48.

    PubMed  Article  Google Scholar 

  5. Poller W, Dimmeler S, Heymans S, Zeller T, Haas J, Karakas M, Leistner DM, Jakob P, Nakagawa S, Blankenberg S, et al. Non-coding RNAs in cardiovascular diseases: diagnostic and therapeutic perspectives. Eur Heart J. 2018;39:2704–16.

    CAS  PubMed  Article  Google Scholar 

  6. Sakornsakolpat P, Prokopenko D, Lamontagne M, Reeve NF, Guyatt AL, Jackson VE, Shrine N, Qiao D, Bartz TM, Kim DK, et al. Genetic landscape of chronic obstructive pulmonary disease identifies heterogeneous cell-type and phenotype associations. Nat Genet. 2019;51:494–505.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  7. Morrow JD, Chase RP, Parker MM, Glass K, Seo M, Divo M, Owen CA, Castaldi P, DeMeo DL, Silverman EK, Hersh CP. RNA-sequencing across three matched tissues reveals shared and tissue-specific gene expression and pathway signatures of COPD. Respir Res. 2019;20:65.

    PubMed  PubMed Central  Article  Google Scholar 

  8. Hardin M, Silverman EK. Chronic obstructive pulmonary disease genetics: a review of the past and a look into the future. Chronic Obstr Pulm Dis. 2014;1:33–46.

    PubMed  PubMed Central  Google Scholar 

  9. CorreiadeSousa M, Gjorgjieva M, Dolicka D, Sobolewski C, Foti M. Deciphering miRNAs’ Action through miRNA Editing. Int J Mol Sci. 2019;20:6249.

    Article  CAS  Google Scholar 

  10. Hobbs BD, Tantisira KG. MicroRNAs in COPD: small molecules with big potential. Eur Respir J. 2019;53:1900515.

    CAS  PubMed  Article  Google Scholar 

  11. Xu H, Sun Q, Lu L, Luo F, Zhou L, Liu J, Cao L, Wang Q, Xue J, Yang Q, et al. MicroRNA-218 acts by repressing TNFR1-mediated activation of NF-κB, which is involved in MUC5AC hyper-production and inflammation in smoking-induced bronchiolitis of COPD. Toxicol Lett. 2017;280:171–80.

    CAS  PubMed  Article  Google Scholar 

  12. Faiz A, Steiling K, Roffel MP, Postma DS, Spira A, Lenburg ME, Borggrewe M, Eijgenraam TR, Jonker MR, Koppelman GH, et al. Effect of long-term corticosteroid treatment on microRNA and gene-expression profiles in COPD. Eur Respir J. 2019;53:1801202.

    CAS  PubMed  Article  Google Scholar 

  13. Quinn JJ, Zhang QC, Georgiev P, Ilik IA, Akhtar A, Chang HY. Rapid evolutionary turnover underlies conserved lncRNA–genome interactions. Genes Dev. 2016;30:191–207.

    PubMed  PubMed Central  Article  Google Scholar 

  14. Gu C, Li Y, Liu J, Ying X, Liu Y, Yan J, Chen C, Zhou H, Cao L, Ma Y. LncRNA-mediated SIRT1/FoxO3a and SIRT1/p53 signaling pathways regulate type II alveolar epithelial cell senescence in patients with chronic obstructive pulmonary disease. Mol Med Rep. 2017;15:3129–34.

    CAS  PubMed  Article  Google Scholar 

  15. Bi H, Zhou J, Wu D, Gao W, Li L, Yu L, Liu F, Huang M, Adcock IM, Barnes PJ, Yao X. Microarray analysis of long non-coding RNAs in COPD lung tissue. Inflamm Res. 2015;64:119–26.

    CAS  PubMed  Article  Google Scholar 

  16. Zhou S, Jiang H, Li M, Wu P, Sun L, Liu Y, Zhu K, Zhang B, Sun G, Cao C, Wang R. Circular RNA hsa_circ_0016070 is associated with pulmonary arterial hypertension by promoting PASMC proliferation. Mol Ther Nucleic Acids. 2019;18:275–84.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  17. Guan R, Wang J, Cai Z, Li Z, Wang L, Li Y, Xu J, Li D, Yao H, Liu W, et al. Hydrogen sulfide attenuates cigarette smoke-induced airway remodeling by upregulating SIRT1 signaling pathway. Redox Biol. 2020;28: 101356.

    CAS  PubMed  Article  Google Scholar 

  18. Ma H, Lu L, Xia H, Xiang Q, Sun J, Xue J, Xiao T, Cheng C, Liu Q, Shi A. Circ0061052 regulation of FoxC1/Snail pathway via miR-515-5p is involved in the epithelial–mesenchymal transition of epithelial cells during cigarette smoke-induced airway remodeling. Sci Total Environ. 2020;746: 141181.

    CAS  PubMed  Article  Google Scholar 

  19. Hrdlickova R, Toloue M, Tian B. RNA-Seq methods for transcriptome analysis. Wiley Interdiscip Rev RNA 2017;8.

  20. Lu S, Zhu N, Guo W, Wang X, Li K, Yan J, Jiang C, Han S, Xiang H, Wu X, et al. RNA-Seq revealed a circular RNA-microRNA-mRNA regulatory network in hantaan virus infection. Front Cell Infect Microbiol. 2020;10:97.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  21. Sheng Z, Wang X, Xu G, Shan G, Chen L. Analyses of a panel of transcripts identified from a small sample size and construction of rna networks in hepatocellular carcinoma. Front Genet. 2019;10:431.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  22. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495:384–8.

    CAS  PubMed  Article  Google Scholar 

  23. Kristensen LS, Andersen MS, Stagsted LVW, Ebbesen KK, Hansen TB, Kjems J. The biogenesis, biology and characterization of circular RNAs. Nat Rev Genet. 2019;20:675–91.

    CAS  PubMed  Article  Google Scholar 

  24. Chen L, Huang C, Shan G. Circular RNAs in physiology and non-immunological diseases. Trends Biochem Sci. 2022;47:250–64.

    PubMed  Article  CAS  Google Scholar 

  25. Pamudurti NR, Bartok O, Jens M, Ashwal-Fluss R, Stottmeister C, Ruhe L, Hanan M, Wyler E, Perez-Hernandez D, Ramberger E, et al. Translation of CircRNAs. Mol Cell. 2017;66:9-21.e27.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  26. Yang Y, Fan X, Mao M, Song X, Wu P, Zhang Y, Jin Y, Yang Y, Chen LL, Wang Y, et al. Extensive translation of circular RNAs driven by N(6)-methyladenosine. Cell Res. 2017;27:626–41.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  27. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495:333–8.

    CAS  PubMed  Article  Google Scholar 

  28. Global, regional, and national incidence, prevalence, and years lived with disability for 354 diseases and injuries for 195 countries and territories, 1990–2017: a systematic analysis for the Global Burden of Disease Study 2017. Lancet. 2018;392:1789–1858.

  29. Greene J, Baird AM, Brady L, Lim M, Gray SG, McDermott R, Finn SP. Circular RNAs: biogenesis, function and role in human diseases. Front Mol Biosci. 2017;4:38.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  30. Kim RY, Sunkara KP, Bracke KR, Jarnicki AG, Donovan C, Hsu AC, Ieni A, Beckett EL, Galvão I, Wijnant S, et al. A microRNA-21-mediated SATB1/S100A9/NF-κB axis promotes chronic obstructive pulmonary disease pathogenesis. Sci Transl Med. 2021;13:eaav7223.

    CAS  PubMed  Article  Google Scholar 

  31. Ambardar S, Gupta R, Trakroo D, Lal R, Vakhlu J. High throughput sequencing: an overview of sequencing chemistry. Indian J Microbiol. 2016;56:394–404.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  32. Hombach S, Kretz M. Non-coding RNAs: classification, biology and functioning. Adv Exp Med Biol. 2016;937:3–17.

    CAS  PubMed  Article  Google Scholar 

  33. Wei JW, Huang K, Yang C, Kang CS. Non-coding RNAs as regulators in epigenetics (review). Oncol Rep. 2017;37:3–9.

    PubMed  Article  Google Scholar 

  34. Chen S, Yao Y, Lu S, Chen J, Yang G, Tu L, Chen L. CircRNA0001859, a new diagnostic and prognostic biomarkers for COPD and AECOPD. BMC Pulm Med. 2020;20:311.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  35. Huang X, Lv D, Yang X, Li M, Zhang H. m6A RNA methylation regulators could contribute to the occurrence of chronic obstructive pulmonary disease. J Cell Mol Med. 2020;24:12706–15.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  36. Seki M, Komuro A, Takahashi M, Nashimoto M. Transcription from the proximal promoter of ELAC1, a gene for tRNA repair, is upregulated by interferons. Biochem Biophys Res Commun. 2021;585:162–8.

    CAS  PubMed  Article  Google Scholar 

  37. Fleet JC, Replogle RA, Reyes-Fernandez P, Wang L, Zhang M, Clinkenbeard EL, White KE. Gene-by-diet interactions affect serum 1,25-dihydroxyvitamin D levels in male BXD recombinant inbred mice. Endocrinology. 2016;157:470–81.

    CAS  PubMed  Article  Google Scholar 

  38. Su L, Hershberger RJ, Weissman IL. LYAR, a novel nucleolar protein with zinc finger DNA-binding motifs, is involved in cell growth regulation. Genes Dev. 1993;7:735–48.

    CAS  PubMed  Article  Google Scholar 

  39. Lu XN, Ju GJ, Wang YX, Wang YL, Wang K, Chen JL, Cai W, Zang QW. LYAR promotes the proliferation of non-small cell lung cancer and is associated with poor prognosis. Folia Histochem Cytobiol. 2021;59:282–90.

    CAS  PubMed  Article  Google Scholar 

  40. Störk T, Nessler J, Anderegg L, Hünerfauth E, Schmutz I, Jagannathan V, Kyöstilä K, Lohi H, Baumgärtner W, Tipold A, Leeb T. TSEN54 missense variant in Standard Schnauzers with leukodystrophy. PLoS Genet. 2019;15: e1008411.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  41. Wain LV, Shrine N, Miller S, Jackson VE, Ntalla I, Soler Artigas M, Billington CK, Kheirallah AK, Allen R, Cook JP, et al. Novel insights into the genetics of smoking behaviour, lung function, and chronic obstructive pulmonary disease (UK BiLEVE): a genetic association study in UK Biobank. Lancet Respir Med. 2015;3:769–81.

    PubMed  PubMed Central  Article  Google Scholar 

  42. Lerner CA, Sundar IK, Rahman I. Mitochondrial redox system, dynamics, and dysfunction in lung inflammaging and COPD. Int J Biochem Cell Biol. 2016;81:294–306.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  43. Kenche H, Baty CJ, Vedagiri K, Shapiro SD, Blumental-Perry A. Cigarette smoking affects oxidative protein folding in endoplasmic reticulum by modifying protein disulfide isomerase. FASEB J. 2013;27:965–77.

    CAS  PubMed  Article  Google Scholar 

  44. Freeman CM, Han MK, Martinez FJ, Murray S, Liu LX, Chensue SW, Polak TJ, Sonstein J, Todt JC, Ames TM, et al. Cytotoxic potential of lung CD8(+) T cells increases with chronic obstructive pulmonary disease severity and with in vitro stimulation by IL-18 or IL-15. J Immunol. 2010;184:6504–13.

    CAS  PubMed  Article  Google Scholar 

  45. Kuo WK, Liu YC, Chu CM, Hua CC, Huang CY, Liu MH, Wang CH. Amino acid-based metabolic indexes identify patients with chronic obstructive pulmonary disease and further discriminates patients in advanced BODE stages. Int J Chron Obstruct Pulmon Dis. 2019;14:2257–66.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  46. Farre-Garros R, Lee JY, Natanek SA, Connolly M, Sayer AA, Patel H, Cooper C, Polkey MI, Kemp PR. Quadriceps miR-542-3p and -5p are elevated in COPD and reduce function by inhibiting ribosomal and protein synthesis. J Appl Physiol. 1985;2019(126):1514–24.

    Google Scholar 

  47. Wang Z, Maschera B, Lea S, Kolsum U, Michalovich D, Van Horn S, Traini C, Brown JR, Hessel EM, Singh D. Airway host–microbiome interactions in chronic obstructive pulmonary disease. Respir Res. 2019;20:113.

    PubMed  PubMed Central  Article  Google Scholar 

  48. Zhang Y, Tao Y, Liao Q. Long noncoding RNA: a crosslink in biological regulatory network. Brief Bioinform. 2018;19:930–45.

    CAS  PubMed  Article  Google Scholar 

  49. Zhang J, Le TD, Liu L, Li J. Inferring and analyzing module-specific lncRNA-mRNA causal regulatory networks in human cancer. Brief Bioinform. 2019;20:1403–19.

    CAS  PubMed  Article  Google Scholar 

  50. Lei K, Bai H, Wei Z, Xie C, Wang J, Li J, Chen Q. The mechanism and function of circular RNAs in human diseases. Exp Cell Res. 2018;368:147–58.

    CAS  PubMed  Article  Google Scholar 

  51. Wilusz JE. A 360° view of circular RNAs: from biogenesis to functions. Wiley Interdiscip Rev RNA. 2018;9: e1478.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  52. Kong S, Tao M, Shen X, Ju S. Translatable circRNAs and lncRNAs: driving mechanisms and functions of their translation products. Cancer Lett. 2020;483:59–65.

    CAS  PubMed  Article  Google Scholar 

  53. Wang Z, Lei X. Prediction of RBP binding sites on circRNAs using an LSTM-based deep sequence learning architecture. Brief Bioinform. 2021;22.

Download references

Acknowledgements

Not applicable.

Funding

This article was supported by Natural Science Research Project of Anhui Universities (Grant No. KJ2021ZD0028; KJ2021A0204), National Key R&D Program of China (Grant No. 2019YFA0802600), Hefei Municipal Natural Science Foundation (Grant No. 2021037), Research Fund of Anhui Institute of Translational Medicine (Grant No. 2021zhyx-c67) and Collaborative Chinese and Western Medicine Research Project for Major Difficult Diseases (Grant No. 2021zdynjb06).

Author information

Authors and Affiliations

Authors

Contributions

Contributions to the conception and design of the work: DZ and LC; Acquisition, analysis, or interpretation of data: PL, YW, XZ and RL; Drafting of the manuscript: LC, PL, NZ and YW; Critical revision of the manuscript for important intellectual content: LC, DZ, YW; Statistical analysis: YW, CC, DW; Administrative, technical, or material support: LC, DZ, XZ; Supervision: DZ and LC. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Liang Chen or Dahai Zhao.

Ethics declarations

Ethics approval and consent to participate

This study was approved by the Ethics Committee of the Second Affiliated Hospital of Anhui Medical University [YX2021-148(F1)], China. Written informed consent was obtained from each patient for this study.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig. S1.

Prediction of the functional mechanism of differentially expressed circRNAs. Functional mechanism prediction of 6 differentially expressed circRNAs. MRE, miRNA response elements. RBP, RNA binding protein. ORF, open reading frame.

Additional file 2: Fig. S2

. GSEA analysis related to COPD from RNA-seq. Enrichment plot of six GO terms associated with COPD pathology, including oxidative phosphorylation, protein folding, cytolysis, amino acid activation, rRNA transcription, ribosomal large subunit biogenesis.

Additional file 3: Table S1

. Primers used in this study.

Additional file 4.

Differentially expressed mRNAs in this study.

Additional file 5.

Differentially expressed lncRNAs in this study.

Additional file 6.

Differentially expressed microRNAs in this study.

Additional file 7.

Differentially expressed circRNAs in this study.

Additional file 8.

444 gene sets significantly enriched at FDR<25% in GSEA analysis.

Additional file 9.

396 gene sets significantly enriched at nominal p-value ≤1% in GSEA analysis.

Additional file 10.

18 miRNAs and 16 mRNAs in circRNA-miRNA-mRNA network.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Liu, P., Wang, Y., Zhang, N. et al. Comprehensive identification of RNA transcripts and construction of RNA network in chronic obstructive pulmonary disease. Respir Res 23, 154 (2022). https://doi.org/10.1186/s12931-022-02069-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12931-022-02069-8

Keywords

  • COPD
  • circRNA
  • lncRNA
  • miRNA
  • RNA network