Identification of lncRNA, MicroRNA, and mRNA-Associated CeRNA Network of Radiation-Induced Lung Injury in a Mice Model.

Radiation-induced lung injury (RILI) can be challenging for thoracic radiotherapy, thus investigating its mechanisms of related pathophysiological process is needed. Long noncoding RNAs (lncRNAs) was found to participate in normal tissue damage induced by ionizing irradiation. Here, we first profiled the dysregulation of lncRNAs, microRNAs (miRNAs), and messenger RNAs (mRNAs) of RILI in mice model receiving 12 Gy thoracic irradiation. The lung tissue was collected 48 hours after irradiation, after which an RNA library was built by RNA sequencing. Compared with the control group, 461 mRNAs and 401 lncRNAs were significantly upregulated, while 936 mRNAs and 501 lncRNAs were significantly downregulated. Then we predicted target miRNAs of the dysregulated lncRNAs and the target mRNAs of these miRNAs. Next, functional annotations of these target mRNAs were performed. Results showed some pathways apparently dysregulated, such as Th1 and Th2 cell differentiation, Th17 cell differentiation, and hematopoietic cell lineage. Through this study, we also highlighted that T helpers could be vital in RILI through lncRNA-miRNA-mRNA network, therefore causing fibrosis, indicating that RNA dysregulation in early stage of RILI may cause severe late complications. Thus, research on the target mechanism and early intervention of lncRNAs with associated competing endogenous RNA network will benefit the treatment of RILI.

Introduction

Radiation therapy constitutes a great part in the treatment of cancer. For those types of cancer which need to receive thoracic radiation, such as non-small cell lung cancer, major limitation exists in delivering higher radiation doses because of the concern for radiation-induced lung injury (RILI). Such dose-limited injury comes from cytokines released after cellular injury, recruitment of inflammatory infiltration, and repairing, which could present in the form of pneumonitis and fibrosis.[1] Many factors, including patient demographic characteristics, tumor characteristics, and dosimetric parameters, may all play vital roles in the risk of RILI.[2]

Noncoding RNAs are identified as RNAs undergoing transcription but not encoding proteins. Among them, long noncoding RNAs (lncRNAs) are RNAs whose length ranges from 200 to 100 000 nt, participating in diverse cellular processes and involved in disease progression.[3] Previous studies have also proved lncRNAs playing a role in DNA damage induced by ionizing irradiation.[4,5] However, changes in lncRNA expression irradiated lungs are less studied, and there is no reports demonstrating the lncRNA-microRNA (miRNA)-messenger RNA (mRNA) interactions in RILI.

In this work, to investigate genes responsive to irradiation in RILI mice model, we used an RNA sequencing (RNA-seq) technique and compared the genome-wide expression between normal lung tissue and the irradiated group. Gene Ontology (GO) and pathway analysis was also performed to identify the function of differentially expressed genes, differentially expressed lncRNAs, and intercellular communication pathways in the RILI process.

Method

Mouse Model of the RILI

Female BALB/c mice, aged 8 weeks, were purchased from Shanghai SLAC Laboratory Animal Co Ltd (Shanghai, China). Mice were housed 5 per cage under standard laboratory conditions (22 ± 2°C, 55 ± 10% humidity, 12-hour light/12-hour dark cycle). Six mice were treated in 2 groups: control group (nonirradiated) and IR group (irradiated). The 3 mice in the control group didn’t receive any irradiation, while the 3 in the irradiated group received 12 Gy thoracic irradiation at room temperature. Radiation of mice was provided by Small Animal Radiation Research Platform (FUSCC, Shanghai, China). Lead shielding was used to protect other body parts from irradiation.

All mice were killed at day 3 after irradiation. Once a mouse was killed, the lung was taken and frozen in −80°C fridge immediately. All experimental procedures and protocols were conducted according to the guidelines of our institutional animal care and use committee.

RNA Library Construction and Sequencing

