Ten male clinical participants who smoked, including five stable COPD patients and five healthy controls, were selected on November 2019. It should be noted that the 10 smokers were in a state of abstinence when they participated in the experiment. The inclusion criteria were as follows: (1) COPD diagnosed using the ratio of postbronchodilator the forced expiratory volume in the first one second to the forced vital capacity of the lungs (FEV1/FVC). Combined with clinical manifestations, if postbronchodilator FEV1/FVC <70%, the patients were diagnosed with COPD. (2) The age of male smokers was between 50 and 80 years old, and the number of years of smoking was between 30 and 50 years. (3) The participants were from Yunnan Province to exclude the impact of different altitudes from different provinces on COPD. (4) There were no other respiratory diseases, including asthma, bronchiectasis, pulmonary fibrosis, lung tumor and tuberculosis. Healthy control individuals had no clinical symptoms of any respiratory disease, such as cough, sputum, and wheezing.

All the data were expressed as the mean ± standard deviation (SD). This study was approved by the Ethics Committee of The Affiliated Hospital of Yunnan University (No. 2020139), and each participant signed the informed consent form before the collection of blood samples.

Peripheral blood samples (20 mL) were collected from the median cubital vein from all participants. Then, they were placed into 5 mL purple-head EDTA-K2 anticoagulation tubes and mixed gently by inversion for subsequent experiments. Based on the manufacturer’s instructions for Human Lymphocyte Separation Medium (DAKEWEI, Beijing, China), PBMCs were isolated and extracted within 4 h.

Total RNA was extracted from PBMCs using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s instructions. After successful RNA quality and integrity testing, the ribosomal RNA removal kit was used to remove ribosomal RNA (rRNA) from the total RNA, and then the RNA was broken into fragments of approximately 300 bp in length by ion interruption. cDNA libraries were generated by polymerase chain reaction (PCR) amplification for the enrich library fragments. Then, the Agilent 2,100 Bioanalyzer was used for quality inspection of the library, including the total concentration and effective concentration. After the quality inspection results were passed, next-generation sequencing (NGS) technology was used to perform paired-end sequencing of these libraries based on the Illumina HiSeq sequencing platform (Personalbio, Shanghai, China).

