Identification and Integrated Analysis of circRNA and miRNA of Radiation-Induced Lung Injury in a Mouse Model

Introduction

Radiation-induced lung injury (RILI) is characterized by acute pneumonitis and chronic fibrosis, both of which could be lethal, occurring in approximately 10–20% of patients who undergo thoracic radiotherapy.¹ Modern radiation therapy has witnessed evolution in radiation planning and delivery, successfully reducing the incidence of lung injury; however, it is still a troubling condition in clinical settings.² RILI is involved in a complex process, including various of molecular and cellular interactions that ultimately lead to large fibroblast proliferation, accumulation, and differentiation, thus causing excessive deposition of the extracellular matrix and pulmonary fibrosis.³

Different from linear RNAs, the structure of circular RNAs (circRNAs) are covalently closed loops. As a result, circRNAs do not contain either polyadenylated tail or 5′–3′ polarity. CircRNAs could act as a regulator of multiple biological processes, including cancer development, metabolism, and immunity regulation.⁴ MicroRNAs (miRNAs) are small, highly conserved non-coding RNA molecules involved in the regulation of gene expression. Recent studies have reported that circRNAs could act as miRNA sponges to compete for the same binding site of miRNAs, thus inhibiting miRNAs combining to their target mRNAs and suppressing the target genes.⁵ Various studies suggested that the circRNA–miRNA–mRNA axis is probably involved in the signaling pathway of many human diseases by regulating pathogenicity-related gene expression. Analysis of circRNA-related ceRNA networks has been performed in cancer,⁵ cardiology analysis,⁶ and even radiation-injury mouse model.⁷ Luo et al successfully investigated the circRNA–miRNA–mRNA pathways associated in a radiation-induced esophageal injury (RIEI) model and revealed that sphingolipid metabolism may be involved in the RIEI process.⁷

However, no report focused on the circRNA–miRNA–mRNA network and pathways of intercellular communication in the RILI process. Understanding the molecular mechanism by revealing circRNA–miRNA–mRNA competitive regulatory networks is of great importance. As the recent improvement of RNA-sequencing (RNA-seq) facilitated obtaining huge and accurate results,⁸ this technique has become a proper method to detect irradiation-responsive circRNAs in RILI. In a previous study, the lncRNA–miRNA–mRNA network was found to contribute to the RILI process by RNA-seq.⁹ In contrast to the findings of this previous study, the role of circRNA–miRNA–mRNA was identified in RILI with different mechanisms from the collagen fibril organization function. Therefore, to spot the potential role of circRNAs in RILI, RNA-seq was first performed in the present study to detect irradiation-responsive genes in a RILI mouse model, and the genome-wide expression between irradiated and normal lung tissues were compared. The function of differentially expressed genes and circRNAs were also analyzed through Gene Ontology (GO) to identify the pathways during intercellular communication in the RILI process. Taken together, the results could provide important resources and information for further RILI research.

Materials and Methods

Animals

An animal model of RILI was established in accordance with previous studies, and it was verified by H&E staining and Masson staining.^10,11 Female Balb/c mice aged 8 weeks old were purchased from Shanghai SLAC Laboratory Animal Co. Ltd (Shanghai, China), In total, 10 mice were housed, with five per cage, and kept under standard laboratory conditions (22 °C ± 2 °C, 55% ± 10% humidity, and 12/12 hours of light/dark cycle). All animal experimental procedures and protocols were conducted in accordance with the guidelines of the Institutional Animal Care and Use Committee and approved under reference number 20150404A114 by the Department of Laboratory Animal Science.

Irradiation Protocol