Lung samples were proceeded by TRIzol (Invitrogen, Carlsbad, California) to isolate and purify total RNA. NanoDrop ND-1000 (NanoDrop, Wilmington, Delaware) was used to quantify RNA amount and purity. RNA integrity was assessed using Agilent 2100 with RIN number >7.0. To deplete ribosomal RNA, Ribo-Zero rRNA Removal Kit (Illumina, San Diego, California) was used. The remaining RNAs were fragmented into small pieces by divalent cations at high temperature. Then, complementary DNAs (cDNAs) were created by reverse transcription of these cleaved RNA fragments. The cDNAs were then used for synthesizing U-labeled second-stranded DNAs with RNase H, Escherichia coli DNA polymerase I, and 2′-deoxyuridine 5′-triphosphate (dUTP). Then, an A-base is added to the blunt ends of each strand, preparing them for ligating to the indexed adapters, which each contains a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single- or dual-index adapters are ligated to the fragments, and AMPureXP beads were used to perform size selection. After the using heat-labile Uracil-DNA Glycosylase (UDG) on the U-labeled second-stranded DNAs, polymerase chain reaction was proceeded with following settings: 3 minutes, 95°C for initial denaturation; 15 seconds, 98°C for 8 cycles of denaturation; 15 seconds, 60°C for annealing; 30 seconds, 72°C for extension; 5 minutes, 72°C for final extension. The insert size of the final cDNA library was proximately 300 bp in average. Finally, Illumina Hiseq 4000 (LC Bio, Hangzhou, China) was used to perform paired-end sequencing.

Different Expression Analysis of MRNA, lncRNAs, and MiRNAs

Cutadapt[6] was used to remove the reads with low-quality or undetermined bases and adaptor contamination. After verifying sequence quality by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), Bowtie2[7] and topaht2[8] were used to map reads to the genome of mice. StringTie[9] was used to assemble these mapped reads of each sample. Then all transcriptomes were merged to reconstruct a comprehensive transcriptome using Perl scripts. The expression levels of these samples were estimated by StringTie and Ballgown.[10]

Transcripts shorter than 200 bp or overlapped with known mRNAs were discarded. Then we predicted transcripts that have coding potential using Coding Potential Calculator (CPC)[11] and Coding-Non-Coding Index (CNCI).[12] After removing transcripts which have CNCI score <0 and CPC score <−1, the remaining transcripts were identified as lncRNAs.

Fragments Per Kilobase of transcript per Million fragments mapped (FPKM)[13] were calculated by StringTie to quantify the expression level of mRNAs and lncRNAs. The differentially expression was defined as an FPKM log2 (fold change) >1 or log2 (fold change) <−1 along with P value <.05 by R package—Ballgown.[10]

As for miRNAs, raw reads were subjected to ACGT101-miR (LC Sciences, Houston, Texas) to remove low complexity, junk, adapter dimers, common RNA families (small nuclear RNA, small nucleolar RNA, rRNA, transfer RNA), and repeats. Subsequently, known miRNAs and novel miRNAs were identified in miRbase 22.0 using BLAST search. Finally, we defined differentially expressed miRNAs according to normalized deep-sequencing counts by chi-square n × n test, Fisher exact test, Student t test, chi-square 2 × 2 test, or analysis of variance with a P value less than .05.

The lncRNA-miRNA-mRNA Network Construction and Functional Analysis

To explore the function of lncRNAs as competing endogenous RNA (ceRNAs), we predicted the target miRNAs of lncRNAs using MiRanda and Targetscan according to the protocol of software with a maximum binding free energy less than −20 along with a Targetscan score more than 50. After the prediction, miRNAs, which could bind with differentially expressed lncRNAs, were selected to further predict the target mRNAs. Then we presented functional analysis of the target genes for lncRNAs by BLAST2GO.[14] And Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to perform pathway analysis. Significance was expressed as a P < .05. Then Cytoscape3.5.1 was used to display the lncRNA-miRNA-mRNA networks.

Result

Expression Profiles of mRNAs and lncRNAs

In total, 39 897 mRNAs and 20 239 lncRNAs were detected in the 6 mice. As shown in Figure 1, the lncRNAs from 6 mice were widely distributed in almost all the chromosomes.

Figure 1.

Long noncoding RNAs (lncRNA) mapping by Circos (www.circos.ca) shows lncRNA distribution in each chromosome. The outer 3 rings are the radiation-induced lung injury group, and the inner 3 are the control group.

To understand the expression profile differences in RILI, we compared the lncRNA and mRNA profiles of the irradiated mice to those in the control group. Compared with the control group, 461 mRNAs were upregulated significantly and 936 were downregulated significantly. For lncRNAs, 401 were upregulated significantly and 501 were downregulated significantly (P < .05, fold change ≥2). The most upregulated lncRNA was Gm37347 with 43-fold change and the most downregulated lncRNA was AC153536.3 with 0.09-fold change. Volcano plots and heat maps of lncRNAs and mRNAs in the 2 groups were displayed to visually show the apparent variations (Figure 2). Besides, the top 5 up- and downregulated lncRNAs and mRNAs are illustrated in Tables 1 and 2.

Figure 2.

