Cigarette smoke extract alters genome‐wide profiles of circular RNAs and mRNAs in primary human small airway epithelial cells

Abstract As a novel kind of non‐coding RNA, circular RNAs (circRNAs) were involved in various biological processes. However, the role of circRNAs in the developmental process of chronic obstructive pulmonary disease (COPD) is still unclear. In the present study, by using a cell model of COPD in primary human small airway epithelial cells (HSAECs) treated with or without cigarette smoke extract (CSE), we uncovered 4,379 previously unknown circRNAs in human cells and 903 smoke‐specific circRNAs, with the help of RNA‐sequencing and bioinformatic analysis. Moreover, 3,872 up‐ and 4,425 down‐regulated mRNAs were also identified under CSE stimulation. Furthermore, a putative circRNA‐microRNA‐mRNA network was constructed for in‐depth mechanism exploration, which indicated that differentially expressed circRNAs could influence expression of some key genes that participate in response to pentose phosphate pathway, ATP‐binding cassette (ABC) transporters, glycosaminoglycan biosynthesis pathway and cancer‐related pathways. Our research indicated that cigarette smoke had an influence on the biogenesis of circRNAs and mRNAs. CircRNAs might be involved in the response to CSE in COPD through the circRNA‐mediated ceRNA networks.

Circular RNAs (circRNAs) are a novel group of endogenous noncoding RNAs, which are originated from linear precursor messenger RNA (pre-mRNA) and generated by back splicing, with covalently joined 3′-and 5′-ends. 4 Unlike liner RNAs, circRNAs are more stable because of their covalently closed loop structures. At first, cir-cRNAs were considered as useless products of abnormal splicing, 5 recent investigations discovered that circRNAs were widespread and could function as gene expression regulators at multiple levels.
For instance, increasing studies revealed that circRNAs could impair miRNA mediated gene silencing by serving as miRNA sponges partially through the competitive endogenous RNA (ceRNA) network. 6,7 Moreover, the intronic circRNAs, majorly accumulating in the nucleus, could interact with the elongating Pol II complex and thereby influence gene transcription. 8 Li et al also discovered that intergenic circRNAs could affect Pol II transcription by interacting with small nuclear ribonucleoprotein to promote transcription. 9 Abundant circRNAs have been discovered in various cells, because of the help of high throughput sequencing technologies and bioinformatic analysis. 8,10 Subsequent reports revealed that cir-cRNAs may be involved in several cellular and developmental process of some diseases, such as prion diseases, neurological disorders, atherosclerotic vascular disease risk, as well as different malignant tumours. 6,7,11,12 It is reported that circRNAs could function as potential disease biomarkers in human saliva and as biomarkers for non-small cell lung cancer and ageing. [13][14][15] However, the possible involvement of circRNAs in the pathogenesis of COPD remains elusive. In this study, we investigated the expression patterns of cir-cRNAs and mRNAs under cigarette smoking stimulation, a major environmental contributor to COPD, in primary human small airway epithelial cells (HSAECs), to explore the possible underlying molecular regulation mechanism of circRNAs in COPD.

| Preparation of cigarette smoke extract
We prepared the cigarette smoke extract (CSE) as previously described with some modifications. 16 Briefly, mainstream smoke derived from five 3R4F Kentucky Research cigarettes was drawn slowly into a 50-mL syringe and bubbled through 10 mL of Dulbecco's modified Eagle's medium (DMEM) which was pre-warmed at 37℃. This preparation was next titrated to pH 7.4 and sterilised with a 0.22 μm filter (Millipore, Bedford, MA, USA). The final solution was considered as 100% CSE. Then the 100% CSE was diluted with serum-free cell culture medium to the required concentrations during use.

| RNA isolation and quantitative real-time PCR
Total RNA was isolated from HSAECs using EZNA Total RNA kit I (Omega Bio-tek, USA), as per manufacturer's instructions. cDNAs were synthesised from total RNA using PrimeiScript RT reagent Kit (TaKaRa) following the manufacturer's instructions. Next, the PCR reaction was performed in triplicates with FastStart Essential DNA Green Master (Roche). The relative PCR amplification was determined using the LightCycler ® 96 PCR system (Roche Molecular Systems, USA) as following: first step is preincubation: 95°C for 10 minutes for 1 cycle; second step is 3-step amplification: 95°C for 10 seconds, Tm for 10 seconds and 72°C for 10 seconds for 40 cycles; third step is melting: 95°C for 10 seconds, 65°C for 60 seconds and 97°C for 1 seconds for 1 cycle. Divergent primers were designed and optimised for circRNAs (Table S1). All data were normalised to GAPDH gene expression.  were identified, at least two distinct reads of a putative circular junction in total (summed for three samples) were required to be identified. Otherwise, at least three distinct reads in total were required.

| DEGs and GO enrichment analysis
To analyse differentially expressed genes (DEGs), we used DESeq R

| CircRNA-associated-ceRNA network prediction
To construct a putative circRNA-mediated competing endogenous RNA (ceRNA) network, we identified all microRNA (miRNA) targeting circRNA and all miRNA targeting gene based on TargetScan and miRanda. Then, if identical miRNA could bind to both circRNA and gene, the targeted circRNAs were defined as candidate ceRNA of gene, and the circRNA-miRNA-mRNA thus represent a candidate ceRNA pairs. 18 Besides, circRNA and gene in candidate ceRNA pairs are required to have the same direction in DEG analysis, because of a positive correlation in expression of ceRNA pairs.

| Expression pattern of circRNA and mRNA in CSE-induced stress
We achieved about 55.2 million (from 50.3 to 60.5 million) and 53.8  Table 2.