The 10 mice were evenly distributed into two groups: control (non‐irradiated) group and IR (irradiated) group. The mice in the control group only received sham irradiation, while those in the IR group received thoracic irradiation. The radiation of mice was provided by Small Animal Radiation Research Platform (FUSCC, Shanghai, China). All mice were induced and maintained by anesthesia with isopentane (increasing from 0.5% to 2.0% in 8min for inducing and 1.5% for maintaining, depending on the status of the animal). Then, they were positioned with a self-designed box to fix their neck and stretch the body. For the IR group, the chest of each mouse was irradiated with a single dose of 12 Gy exposure by using a 6 MV linear accelerator (Siemens Primus-Hi). For the control group, each mouse received 0 Gy exposure. The dose rate was 2 Gy/min. The source–surface distance was 1 m, and the size of the radiation field was 2.5×15 cm. The head and abdomen were shielded by 2.5 cm-thick lead blocks to protect other body parts from irradiation. Three mice from each group were randomly selected, and all these six mice were killed by cervical dislocation at 48 h after irradiation. Once a mouse was killed, the lung was taken and frozen in −80 °C fridge immediately. The same protocol was followed for PCR validation using another six mice.

H&E Staining and Evaluation

A part of the lung samples from three irradiated and three control mice was processed to make paraffin blocks, and then the tissue sections were deparaffinized by xylene and rehydration in a series of alcohol wash. Afterwards, the tissue sections were processed via the H&E staining reported previously as follows:¹² 10 min of hematin, 1–3 s of washing, 1–3 s of 1% ethanol hydrochloride, 10–30 s of washing, 1–3 min of 0.5% eosin, 1–2 s of washing, 1–2 s of 80% alcohol, 2–3s of 95% alcohol, 3–5 s of 95% alcohol, 5–10 min of 100% alcohol, 5–10 min of xylene carbonate, 2 min of xylene (I), 2 min of xylene (II), 2 min of xylene (III). After drying for some time, neutral balsam was used for mounting. The H&E results were evaluated in accordance with the previously reported scoring system.¹³ Pneumonia alterations were defined as thickened alveolar walls, widening of the interlobular septa, alveolar septal edema and massive inflammatory cell infiltration in the alveoli. Three randomly selected fields of each mouse were scored by a pathologist in Fudan University Shanghai Cancer Center and scored as follows: 0, no pneumonia; 1+, pneumonia alterations in 1–20% of the field; 2+, pneumonia alterations in 21–50% of the field; and 3+, pneumonia alterations in 51–100% of the field.

RNA Library Construction and Sequencing

The lung samples were treated with TRIzol (Invitrogen, Carlsbad, CA, USA) to isolate and purify the total RNA within a month after irradiation. NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA) was used to quantify the RNA amount and purity. RNA integrity was assessed using Agilent 2100 with RIN number >7.0. Ribo-Zero rRNA Removal Kit (Illumina, San Diego, USA) was used to deplete ribosomal RNA. The remaining RNAs were fragmented into small pieces by divalent cations at high temperature. Afterwards, cDNAs were created by reverse transcription of the cleaved RNA fragments. The cDNAs were then used for synthesizing U-labeled second-stranded DNAs with RNase H, Escherichia coli DNA polymerase I, and dUTP. Subsequently, an A-base was added to the blunt ends of each strand to prepare them for ligating to the indexed adapters, with each containing a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single-or dual-index adapters were ligated to the fragments, and AMPureXP beads were used to perform size selection. After a heat-labile UDG enzyme was used on the U-labeled second-stranded DNAs, PCR was conducted under the following settings: 3 mins, 95 °C for initial denaturation; 15 secs, 98°C for 8 cycles of denaturation; 15 secs, 60°C for annealing; 30 secs, 72°C for extension; and 5 mins, 72°C for final extension. The insert size of the final cDNA library was approximately 300 bp in average. Finally, Illumina Hiseq 4000 (LC Bio, China) was used to perform paired-end sequencing.

The quantity and purity of each RNA sample was verified using Bioanalyzer 2100, with RIN number >7.0, and RNA 6000 Nano LabChip Kit (Aglient, Santa Clara, CA, USA), with 28S/18S ratio between 1.8 and 2.2. Approximately 1 μg of the total RNA was used to prepare miRNA library in accordance with the protocol of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA).

Differential Expression Analysis of circRNAs and miRNAs