The expression profiling in the radiation-induced lung injury group and control group. A total of 461 upregulated messenger RNAs (mRNAs) and 936 downregulated mRNAs are displayed in (A) the heat map and (C) the volcano plot; 401 upregulated long noncoding RNAs (lncRNAs) and 501 downregulated lncRNAs are displayed in (B) the heat map and (D) the volcano plot. Red represents upregulation and blue represents downregulation.

Table 1.

Top 5 Significantly Up- and Downregulated Long Noncoding RNA.

Transcript ID	Fold Change	Regulation	P Value
MSTRG.15609.2	438.73	Up	<.01
MSTRG.31845.13	288.59	Up	<.01
MSTRG.14639.19	260.29	Up	<.01
MSTRG.7054.6	237.58	Up	<.01
MSTRG.5690.1	228.31	Up	<.01
MSTRG.14639.13	0.01	Down	<.01
ENSMUST00000174808	0.01	Down	<.01
MSTRG.28131.2	0.01	Down	<.01
ENSMUST00000174808	0.01	Down	<.01
ENSMUST00000181447	0.01	Down	<.01

Table 2.

Top 5 Significantly Up- and Downregulated Messenger RNAs.

Gene Name	Fold Change	Regulation	P Value
Gm5483	4266.42	Up	<.01
Stfa2l1	444.83	Up	<.01
Ifitm6	188.36	Up	<.01
Asprv1	133.76	Up	<.01
Mmp8	117.05	Up	<.01
Acta1	<.01	Down	<.01
Ighv11-2	<.01	Down	<.01
Ighv12-3	.01	Down	<.01
Pla2g2d	.01	Down	<.01
Cd207	.01	Down	<.01

Function of lncRNAs as ceRNAs via lncRNA-miRNA-mRNA Networks

To further investigate the potential functions of the dysregulated lncRNAs in the RILI via ceRNA networks, we selected all the dysregulated lncRNAs to predict their respective sponge miRNAs and target mRNA. A Venn diagram was used to show the relationship between dysregulated mRNAs and target mRNAs of dysregulated lncRNAs (Figure 3). As shown, 424 target mRNAs of dysregulated lncRNAs were dysregulated. Then a GO analysis of all the target mRNAs of significantly differently expressed lncRNAs was performed. As was shown in Figure 4A, interstitial matrix, collagen fibril organization, cellular defense response, chemokine receptor activity, and B-cell receptor signaling pathway were the main functions associated with dysregulated lncRNAs. A GO analysis of differentially expressed mRNA was also performed and the results are shown in Supplementary Figure 1A.

Figure 3.

Differentially expressed messenger RNAs (mRNAs) and long noncoding RNAs (lncRNAs) and their relationships. The red circle is the 1397 dysregulated mRNAs, the blue one is the 1166 target mRNAs of dysregulated lncRNAs.

Figure 4.