| Validation of differentially expressed circRNA and time series analysis
To verify the RNA sequencing data, we selected six identified cir-cRNAs for the analysis of RT-qPCR, including two up-regulated circRNAs and four down-regulated circRNAs. The results showed that two circRNAs (hsa_circ_0010023 and hsa_circ_0040929) were overexpressed, whereas four circRNAs (hsa_circ_0060927, hsa_circ_0002472, hsa_circ_0002153, hsa_circ_0001573) were underexpressed. The data from RT-qPCR were consistent with the RNA-seq data ( Figure 3).
Although, circRNA expression was measured at 24, 48 and 72 hours after CSE treatment, the 48 hours measurements gave the strongest enrichment for the majority of selected circRNAs ( Figure S1). LGALS3BP

| The prediction of circRNA-mediated ceRNA networks
CircRNAs could function as miRNA sponges and inhibit miRNA activity mediated by ceRNA network. In the present study, we pre-

| D ISCUSS I ON
As a newly identified group of widespread endogenous non-coding Our results indicated that the differential expression patterns of circRNAs may be related to their involvement in the process of COPD with smoke inhalation.
Data from mRNA sequencing were also informative. Plenty of studies performed RNA-seq to explore non-coding RNA and mRNA expression profiles in COPD patients' lungs and blood samples. We compared our list of differentially expressed mRNA with those previously reported in COPD lung tissues and blood samples. The mRNAs including IL6, S100A9 and SPP1, which were differentially expressed in our study, were also identified in COPD patients' lung tissues. 23 Furthermore, it is reported that FAM13A in the top up-regulated mRNAs could increase emphysema susceptibility of mice exposing to cigarette smoke by promoting the degradation of β-catenin, 24 a protein involved in cigarette smokeinduced inflammatory cytokine production and airway inflammation in COPD. 16 A previous study also suggested that FAM13A polymorphisms might be associated with COPD susceptibility in tobacco-smoke exposure. 26 In the top GO processes of increased mRNAs, autophagy identified in BP was reported to promote lung epithelial cell death, airway dysfunction in the experimental models of COPD exposed to cigarette smoke in vivo and in vitro. 27 In addition, as a part of endoplasmic reticulum stress, IRE1-mediated unfolded protein response was suggested involved in CSE-induced apoptosis in alveolar type II epithelial cells. 28 In the GO processes of decreased mRNAs, regulation of signal transduction by p53 class mediator process and extracellular matrix organisation process were identified in BP whereas histone serine kinase activity and RAGE receptor binding were found in MF, which is corresponding to the previous studies. 29 Concerning KEGG pathway,

ATP-binding cassette (ABC) transporters in increased mRNAs 30
showed their importance in CSE-induced stress. All these findings provided novel clues for COPD study.
Studies using different models have shown that miRNAs expression could be affected by direct exposure to cigarette smoke or exposure to CSE in the lungs. For instance, a decreased miR200c was found to regulate an increase in epithelial to mesenchymal transition (EMT) during NF-κB-mediated inflammation in human bronchial epithelial cells induced by CSE, and the increased EMT is associated with pathological processes and tissue remodelling. 31 Another miRNA molecular, miR135b was reported to serve as a counter-regulatory mechanism by regulating the IL-1 receptor expression for cigarette smoke-induced inflammation in vivo. 32 The possible relationship between miRNAs and the related target genes with COPD has been supported by large numbers of studies using both animal models and human samples, which were recently well summarised. 33 CircRNAs were reported to serve as miRNA sponges and keep miRNA away from its target genes through ceRNA network, 11 suggesting that they might also play important roles in the development of COPD associated with smoke inhalation, considering the close linkage between miRNA and COPD. The circRNA-mediated ceRNA networks were predicted in several types of human cells. [34][35][36] In the present study, we also constructed a putative circRNA-mediated ceRNA network in HSAECs under cigarette smoke stimulation. As shown in Figure 6, circRNAs could influence expression of many genes which were required for pentose phosphate pathway, ATP-binding cassette (ABC) transporters, glycosaminoglycan biosynthesis pathway and other pathways, by interacting with miRNAs. Our results suggested that circRNA-mediated post-transcriptional regulations might be involved in the process of COPD through such miRNAs in airway epithelial cells.
Although we could not conclude whether these circRNAs are specific responses to CSE exposure, comparison of the circRNAs expression profiles among different stimulus such as CSE and LPS is well worth being performed in the future study.
Some limitations in this study should be acknowledged. First, it is reported that different circRNA detection tools result in different sets of circRNAs, and several algorithms combined could contribute to ideally reliable predictions. 37 In the present study, only CIRCexplore2 was used to identify candidate circRNAs, and combination of multiple tools may increase the robustness of cir-cRNA detection. Second, this study only identified the circRNAs and mRNAs expression profiles in small airway epithelial cells.
Although these findings provide an important research basis for the future investigation of the pathological and molecular mechanisms underlying the regulation of the development of COPD by these specific circRNAs and mRNAs molecules, further verification in animals and patients is required.

| CON CLUS ION
In conclusion, a unique set of circRNAs and mRNAs expression patterns were found in primary HSAECs treated with CSE and these

ACK N OWLED G EM ENTS
We thank Deqing Yang for assistance with qRT-PCR. Thanks to Yuanzhi Cheng from Chengdu Basebiotech Co., Ltd for providing assistances on bioinformatic analysis. Sequencing was performed at Novogene Technology Co., Ltd.

CO M PE TI N G I NTER E S TS
All authors declare that they have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data used to support the findings of the current study are available from the corresponding authors upon request.