Cutadapt¹⁴ was used to remove the reads with low quality or undetermined bases and adaptor contamination. After the sequence quality was verified by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), Bowtie2¹⁵ and topaht2¹⁶ were used to map reads to the mouse genome. The remaining reads (unmapped reads) were still mapped to the genome by using tophat-fusion. CIRCExplorer was used to denovo assemble the mapped reads to circRNAs at first. Then, back splicing reads were identified in the unmapped reads by tophat-fusion and CIRCExplorer. All samples generated unique circRNAs. The circRNA expression levels from different samples or groups were calculated by scripts in house. Differential expression was defined as an srpbm log2 (fold change) >1 or log2 (fold change) < −1, with p<0.05, by R package – EdgeR.

As for miRNAs, raw reads were subjected to ACGT101-miR (LC Sciences, Houston, Texas, USA) to remove low complexity, junk, adapter dimers, common RNA families (snRNA, snoRNA, rRNA, and tRNA), and repeats. Subsequently, knownand novel miRNAs were identified in miRbase 22.0 by using BLAST search. Finally, differentially expressed miRNAs were defined in accordance with normalized deep-sequencing counts by using Chi-squared nXn test, Fisher exact test, Student t test, Chi-squared 2×2 test, or ANOVA, with p < 0.05.

CircRNA–miRNA–mRNA Network Construction and Functional Analysis

The target miRNAs of circRNAs were predicted to explore the function of circRNAs as ceRNAs by using MiRanda (3.3a) and Targetscan (5.0) in accordance with the software protocol, with a maximum binding free energy of less than −20 and a Targetscan score of more than 50. Afterwards, the miRNAs that could bind with differentially expressed circRNAs were selected to further predict the target mRNAs. DAVID 6.8 Bioinformatics Resources was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations. Then, the significantly differentially expressed gene symbols were mapped into the terms in GO database and pathways in KEGG database. On this basis, four parameters were calculated, including S (significantly differentially expressed gene symbols of certain GO term or pathway), TS (total gene symbols that are significantly differentially expressed), B (total gene symbols of certain GO term or pathway), and TB (total gene symbols). Then, hypergeometric test was performed to obtain the most significantly enriched terms or pathways compared with the genome background. Then, Cytoscape3.5.1 was used to display the circRNA–miRNA–mRNA networks.

RT-PCR

Following the manufacturer’s instructions, total RNA was extracted with TRIzol reagent (Invitrogen, USA). For RNA RT-PCR, a First-strand cDNA Synthesis Kit (Thermo Scientific, USA) was used for reverse transcription in accordance with the manufacturer’s instructions, and then the reactions were processed using a Realtime PCR System (Applied Biosystems 7500, Carlsbad, CA, USA) with SYBR Premix Ex Taq Kit (TaKaRa, Japan). Expression was normalized to internal controls GAPDH (RiboBio, Guangzhou, China). The primer sequences are listed in Table S1.

Statistics