After data filtering (Cutadapt) and quality evaluation (FastQC), the RNA-Seq data were compared to the reference genome (TopHat2), and the data source for the reference genome was the Ensembl database (http://www.ensembl.org/, Homo_sapiens.GRCh38.dna.primary assembly.fa). Then, we used HTSeq to compare the read count value of each gene as the original expression of the gene and Fragments Per Kilobase of transcript per Million mapped reads (FPKM) to normalize the expression. These RNA-Seq data include circRNAs and mRNAs.

After filtering the original offline data, the sequences from these data were compared to that of the reference genome of humans, and the expression of each gene was calculated. R software (version number: 3.4.1, http://www.r-project.org/) was used to normalize and standardize the data, and the DESeq R package was used for differentially expressed gene (DEG) analysis (Sun et al., 2013). Under the conditions of a |log2FoldChange| value > 1 and a significant p-value < 0.05, differentially expressed (DE) circRNAs (DEcircRNAs) and mRNAs (DEmRNAs) were identified. To show the characteristics of DEGs, volcano maps and hierarchical bidirectional clustering analysis were generated using the ggplots2 R package. In addition, we used the circlize R package to label DEGs in the genome in a genomic circle diagram based on the genomic information and the RNA differential expression analysis results.

Gene Ontology (GO, http://www.geneontology.org/) functional enrichment analysis of DEGs covers three domains, i.e., biological processes (BPs), cellular components (CCs), and molecular functions (MFs). The Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.kegg.jp/) contains a gene information database, a system information database, and a chemical information database. The topGO and pathview R packages were used for GO enrichment analysis, and a p-value <0.05 was defined as statistically significant.

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org/) was used to construct the protein–protein interaction (PPI) networks of integrated DEGs (Szklarczyk et al., 2019). The minimum required interaction score was 0.15. The obtained results were imported into Cytoscape version 3.6.1. The cytoHubba plugin provides a user-friendly interface to explore important nodes in PPI networks by the multiscale curvature classification (MCC) algorithm (score ≥5), and it was used to select the hub node, which is named the core gene or core protein (Chin et al., 2014).

The Connectivity Map (CMap), a database of biological applications, revealed the relationship of small-molecule drugs, gene expression levels and interrelated diseases (Gao et al., 2019). In other words, the gene expression profile was used to establish the association between genes, diseases and drugs, which should help scholars to quickly identify highly associated genes with a disease, identify the main chemical structure of a molecule, and summarize the possible directions of the mechanism of drug molecules.

The regulatory relationships of the DEcircRNA-DEmiRNA pairs and DEmiRNA-DEmRNA pairs were predicted based on the CircInteractome (https://circinteractome.nia.nih.gov/), StarBase (http://starbase.sysu.edu.cn/), TargetScan (http://www.targetscan.org/), miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/), miRDB (http://www.mirdb.org/), miRanda (http://www.microrna.org/) and miRBase (http://www.mirbase.org/) databases (Fan et al., 2016; Fan et al., 2018). Then, the intersection of predicted target genes and DEmRNAs was analysed by the custom Venn diagram tool. The circRNA acts as a competitive endogenous RNA (ceRNA) to competitively bind with miRNA and to regulate the expression of downstream genes and proteins. According to the ceRNA theory, a circRNA-miRNA-mRNA ceRNA regulatory network was constructed by Cytoscape version 3.6.1 (Shannon et al., 2003), integrating circRNA-miRNA pairs and miRNA-mRNA pairs.

Total RNA was extract from PBMCs by TRIzol (Invitrogen, USA). cDNA was synthesized using a reverse transcription kit (Vazyme Biotech Co., Ltd., Nanjing, China). qPCR was conducted using the SYBR Master Mix Kit (Vazyme, Nanjing, China) in an ABI 7500 PCR System (Applied Biosystems, CA, USA). The PCR conditions were 95°C for 5 min, followed by 40 cycles of 95°C for 15 s and 60°C for 30 s. Data were normalized using small nuclear RNA 6 (U6) for circRNAs, GAPDH for mRNAs as a reference gene. The relative expression was calculated using the 2^−ΔΔCT method. The primer for RT-qPCR in Table 1.

The statistical analyses were performed using R software and SPSS. Among them, R software was used for the visualisation of the results of GO enrichment, KEGG function and differential gene analysis, and t-tests for RT-qPCR were performed with SPSS version 24.0 (SPSS, IL, USA). Additionally, independent samples t-tests were carried out using SPSS for age and years of smoking, and p > 0.05 was deemed to be not statistically significant.

The characteristics of the ten eligible participants are presented in Table 2, including five COPD patients and five healthy controls. All participants were male smokers with good general physical fitness from Yunnan Province. In the COPD group, the age range was 53–76 years, the pack-years of smoking ranged from 40 to 52 years, and the FEV1/FVC values were all <70%. In the control group, the age range was 58–76 years, the pack-years of smoking ranged from 30 to 48 years, and the FEV1/FVC values were all >70%. There were no statistically significant differences in age (66.00 ± 8.60 years vs. 65.40 ± 6.73 years, p > 0.05) or number of years of smoking (44.00 ± 4.95 years vs. 38.00 ± 8.03 years, p > 0.05) between the two groups.

Correlation analysis suggested that the correlation coefficient was highly strong (between 0.8-1), indicating that the similarity of the expression patterns between samples was high (Figures 1A, D). Principal component analysis (PCA) determined patterns between similar samples, and if the distance is closer, the similarity between the samples is higher. Our results showed good similarity between samples in Figures 1B, E. In addition, the distributions of the datasets for the circRNAs/mRNAs of the ten samples were nearly identical after normalization, as shown in violin plots (Figures 1C, F).

By comparing the COPD group with the healthy control group using the screening criteria, there were 114 significant DEcircRNAs, in which 57 were upregulated (including the top 10 genes: hsa_circ_032857, hsa_circ_019208, hsa_circ_024906, hsa_circ_035125, hsa_circ_010451, hsa_circ_019753, hsa_circ_031856, hsa_circ_040189, hsa_circ_021622, and hsa_circ_003470), and 57 were downregulated (including the top 10 genes: hsa_circ_007446, hsa_circ_016302, hsa_circ_019528, hsa_circ_019015, hsa_circ_038146, hsa_circ_026711, hsa_circ_026970, hsa_circ_001331, hsa_circ_012036, and hsa_circ_012997), and 58 significant DEmRNAs, in which 21 were upregulated (including the top 10 genes: KCNJ12, SHE, GCOM1, NRCAM, MYZAP, TMEM178B, SEMA5A, VSTM4, MYOM2, and MYL4), and 37 were downregulated (including the top 10 genes: GDF15, RNF17, HYDIN, HLA-DQA1, RSPO4, MUC19, TCF7L1, COL26A1, HOXA10, and C7orf61) that were identified (Supplementary Table S1). DEcircRNAs and DEmRNAs were confirmed by a heatmap (Figures 2A, D) and volcano plot, respectively (Figures 2B, E). In addition, we marked DEGs on the genomic circle diagram to display the position of the DEGs on the chromosome in detail (Figures 2C, F).

GO functional enrichment analysis showed that the main pathways enriched for biological process primarily were the response to resveratrol, response to metformin, positive regulation of interleukin (IL)-18 production, positive regulation of IL-5 secretion, multicellular organismal signalling, positive regulation of NLRP3 inflammasome complex assembly, positive regulation of endothelial cell chemotaxis, response to nitrite, cellular response to nitrite, cell-substrate adhesion, and axonal fasciculation. The cellular component included the cytoplasmic side of the plasma membrane, cell surface, axonal surface, symmetric synapse, side of the membrane, and myosin complex, major histocompatibility complex (MHC) class II protein complex, and MHC protein complex. The MFs included IL-2 binding, talin binding, protein-coupled neurotensin receptor activity, triplex DNA binding, transketolase, transferase activity, and transferring aldehyde or ketones groups Figure 3A shows the top 10 clusters of GO functional enrichment analysis of DEmRNAs.

KEGG pathway enrichment analysis found that the signalling pathways related to COPD included the differentiation of Th1, Th2 and Th17 cells, apelin signalling pathway, AMPK signalling pathway, Wnt signalling pathway, JAK-STAT signalling pathway, and PI3K-Akt signalling pathway. Figure 3B shows the top 20 clusters of KEGG pathway enrichment analysis of DEmRNAs.

As shown in Figure 4A, the DEmRNAs are plotted as a PPI network showing the interaction between multiple proteins. In the PPI map, the core genes in the network are MYO16, MYL4, SCN4A, NRCAM, HMCN1, MYOM2, and IQSEC3.

Based on DEGs between the COPD patients and healthy controls, the CMap database was utilized to screen candidate small-molecule drugs. Suitable parameters included a significant p-value < 0.05 and an absolute enrichment value > 0.80. As Table 3 shows, there were negative correlations between hycanthone, ascorbic acid, canavanine, puromycin and lincomycin; therefore, they were considered to be potential drugs for the treatment of COPD.

The ceRNA regulatory network of circRNAs-miRNAs-mRNAs (Figure 4B) contains 60 DEcircRNAs (23 upregulated, 37 downregulated), 39 predicted DE miRNAs (DEmiRNAs; 28 upregulated, 11 downregulated) and 13 DEmRNAs (5 upregulated, 8 downregulated). There were 225 ceRNA pairs in the network. As shown in Supplementary Table S2, the negative regulatory relationships among the DEcircRNAs, DEmiRNAs, and DEmRNAs in the ceRNA network are listed in detail.

To verify the RNA sequencing data, we selected six identified circRNAs for the analysis of RT-qPCR including three upregulated circRNAs and three downregulated circRNAs, and six candidate mRNAs including three upregulated mRNAs and three downregulated mRNAs. The results showed that hsa_circ_006155, hsa_circ_031080, hsa_circ_041434, MYZAP, NRCAM and SHE were increased, whereas hsa_circ_007446, hsa_circ_019015, hsa_circ_038146, HYDIN, MUC19, and RNF17 were decreased in COPD patients versus normal control samples. The results were consistent with the RNA-Seq data (Figure 5).