(A) Gene Ontology (GO) annotations and (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for the messenger RNAs (mRNAs) regulated by long noncoding RNA-microRNA-mRNA network. The top 20 according to P value of each analysis are displayed.

We next performed KEGG pathway analysis for these target mRNAs of dysregulated lncRNAs. The main enriched pathways and related genes include Th1 and Th2 cell differentiation (Cd247, Cd3e, Cd3 g, Cd4, H2-Aa, H2-Ab1, H2-DMa, H2-DMb1, H2-DMb2, H2-Eb1, H2-Oa, H2-Ob, Il1b, Il21r, Il23r, Il27ra, Il2ra, Il2rb, Il2rg, Gm20489, Gm614, Irf4, Lat, Lck, Mir8119, Mapk13, Nfatc2, Nfkbie, Prkcq, Tbx21, and Zap70), Th17 cell differentiation (Cd247, Cd3e, Cd3 g, Cd4, H2-Aa, H2-Ab1, H2-DMa, H2-DMb1, H2-DMb2, H2-Eb1, H2-Oa, H2-Ob, Il1b, Il21r, Il23r, Il27ra, Il2ra, Il2rb, Il2rg, Gm20489, Gm614, Irf4, Lat, Lck, Mir8119, Mapk13, Nfatc2, Nfkbie, Prkcq, Tbx21, and Zap70), hematopoietic cell lineage (Cd19, Cd1d2, Cd22, Cd2, Cd33, Cd3e, Cd3 g, Cd4, Cd5, Cd7, Cd8a, Cd8b1, Cr2, Csf2ra, Csf3r, Fcer2a, Flt3, Gp5, Gp9, H2-Aa, H2-Ab1, H2-DMa, H2-DMb1, H2-DMb2, H2-Eb1, H2-Oa, H2-Ob, Il1a, Il1b, Il1r2, Il2ra, Il5ra, Il7r, Il7, Il9r, Itga2b, Itga4, Itgam, Gm45700, Ms4a1, and Siglech), primary immunodeficiency, and leishmaniasis (Figure 4B). A KEGG pathway analysis for differentially expressed mRNA was performed as well (Supplementary Figure 1B).

Competing Endogenous RNA Networks Were Constructed Based on the Screened mRNAs and Bioinformatics Prediction

To better understand the relationship of these lncRNAs and RILI, differentially expressed mRNAs were filtered from mRNAs of collagen fibril organization function in GO analysis to construct an lncRNA-miRNA-mRNA network. After filtering mRNAs, a total of 2735 lncRNA-miRNA-mRNA pathways were constructed including 350 lncRNAs, 15 miRNAs, and 12 mRNAs. Considering the complexity in this ceRNA network, we constructed a visual lncRNA-miRNA-mRNA network according to the filtered RNAs (Figure 5). It may help us better understand the interaction of RNAs and further investigate the mechanisms in the process of RILI.

Figure 5.

Long noncoding RNA (lncRNA)-microRNA (miRNA)-messenger RNA (mRNA) network based on the filtered mRNA associated with collagen fibril organization function in Gene ontology analysis. The yellow squares represent mRNAs, the green circles represent miRNAs, and the small red diamonds represent lncRNAs. Their regulatory relationships are displayed as lines between them.

Discussion

Radiation-induced lung injury is a complex process, including varieties of molecular and cellular interactions, thus leading to large fibroblast proliferation differentiation and accumulation, and ultimately resulting in excessive extracellular matrix deposits and pulmonary fibrosis.[15] As a result, researches are needed to figure out the appropriate biomarkers and molecular mechanisms of RILI.

The improvement in RNA-seq method and microarray method both gave us chances to achieve more accurate results. The RNA-seq method provide the chance to detect more genes when it well correlates with the microarray method. Also, oligo-dT primed RNA also can help to avoid sequencing bias in the 3′ end of the transcript compared to oligo-dT primed cDNA fragmented.[5]

By RNA-seq, we investigated radiation-induced whole transcriptome alterations in nonirradiated and irradiated lung tissues. By comparison of the transcriptome profiles in response to irradiation, we identified 1397 differentially expressed genes between control and IR groups, including 461 upregulated and 936 downregulated genes.

The functional annotation showed that the significantly dysregulated pathways include Th1 and Th2 cell differentiation, Th17 cell differentiation, and hematopoietic cell lineage. Similar results were shown in the KEGG pathway analysis of all the differentially expressed mRNAs, indicating that these pathways were actually influenced by the RILI not only with dysregulated lncRNAs but also with the dysregulation of downstream mRNAs. Notably, according to our results, T helper subtypes seem to play an important role in the process of RILI, which is consistent with previous studies.[16] Meanwhile, other types of lung injury could also cause an increase of T helpers.[17,18] T helpers are CD4+ T cells that contain various subtypes with distinct cytokine production profiles which lead to different functions. Evidences[19,20] have shown that T helpers could contribute to the process of lung fibrosis. Different subtypes of T helpers may play distinct roles in RILI. Interferon-γ−/− mice, which are deficient in Th1 cells, could develop a significantly enhanced fibrosis response, comparing to C57BL/6J mice. In contrast, irradiated Il17−/− mice, which lacked Th17 cells, spared both pneumonitis and fibrosis. These studies showed that the severe fibrosis response to irradiation could be explained by the activation of the Il6-Tgfβ-Il17 and T17 pathway, with infiltration of neutrophil to the lung.[21] It is confirmed by the functional annotation in our study that Th1 and Th17 differentiation is apparently dysregulated through lncRNA-miRNA-mRNA network, which means targeting related lncRNAs or miRNAs could be a potential way to balance the differentiation of T helper subtypes and therefore leads to a better outcome of RILI.

In conclusion, by RNA-seq, we detected the differentially expressed lncRNA and mRNA and constructed the ceRNA network. The functional annotation highlighted that T helpers could play important roles in RILI through lncRNA-miRNA-mRNA network and therefore cause fibrosis in the future, which indicated that RNA dysregulation in early stage of RILI may cause severe late complications. Thus, research on the target mechanism and early intervention of lncRNAs with associated ceRNA network will benefit the treatment of RILI.

Click here for additional data file.

Supplementary_Figure_1 for Identification of lncRNA, MicroRNA, and mRNA-Associated CeRNA Network of Radiation-Induced Lung Injury in a Mice Model by Yida Li, Liqing Zou, Xi Yang, Li Chu, Jianjiao Ni, Xiao Chu, Tiantian Guo and Zhengfei Zhu in Dose-Response