FPKM¹⁷ was calculated by StringTie to quantify the expression of mRNAs and circRNAs. Differential expression was defined as FPKM log2 (fold change) >1 or log2 (fold change) <-1, with p<0.05, by R package – Ballgown.¹⁸ Differentially expressed miRNAs were defined in accordance with normalized deep-sequencing counts by using Chi-squared nXn test, Fisher exact test, Student t test, Chi-squared 2×2 test, or ANOVA, with p < 0.05. The expression levels of circRNAs and mRNAs examined by RT-PCR were calculated via delta-delta-ct method and normalized as ratios between the IR group and the control group. Z-score was used for heat map: Z sample-i = [(log2(signal sample-i)-mean (log2(signal) of all samples)][standard deviation (log2(signal) of all samples)]. Student t test or Mann–Whitney U-test was used to compare the expression levels of circRNAs and mRNAs examined by RT-PCR depending on whether the data obeyed the normal distribution or not.

Results

Histologic Evaluation of Acute Radiation-Induced Lung Injury

By using H&E staining, the lung tissue from radiated mice at 48 h after radiation was compared with that of the control group to identify the model of RILI. As described in the method, three randomly selected fields of each mouse were evaluated. For the nine fields of three mice in the IR group, seven scored 3+, while the other two fields had a score of 2+. As for the control group, all the nine fieldsscored 0. As shown in the representative images in Figure 1, compared with the control group, the IR group showed thickened alveolar walls, alveolar septal edema, and massive inflammatory cell infiltration in the alveoli. These pathological changes is consistent of previous report of acute RILI induced by a single dose of 12 Gy chest radiation.¹⁹

Figure 1 H&E staining of lung tissue from: (A) control group and (B) radiation-induced lung injury (IR) group at 48 h after radiation. Compared with the control group, the irradiated lungs showed thickened alveolar walls, alveolar septal edema, and massive inflammatory cell infiltration in the alveoli.

Expression Profile of circRNAs

In accordance with the methods previously described, 7103 circRNAs were detected in the six mice in total. Differently expressed circRNAs with a cutoff fold change ≥2 or ≤0.5, along with p < 0.05, were analyzed to determine the changes in circRNAs after RILI. As a result, 27 significantly dysregulated circRNAs were found. Among these circRNAs, 17 were up-regulated significantly, while 10 were down-regulated significantly. The results were shown as volcano plots and heatmaps to visually illustrate the expression profile differences between RILI and normal lung (Figure 2).

Figure 2 CircRNA profiling in radiation-induced lung injury (IR) group and control group. The (A) heatmap (the colored scale is Z-score) and (B) volcano plot displayed 17 up-regulated and 10 down-regulated circRNAs. Red represents up-regulation and blue represents down-regulation. The top five fold-change gene names were marked. CircRNAs only detected in one group are not shown in the volcano plot.

Differentially Expressed miRNAs

After the sequencing data of miRNAs were acquired, 785 miRNAs were detected in total in the two groups. Seventy-three miRNAs were only detected in the control group, while 60 miRNAs were only detected in the RILI group (Supplementary Figure 1). Then, the differently expressed miRNAs were analyzed as described previously, with p<0.05. Among these dysregulated miRNAs, 21 were significantly up-regulated, while 33 were significantly down-regulated. Volcano plots and heatmaps were also displayed to visually understand the expression profile differences in RILI from normal lung (Figure 3). The top five up- and down-regulated circRNAs and miRNAs are listed in Tables 1 and 2.

Table 1 Top Five Significantly Up- and Down-Regulated circRNAs

Table 2 Top Five Significantly Up- and Down-Regulated miRNAs

Figure 3 MiRNA profiling in radiation-induced lung injury (IR) group and control group. The (A) heatmap and (B) volcano plot displayed 21 up-regulated miRNAs and 33 down-regulated miRNAs. Red represents up-regulation and blue represents down-regulation. The miRNAs only detected in one group are not shown in the volcano plot.

GO and KEGG Pathway Analysis in circRNAs Functioning as ceRNAs

A series of analysis was performed to predict the respective sponge miRNAs and target mRNA of all the dysregulated circRNAs to further investigate the potential functions of the dysregulated circRNAs via ceRNA networks. A GO analysis of the target mRNAs of these significantly differentially expressed circRNAs was performed. As shown in Figure 4A, positive regulation of alpha-beta T cell proliferation, interstitial matrix, collagen fibril organization, chemokine receptor activity, cellular defense response, and B cell receptor signaling pathway were the main functions associated with dysregulated circRNAs.

Figure 4 (A) Gene Ontology (GO) annotations and (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for mRNAs regulated by circRNA–miRNA–mRNA network. The top 20 according to p value of each analysis are displayed.

Next, KEGG pathway analysis was performed for these target mRNAs of dysregulated circRNAs (Figure 4B). The main enriched pathways include: Th17 cell differentiation, Th1 and Th2 cell differentiation, primary immunodeficiency, leishmaniasis, and hematopoietic cell lineage.

CeRNA Networks Constructed Based on Screened mRNAs and Bioinformatics Prediction

Differentially expressed mRNAs were filtered from the mRNAs of Th1 and Th2 differentiation pathways in KEGG analysis to construct a circRNA–miRNA–mRNA network to further understand the relationship between these circRNAs and RILI. After filtering, a total of 105 cirRNA–miRNA–mRNA pathways were constructed, including 16 cirRNAs, 26 miRNAs, and 22 mRNAs. Considering the complexity in this ceRNA network, a visually cirRNA–miRNA–mRNA network was constructed in accordance with the filtered RNAs (Figure 5). This network may help enhance the understanding on the interaction among RNAs and further investigate the mechanisms in the process of RILI.

Figure 5 CircRNA–miRNA–mRNA network based on the filtered mRNAs associated with collagen fibril organization function in Gene Ontology (GO) analysis. The yellow square represents mRNAs, the green circle represents miRNAs, and the red diamond represents circRNAs. Their regulatory relationships are displayed as lines between them. The red “*” marked the verified genes by qRT-PCR.

Verification of Differentially Expressed circRNAs and mRNAs in Constructed ceRNA Network

Six circRNAs that were up-regulated significantly (p < 0.05) and substantially (fold changes ≥ 2) in RILI were selected; they take part in the ceRNA network related with Th1 and Th2 differentiation pathways, including circRNA4146, circRNA4584, circRNA5229, circRNA544, circRNA1092, circRNA3340. Then, they were verified in radiated and non-radiated mice by RT-PCR. The results showed that circRNA5229, circRNA544, and circRNA3340 were significantly up-regulated in radiated mice (Figure 6A, p = 0.024, p = 0.036, and p = 0.028, respectively), consistent with the sequencing results. Besides, although not significantly, the other circRNA verified by RT-PCR showed the same trend with sequencing results. In addition, to verify the effect of RILI on the ceRNA network associated with Th1 and Th2 differentiation pathways, the ceRNA-predicted mRNAs and proteins associated with the network related to Th1 and Th2 differentiation pathways were examined. The mRNA of Il2ra, Rbpj, Mapk14, CD4, Il12a, Stat1, and Nfatc3 was found to be significantly up-regulated (Figure 6B, p = 0.004, p = 0.025, p = 0.025, p = 0.006, p = 0.01, p = 0.004, and p = 0.004, respectively). Similarly, the IHC results identified the up-regulation of CD4 and ZAP70 which play major roles in the Th1 and Th2 differentiation pathways (Figure 7). Taken together, these results verified the sequencing data and furthermore identified a major regulation of Th1 and Th2 differentiation pathways after RILI.

Figure 6 Verification of the expression level of (A) circRNA and (B) mRNA by RT-PCR. The bars show mean with SD. Significant results are noted by * for 0.01 ≤ p < 0.05 and ** for p < 0.01.

Figure 7 Representative images of mouse lung to examine CD4 and ZAP70, two proteins that play roles in Th1 and Th2 differentiation pathways, in irradiated or control group (original magnification: ×100). For CD4 and ZAP70, higher expression level of proteins were identified in the irradiated group than in the control group. The 50 μm scale bars are marked in the lower left corner.

Discussion

In this work, irradiation-induced changes at the transcription level of miRNAs and circRNAs were first investigated in irradiated and non-irradiated lung tissues by RNA-seq. By comparison of the transcriptome profiles in response to irradiation, 21 significantly up-regulated and 33 down-regulated miRNAs were identified. Among the 27 differentially expressed circRNAs, 10 were down-regulated and 17 were up-regulated. GO analysis of the target mRNAs of these significantly differently expressed circRNAs was then performed. The differentially expressed miRNAs took part in a series of cellular processes, such as positive regulation of alpha-beta T cell proliferation, interstitial matrix, collagen fibril organization, chemokine receptor activity, cellular defense response, and B cell receptor signaling pathway.

Previous studies have revealed that after inducing lung injury, T cell proliferation could be inhibited. In-vitro studies have found that lymphocyte proliferation was limited, including T cells, in alveolar liquid extracted from patients with acute hypersensitivity pneumonitis.²⁰ In addition, some studies have shown part of the underlying mechanisms. An in-vivo study has demonstrated that after physical injury, the intracellular arginine in T cells decreased, thus causing inhibition of proliferation.²¹ As for collagen fibril organization, many studies have focused on the relationship between lung injury and collagen fibril formation or remodeling. For example, Ley K et al have provided evidence that lung fibrosis is mediated by lysophosphatidic acid after injury.²² Meanwhile, injury of alveolar epithelial type II cells (AE2 cells) is tightly correlated with fibrotic remodeling.²³ Interstitial matrix is associated with lung injury as well. Some mechanisms after lung injury have also been reported to be associated with changes in the interstitial matrix, including collagen organization mentioned above. In acute lung injury, restoration and thickening of the interstitial matrix have been reviewed.²⁴ The pathological characteristics of the interstitial compartment following irradiation were also listed and reviewed.²⁵ On the basis of GO results of the differentially expressed miRNAs in the present study, these three functions stood out and interestingly, they are related to lung injury, consistent with the previous studies above. With the help of RNA profile results, the underlying mechanisms could be further investigated.

CircRNAs were reported to have functions such as sequestration of microRNAs or proteins, modulation of transcription and interference with splicing, and even translation to produce polypeptides.²⁶ The functions and clinical effects of ceRNAs were reported in various diseases and thus may present opportunities for new therapeutic approaches.²⁷ In recent years, microRNA and circRNA sequencing techniques were broadly used to recognize potential targets from the complex RNA networks based on the models of various disease.^28–30 CircRNAs sequencing was even performed in an RIEI mouse model and sepsis-induced acute lung injury.³¹ These previous studies have successfully provided meaningful potential target or mechanism for the diseases. In the present study, verification of the results showed that several circRNAs and mRNAs associated with Th1 and Th2 cell differentiation were up-regulated significantly, which is a preliminary confirmation of the immunomodulatory role of related pathways in RILI such as circRNA3340/miR-146a-5p/CD4, circRNA544/miR-188-3p/Il2ra, and circRNA439/miR-702-3p/Il12a. Blocking of these RNAs could thus be a potential targets to modulate immune microenvironment in a beneficial manner after RILI.

Although limitations exist, such as that the timepoint selected and the animals used were few to clearly obtain a landscape of lungs, which consist of various structures and cell types, to the knowledge of the authors, this study was the first to analyze the relationship between RILI and circRNAs by using RNA-seq data and computational approaches to identify the differentially expressed miRNAs and circRNAs. Previous studies have used similar methods to investigate the mechanism of RIEI. The results also provided potential clues to the mechanism of RILI. As the roles of circRNAs in RILI have not yet been fully studied, this analysis could provide valuable data and information for future studies. Unlike many previous RILI studies on pulmonary fibrosis, the present study found that T cell proliferation and differentiation pathway may be particularly important in the study of radiation pneumonia. Its association with circRNAs may be the focus of future research on circRNAs and radiation pneumonia. As reported previously,³² several circRNAs are highly conserved in human and mice. Therefore, future research must focus on finding the conservative circRNAs and verify their function in human. Meanwhile, with the increasingly development of immunotherapy in thoracic tumors,³³ immune-associated pneumonia and radiation pneumonia may play an important role in the study of RILI. The possible interaction of RILI also makes the study of immunity more important in RILI than before.² The results of the present study exactly demonstrated the dysregulation of ceRNAs after a series of crosstalk in the process of RILI, and these ceRNAs have a potential in modulating immunology. However, more in-depth studies are needed due to the limited sample size in the current research. Besides, as various lineage cells exist in lung tissues, further exploration of the relationship and molecular mechanism between cells, such as fibrocytes and immune cells, are needed in future studies.

Conclusions

In conclusion, unique irradiation-induced miRNA and circRNA profiles were primarily identified in RILI. The data showed that they have potential immunological functions during the pathogenesis of RILI. Moreover, the data could provide valuable resources and information for further studies, especially the roles of circRNAs in RILI. The results also revealed the molecular function of associated circRNAs, and further in-depth research, such as immunomics analysis and single-cell sequencing, should provide more valuable resource and information for RILI research.