An RNA-seq-Based Expression Profiling of Radiation-Induced Esophageal Injury in a Rat Model

Abstract

Radiation-induced acute injury is the main reason for the suspension of radiotherapy and unsuccessful treatment of cancer. It is of great importance to understand the molecular mechanism of radiation-induced esophageal injury. We used RNA-seq data from normal esophageal tissue and irradiated esophageal tissues and applied computational approaches to identify and characterize differentially expressed genes and detected 40 059 messenger RNAs (mRNAs) previously annotated and 717 novel long noncoding RNAs (lncRNAs). There were 14 upregulated and 32 downregulated lncRNAs among the differentially expressed lncRNA group. Their target genes were involved in the mRNA surveillance pathway, pathological immune responses, and cellular homeostasis. Additionally, we found 853 differentially expressed mRNAs, and there were 384 upregulated and 469 downregulated mRNAs. Notably, we found that the differentially expressed mRNAs were enriched for steroid biosynthesis, the tumor necrosis factor signaling pathway, focal adhesion, pathways in cancer, extracellular matrix–receptor interaction, and so on. The response of normal esophageal tissues to ionizing radiation is multifarious. The radiation-induced cell damage response by multiple pathways followed by pathological immune responses activated. Studies on the dynamic network of molecules involved in radiation-induced esophageal injury are under way to clarify the regulatory mechanisms and identify the candidate targets.

Introduction

Lung and esophageal cancers are the most common cancer. ^1,2 Lung and esophageal cancers are rapidly metastasizing malignancies that are associated with poor patient prognosis and a 5-year survival rate of less than 20%. ^3,4 Treatment includes surgery, adjuvant chemotherapy, and radiation therapy; however, chemoradiotherapy is burdened by many adverse effects. ⁵ Radiation-induced acute injury is the main reason for the suspension of radiotherapy and the unsuccessful treatment of cancer. ^6,7 Esophageal epithelial cells are extremely sensitive to ionizing radiation (IR). Different doses of IR may cause different degrees of esophageal epithelium damage. ⁸ However, IR is also identified as an important risk factor for esophageal cancer cells (ESCC), as noted by the increased incidence of ESCC as a secondary cancer after radiotherapy for primary carcinomas of the head and neck and mediastinal regions. ^{9 -11} The molecular mechanism of radiation-induced injury is still unclear. A comprehensive understanding of the response of normal esophageal tissue to IR and detailed analysis may help identify biomarkers for the diagnosis of esophagitis or targeted drugs for the prevention of esophagitis. Thus, it is of great importance to understand the molecular mechanism of radiation-induced esophagitis.

Ionizing radiation–induced DNA damage and reactive oxygen species (ROS) can kill cells. The IR-induced DNA damage response, such as the activation of checkpoint pathways, can repair damage. ¹² The ROS play central roles in the determination of cell fate as second messengers and modify various signaling molecules (nuclear factor-κB, mitogen-activated protein kinases, Keap1-Nrf2-ARE, and PI3K-Akt), ion channels and transporters (Ca²⁺ and mPTP), and components of the protein kinase and ubiquitination/proteasome system. ^{13 -15} Ionizing radiation-induced esophagitis exhibits a multistep and multifactor process. Epperly et al reported that radiation of the esophagus induces significantly increased levels of messenger RNA (mRNA) for proinflammatory cytokines, such as transforming growth factor β1 (TGF-β1), interleukin-1 (IL-1), tumor necrosis factor α (TNF-α), IL-18, and interferon γ via ROS, and superoxide dismutase-plasmid/liposome decreases the acute and chronic side effects of radiation-induced damage to the esophagus. ¹⁶ Furthermore, the level of epidermal growth factor (EGF) was decreased in the irradiated esophagus in both in vivo and in vitro models, and EGF may be a potential therapeutic intervention strategy to treat radiation-induced esophagitis. ¹⁷ Moreover, a report by Patel et al found that IR enhances esophageal epithelial cell migration and invasion through a paracrine mechanism involving stromal-derived hepatocyte growth factor using cytokine arrays. ¹⁸ Although multiple mechanisms have been proposed for the development of esophagitis, analysis of the molecular signaling events involved is still not comprehensive. To date, advances in high-throughput sequencing methodology have provided a large amount of information regarding gene expression at the transcriptome level as well as the underlying molecular events in response to irradiation. We focused this work on the identification of pathways of intercellular communication and the lncRNA network in the esophagus tissue that are impacted by IR.

In this work, we used an RNA-seq technique to investigate irradiation-responsive genes in normal esophageal tissue. We compared the genome-wide expression between the normal esophageal tissue and the IR group. The functional categories of differentially expressed genes and differentially expressed lncRNAs (DELs) were also analyzed based on gene ontology (GO). Finally, we focused on the TNF signaling pathway and PI3K-Akt signaling pathway, which are activated by irradiation. Together, these data should provide prime information for research on radiotherapy of ESCC. To elucidate the role of radiation exposure in the development of esophagitis, we focused this work on the identification of pathways involved in the intercellular communication in the esophagus that are impacted by IR.

Materials and Methods

Animals and Treatments

Ten male Sprague Dawley rats (4 weeks old) were purchased from the Shanghai SLAC Laboratory Animal Co, Ltd (Shanghai, China). The animals were housed using a 12-hour light–dark cycle and had free access to food and water. The rats were anesthetized with an intraperitoneal injection of ketamine (75 mg/kg) and xylazine (10 mg/kg), and the hair on the rat’s neck–chest region was shaved using a razor. The rats were immobilized with adhesive tape on a plastic plate to minimize motion during radiation exposure. A 3-cm-thick piece of lead was used to shield the rats and localize the radiation field (3 cm × 4 cm). A single dose of 20-Gy irradiation was administered to the esophageal area at a dose rate of 2 Gy/min using 6-MeV x-ray irradiation (Clinac 2100EX; Varian Medical Systems, Inc, California). The control group of rats was sham irradiated. Seven days after irradiation, the esophageal tissues were resected for analysis. All animal experiments were conducted according to legal regulations in China and were carried out with permission and under the regulation of the Institutional Animal Care and Use Committee of the National Institute of Radiological Sciences.

Sample Preparation for RNA-seq

Total RNA was extracted from each sample of esophageal tissue (control group or irradiation group) using TRIzol reagent according to the manufacturer’s protocol (Invitrogen, California). The RNA was quantified using a NanoDrop ND-2000C spectrophotometer (Thermo, California). For RNA high-throughput sequencing, RNA libraries were created from each group using the NEBNext Ultra Directional RNA Library Preparation kit from Illumina (San Diego, CA). The main steps in the workflow involved the removal of ribosomal RNA, the fragmentation of total RNA, reverse transcription and second-strand complementary DNA (cDNA) synthesis, end repair, dA tailing, and adaptor ligation. The products of these reactions were purified and enriched by polymerase chain reaction to create the final cDNA library. The libraries were then sequenced using an Illumina HiSeq2500 (Illumina, paired-end sequencing).

RNA-seq Data Acquisition and Quality Control

We obtained the RNA sequence row data for the control group and irradiated esophageal tissue group by RNA-seq on the Illumina platform. Then, we clipped and trimmed the reads to avoid low-quality data using Trim Galore (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Trimmed sequence files underwent quality control analysis using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We utilized several parameters to evaluate the read quality, including the number of the reads, guanine-cytosine (GC) contents, and average length of the reads.

Differential Expression Analysis of Genes

For the analysis of differentially expressed genes, the clean data for each sample were aligned to the rat reference genome (ftp://ftp.ensembl.org/pub/release-83/fasta/rattus_norvegicus/dna/) using TopHat (version 2.0.10) software. The alignment files were assembled using Cufflinks and the StringTie software based on the location of the known transcript, and the transcripts of all samples were assembled again using Cuffmerge. ¹⁹ The differential gene expressed values of each sample were normalized using the fragments per kilobase of transcript per million fragments mapped (FPKM). We calculated the differential expression levels and evaluated the statistical significance of detected genes. Genes were considered as significantly differentially expressed with false discovery rates (FDRs) ≤0.05 and |log2 fold-changes| ≥1.0.

Identification of lncRNA

According to the definition of lncRNA, the unannotated reads with lengths more than 200 nt were selected as the object for further analysis by Cuffcompare from the assembled transcripts using the Cufflinks and StringTie software. The selected transcripts were also evaluated using Cuffcompare to compare with the transcripts in the ENSEMBL database. The identified transcripts were classified in different categories. These categories were considered as preliminary candidate lncRNAs, including the new intergenic transcripts, the transcripts that were exactly matched with a reference intron, and the antisense transcripts that overlapped with a reference exon or intron. To avoid the transcripts that have the ability to encode the protein, we used 4 different software programs to analyze the data, including HMMer+Pfam (based on the Pfam protein domain sequence characteristic sequence), CPC (based on Blastx comparison result and support vector machine), CNCI codon frequencies, and CPAT (based on nonaligned statistical characteristics of the sequence) with the intersection as the predicted lncRNA. The analysis of DELs was described in detail for the differential expression analysis of genes.

Functional Characterization of the Identified lncRNAs

To better understand the functional ramifications of the identified lncRNAs, we predicted the lncRNA targets from the view of bioinformatics. Differentially expressed lncRNAs were selected for target prediction. The targets of DELs were identified via cis- or trans-regulatory effects that used 2 independent algorithms. First, according to the classification and annotation information of lncRNAs, the neighbor protein genes of sense lncRNA, antisense lncRNA, intronic lncRNA, enhancer lncRNA, and bidirectional lncRNA were considered as cis-target genes. The second algorithm is based on the complementarity of the 3′UTR region sequences of the coding protein genes and ALU repeat sequences of lncRNAs. We used the Repeat Masker software to search for ALU repeat sequences, BLAST software for sequence alignment of complementary regions, and RNAplex software for prediction of the thermodynamic stability and binding ability of lncRNA and 3′UTR complexes to choose trans-acting target genes. ²⁰ Furthermore, we integrated the predicted potential lncRNA targets with the differently expressed mRNAs. The lncRNA and predicted target correlation networks showed that the lncRNAs regulated the expression of targets mRNAs.

Functional Enrichment Analysis

We used Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/), which leveraged the GO to determine the most functional annotation and classification of significant differently expressed genes and DEL targets. To demonstrate GO or molecular pathway enrichment, DAVID calculates a modified Fisher exact P value. Values <.05 are considered to be strongly enriched in the annotation category. In addition, we used the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.ad.jp/kegg/) to analyze the roles of differently expressed genes and DEL targets in the pathways.

Statistical Analysis

Transcript expression was calculated as FPKM. Known differentially expressed transcripts were screened for using the Cuffdiff program. The amount of expression was calculated using the Cuffnorm program. A transcript was determined to be differentially expressed if the difference in expression levels had a P < .05 and those showing a ≥2-fold difference in the expression between the 2 groups were screened.

Results

RNA-seq Data Quality Test and Alignment Analysis

According to the criteria for RNA-seq data acquisition mentioned in “Materials and Methods” section, clean reads were obtained after removal of the base with quality less than Q25, the adaptor of the 3′ end and sequences with a length of less than 50 bp by Trim Galore. As shown in Table 1, the total numbers of clean reads acquired from control and irradiated esophageal tissues were 149 880 052 (98.1%) and 135 995 954 (98.0%), respectively. Also, parameters such as the average length and GC content of the reads from FASTQ indicating the quality of the sequencing clean data were great, and this was vital for further analysis. In addition, we aligned the sequencing reads to a rat genome reference using TopHat. In total, at least 1 end of 130 258 395 (86.9%) and 116 763 816 (85.9%) reads for the control and irradiated esophageal tissues, respectively, was successfully mapped back to the reference genome, showing a good utilization of the sequencing reads as well as excellent credibility of the results in the downstream analysis.

Table 1.

Statistical Results of the RNA-seq Data Quality Test and Alignment Analysis of the Irradiated and Nonirradiated Esophageal Tissue Samples.

Sample	RNA-seq Data				Alignment Data
	Raw Reads	Clean Reads	Average Length	GC Content (%)	Total Mapped	Multiple Mapped	Unique Mapped
Control	152 826 434	149 880 052 (98.1%)	140.827	47	130 258 395 (86.9%)	12 750 968 (8.5%)	117 507 427 (78.4%)
Irradiated rat	138 792 138	135 995 954 (98.0%)	137.681	46	116 763 816 (85.9%)	13 197 512 (9.7%)	103 566 304 (76.2%)

Abbreviation: GC, guanine-cytosine.

Differentially Expressed Genes

In the present study, the transcript expression levels were standardized using the FPKM values with the Cufflink software. The results showed that the normalized FPKM densities for all genes were similar between control and irradiation-treated esophageal tissues (Figure 1A). The FPKM value of transcripts in irradiated samples was slightly lower than controls. After that, we used Cuffdiff to evaluate the differential expression at both the gene and the transcript levels between the control and the irradiation-treated groups using the fold-change of expression as well as the statistical significance of the values. The volcano plots (Figure 1B) show the differential expression of the transcript in 2 samples.

Figure 1.

Overview of the RNA-seq data analysis. A, Box plot of log2 (FPKM) values across control and IR expressed transcripts. Control represents normal esophageal tissue sample, and IR represents ionizing radiation-treated sample. B, Volcano plot of 2 samples (control and IR) of different genes. C, Venn diagram of the significantly differentially expressed genes (lncRNAs and mRNAs) in the 2 samples (control and IR). FPKM indicates fragments per kilobase of transcript per million fragments mapped; lncRNAs, long non-coding RNAs; mRNA, messenger RNA.

With the cutoff of FDRs ≤0.05 and |log2 fold-changes| ≥1.0, we identified 899 differentially expressed genes between the control and the irradiation-treated samples, including 398 upregulated and 501 downregulated genes (Figure 1C). The log2 (fold-changes) of the differentially expressed genes in irradiation-treated samples compared to controls ranged from −13.82 to 15.32.

Identification of lncRNAs

In order to identify the lncRNAs from the sequencing reads, as described previously, we have 3 criteria in the methods. First, we obtained 25 486 unannotated sequences and then we choose1692 primary lncRNAs with multiexons and lengths greater than 200 nt. Finally, after analysis using HMMer+Pfam, CPC, CPAT, and CNCI, we identified 717 lncRNAs in this study (Figure 2A). The droplet plots (Figure 2B) show that the lncRNAs displayed lower expression levels than identified mRNAs. The average values of FPKM of the identified lncRNAs were 2.37. In order to study the differences in lncRNA expression in irradiated esophageal tissues compared to normal ones, we also compared the expression levels of those lncRNAs in both the control and the irradiation-treated set. We found 46 DELs, 14 transcripts of which matched with the known lncRNAs from ENSEMBL. The remaining 32 were identified as novel lncRNAs. There were 14 upregulated and 32 downregulated lncRNAs in those 46 DELs (Figures 2C and 3A and Supplementary Table S1). Additionally, we found 853 differentially expressed mRNAs, of which 384 upregulated and 469 downregulated mRNAs within those differentially expressed mRNAs (Figures 2D and 3B and Supplementary Table S2).

Figure 2.

Identification of lncRNAs and the distribution of differentially expressed transcripts. A, Venn diagram of the predicted lncRNAs using HMMer+Pfam, CPC, CPAT, and CNCI. B Violin plot of log2 (FPKM) values across lncRNA and mRNA expressed transcripts. C, Venn diagram of the mRNAs in 2 samples (control and IR). D, Venn diagram of lncRNAs in 2 samples (control and IR). LncRNA indicates long non-coding RNA, and mRNA, messenger RNA.

Figure 3.

A heat map of transcriptomic expression patterns. A, A heat map of the significant mRNA expression patterns in control and IR group populations. B, A heat map of the significant lncRNA expression patterns in control and IR group populations. LncRNA indicates long non-coding RNA, and mRNA, messenger RNA.

Functional Annotation of the DELs

To investigate the possible functions of the lncRNAs, we predicted the potential targets of lncRNAs in cis-regulatory relationships. We found 2 significantly DEL transcripts with their predicted cis-regulated protein-coding genes through accurate genomic mapping and using the previously established criteria. ²¹ As shown in Figure 4A, the lncRNAs TCONS_00062404 and TCONS_00012350 are predicted to cis-regulate the genes Serpina9 and Rnps1, respectively. They were associated with regulation of the mRNA processing and mRNA surveillance pathways (Supplementary Figure S1A and B).

Figure 4.

The network between lncRNAs and their target genes. A, The coexpression network between cis-lncRNAs (TCONS_00062404 and TCONS_00012350) and their target genes. B, The coexpression network between trans-lncRNAs (TCONS_00076064, TCONS_00008311, and TCONS_00012350) and their target genes. The big box indicates novel lncRNAs, and the small box indicates the target genes. The red box indicates upregulated novel genes, and the green box indicates upregulated genes. LncRNA indicates long non-coding RNA, and mRNA, messenger RNA.

We also predicted the potential targets of lncRNAs in trans-regulatory relationships using coexpression analysis. As shown in Figure 4B, 40 interaction relationships were detected between 3 DEL transcripts (TCONS_00076064, TCONS_00008311, and TCONS_00012350) and 34 protein-coding genes in the rat reference genome. Functional analysis showed that the coexpressed genes were enriched in 111 GO terms (74 GO under biological process, 6 GO under cellular component, and 31 GO under molecular function) that encompassed a variety of biological processes (Supplementary Figure S3A). Importantly, some of the terms were nucleolar fragmentation-related terms. In addition, the coexpressed genes were enriched in 3 KEGG pathways, several of which were related to tyrosine metabolism (Supplementary Figure S2A and B).

The position of the genome, the cis-DELs (blue), and target genes (red) annotate the genomic location in the same position. The position of the genome of trans-lncRNAs and their target genes was distributed in different chromosomes (Figure 5).

Figure 5.

DEG distribution in rat chromosomes. Black bars above the horizontal line represent all DEGLs; red and blue bars represent upregulated and downregulated DEGs, respectively. DEG indicates differentially expressed gene.

Functional Enrichment Analysis of Differentially Expressed mRNAs

To analyze our RNA-seq data based on groups of functionally related genes instead of individual genes, we used DAVID to identify significantly enriched GO terms and characterized the radiation responses of esophageal tissues. It includes 3 ontology categories: biological process, molecular function, and cellular component. In our study, we obtained 1608 subcategories under biological process, 71 subcategories under cellular component, and 265 subcategories under molecular function that are significantly enriched GO terms with P values <.01. In the biological process ontology, the results showed that the genes that were affected by x-ray irradiation mainly involved the following biological processes: response to wounding, lipid metabolic process, lipid biosynthetic process, response to oxygen-containing compound, tissue development, response to alcohol, cellular lipid metabolic process, response to organic substance, cell proliferation, and response to stress. Cellular component analysis showed that the genes that were affected by x-ray irradiation were mainly located in the following groups: extracellular region, extracellular region part, extracellular space, extracellular vesicular exosome, extracellular organelle, extracellular membrane-bounded organelle, membrane-bounded vesicle, vesicle, cell surface, and extracellular matrix (ECM). The molecular function analysis showed that the functions of the genes that were affected by x-ray irradiation mainly included the following functions: catalytic activity, oxidoreductase activity, peptidase activity, acting on l-amino acid peptides, peptidase activity, oxidoreductase activity, acting on the CH–OH group of donors, NAD or NADP as acceptor, endopeptidase activity, hydrolase activity, calcium ion binding, cadmium ion binding, and carboxylic acid binding (Table 2 and Supplementary Figure S3A).

Table 2.

GO Category Enrichment of the Differentially Expressed mRNAs From Rat Esophageal Tissues After Treatment With IR (P value, .01).

GO_term	ID	P Value	Num_of_symbols in List in GO
Biological process
Response to wounding	GO:0009611	1.15E−19	93
Lipid biosynthetic process	GO:0008610	4.97E−16	98
Response to oxygen-containing compound	GO:1901700	9.46E−16	57
Cell proliferation	GO:0008283	7.07E−14	110
Response to stress	GO:0006950	1.07E−13	122
Inflammatory response	GO:0006954	2.47E−12	57
Positive regulation of cellular component movement	GO:0051272	2.67E−12	78
Regulation of cell migration	GO:0030334	4.74E−12	165
Cell adhesion	GO:0007155	5.80E−09	126
Vasculature development	GO:0001944	6.16E−08	185
Cellular component
Extracellular region	GO:0005576	1.02E−30	261
Extracellular region part	GO:0044421	3.63E−29	239
Extracellular space	GO:0005615	4.25E−24	116
Extracellular vesicular exosome	GO:0070062	2.90E−13	152
Extracellular organelle	GO:0043230	3.27E−13	152
Extracellular membrane-bounded organelle	GO:0065010	3.27E−13	152
Membrane-bounded vesicle	GO:0031988	1.05E−12	178
Vesicle	GO:0031982	1.06E−12	183
Cell surface	GO:0009986	2.08E−08	55
Extracellular matrix	GO:0031012	4.94E−07	36
Molecular function
Catalytic activity	GO:0003824	2.12E−10	325
Oxidoreductase activity	GO:0016491	8.21E−08	59
Peptidase activity, acting on l-amino acid peptides	GO:0070011	2.48E−07	66
Peptidase activity	GO:0008233	3.34E−07	68
Oxidoreductase activity, acting on the CH–OH group of donors, NAD or NADP as acceptor	GO:0016616	2.01E−06	16
Endopeptidase activity	GO:0004175	2.27E−06	56
Hydrolase activity	GO:0016787	2.89E−06	169
Calcium ion binding	GO:0005509	5.61E−06	46
Cadmium ion binding	GO:0046870	1.44E−05	4
Carboxylic acid binding	GO:0031406	1.66E−05	21

Abbreviations: GO, gene ontology; IR, ionizing radiation; mRNA, messenger RNA.

In addition, the underlying functions of differentially expressed genes were predicted using the KEGG pathway enrichment analysis. As shown in Table 3 and Figure S3B, we found10 pathways that were significantly enriched, including steroid biosynthesis, TNF signaling pathway, focal adhesion, pathways in cancer, ECM–receptor interaction, PI3K-Akt signaling pathway, cytokine–cytokine receptor interaction, arachidonic acid metabolism, chemical carcinogenesis, and steroid hormone biosynthesis.

Table 3.

KEGG Pathway Enrichment of Differentially Expressed mRNAs From Rat Esophageal Tissues After Treatment With IR (Top 10).

Pathway_name	ID	P Value	Symbols_in_list
Steroid biosynthesis	rno00100	1.56E−10	Dhcr7, Sqle, Soat1, Sc5d, Hsd17b7, Cyp24a1, Msmo1, Lss, Cyp51, Nsdhl, Ebp
TNF signaling pathway	rno04668	1.96E−06	Csf1, Traf1, Fos, Pik3cb, Ccl2, Mmp14, Icam1, Junb, Ripk3, Sele, Tnfrsf1b, Atf6b, Map2k6, Ptgs2, Socs3, Casp3, Cxcl1, Edn1
Focal adhesion	rno04510	3.07E−05	Igf1, Pdgfc, Lama3, Vav3, Pik3cb, Igf1r, Tnc, Col6a6, Lamb3, Thbs1, Spp1, Col4a6, Reln, Actn1, Flt4, Myl9, Itga3, Col6a1, Lamc2, Itgb4, Hgf, Erbb2, Col6a5
Pathways in cancer	rno05200	1.17E−02	Igf1, Lama3, Traf1, Fzd4, Fos, Pik3cb, Igf1r, Cblb, Fgfr3, Lamb3, Arnt2, Tgfb1, Ednra, Runx1, Col4a6, F2r, Zbtb16, Itga3, Ptgs2, Ednrb, Lamc2, Fgf10, Casp3, Rasgrp1, Hgf, Erbb2, Wnt5b
ECM–receptor interaction	rno04512	3.65E−05	Lama3, Tnc, Col6a6, Lamb3, Thbs1, Spp1, Col4a6, Reln, Itga3, Col6a1, Lamc2, Itgb4, Col6a5
PI3K-Akt signaling pathway	rno04151	1.81E−04	Epha2, Igf1, Csf1, Pdgfc, Lama3, Pik3cb, Igf1r, Tlr2, Tnc, Col6a6, Fgfr3, Lamb3, Thbs1, Spp1, Col4a6, F2r, Atf6b, Reln, Flt4, Itga3, Col6a1, Osmr, Lamc2, Fgf10, Chrm2, Eif4ebp1, Itgb4, Hgf, Col6a5
Cytokine–cytokine receptor interaction	rno04060	3.02E−04	Tslp, Cxcl9, Csf1, Ccl2, Cxcr2, Ccl11, Il1r1, Tgfb1, Il1r2, Tnfrsf1b, Il20rb, Tnfrsf19, Ccl21, Flt4, Csf2rb, Ccl7, Osmr, Il20ra, Tnfsf8, Tnfrsf12a, Hgf, Cxcl14
Arachidonic acid metabolism	rno00590	3.52E−04	Pla2g2f, Ptgs1, Cyp2b21, Pla2g3, Ephx2, Ptgs2, Alox12b, Alox12, Pla2g2a, Cyp2e1, Cbr3
Chemical carcinogenesis	rno05204	1.64E−03	Mgst1, Cyp1a1, Adh7, Cyp2b21, Ptgs2, RGD1562107, Ephx1, Aldh3a1, Cyp2e1
Steroid hormone biosynthesis	rno00140	3.54E−03	Srd5a1, Cyp1a1, Hsd17b7, Hsd17b2, Cyp2b21, Sts, Hsd11b2, Cyp2e1

Abbreviations: ECM, extracellular matrix; IR, ionizing radiation; KEGG, Kyoto Encyclopedia of Genes and Genomes; mRNA, messenger RNA; TNF, tumor necrosis factor.

Discussion

Radiation-induced esophageal injury has been found after mediastinal radiation in the treatment of lung and esophageal cancers. Radiation therapy is a risk factor for esophageal malignancy and hence an initiator of this carcinogenic sequence. Furthermore, most patients with ESCC respond poorly to radiotherapy; thus, it is necessary to figure out biomarkers for radiotherapy sensitivity or resistance to perform individualized therapy.

The RNA-seq sequencing method and the improvement in experimental technologies both facilitate obtaining more accurate results. The RNA-seq sequencing method correlates well with the microarray method and can detect more genes when compared with microarray data. ²² Also, oligo-dT primed RNA instead of oligo-dT primed cDNA was fragmented to avoid sequencing bias in the 3′ end of the transcript.

Using RNA-seq, we investigated radiation-induced whole transcriptome alterations in irradiated and nonirradiated esophageal tissues. Comparing the transcriptome profiles in response to IR, we identified 899 differentially expressed genes between the control and the irradiation-treated samples, including 398 upregulated and 501 downregulated genes. Also, we investigated the expression landscape across whole chromosomes. We found that the genes showing IR-altered expression were evenly located on every chromosome except sexual chromosome Y, and chromosomes 1 and 2 revealed the highest gene density (Figure 5). The functional annotation of IR-altered expression genes shows that the significantly upregulated genes included the inflammatory and immune response (Ccl11, Ccl21, Ccl2, Ccl7, Cxcl9, Csf1, Cxcl14, Cxcl1, Il36a, Tnfsf8, Ptgs1, Ptgs2, Tgfb1, Ifitm1, Ifi27l2b, Socs3, Cd14, Icam1, Il1r1, Tlr2, Tlr7, Tnfrsf12a, and Tnfrsf1b), cell growth and proliferation (Igf1, Igf1r, Igfbp4, Igfbp6, Hbegf, Fgf10, Hgf, Pdgfc, Sat1, Csf1, Nrg1, and Bcl3), and cell apoptosis (Casp3, Btg2, Epha2, Phlda1, and Plk2). Aberrant expression of some of these genes was previously reported to be important in the development of IR-induced ESCC. This evidence suggested that esophageal tissue that was exposed to IR shows a heavy inflammatory response and slight oncogenicity that has also been described in previous studies concerning irradiated non-small cell lung cancer A549 cells. ²³ However, Fgfr3, Emp2, Traf1, Frzb, Bdnf, Stk26, Tmem109, Bnipl, Bcl2l14, Bcl2l15, Btc, Traf1, and Hgf, which are involved in tumorigenesis and development progress, were downregulated after IR treatment, indicating that those cancer-associated genes might be highly sensitivity to IR in the esophageal tissue. In addition, the annotation of the RNA-seq data also revealed IR-altered expression genes implicated in epigenetic regulation, including lncRNA regulation, which is described in detail below, defined as genetic control through factors other than DNA sequence. ²⁴ Studies of epigenetic regulation to potentiate esophageal squamous cell carcinoma have recently emerged. ²⁵ In our study, we found that histone deacetylase 11, which is responsible for the deacetylation of lysine residues on the N-terminal part of core histones, was significantly downregulated in IR-stimulated esophageal tissues, suggesting that histone acetylation might be involved in response to IR.

Next, the functional enrichment analysis of the differentially expressed genes was performed by GO and KEGG pathway analysis. We found that the TNF signaling pathway, ECM–receptor interaction, PI3K-Akt signaling pathway, inflammatory processes, focal adhesion, cell proliferation, and migration could be meaningfully related to the regulation of IR treatment in esophageal tissues. In our previous study, we analyzed microRNAs (miRNAs) and circular RNAs (circRNAs) transcription alteration in radiated and nonradiated esophageal tissues. The differential expressed miRNAs and circRNAs are involved in many cellular processes, such as cell proliferation, cell migration, lipid metabolism, cellular macromolecule metabolic processes, ion binding, enzyme binding, nucleotide binding, and cellular component. ²⁶ Complemented with our previous research, our data provided new clues for the understanding of the molecular mechanisms of radiation-induced esophageal injury pathogenesis.

The lncRNAs play important functional roles in carcinogenesis by forming networks of ribonucleoprotein complexes with chromatin regulators ²⁷ and targeting their action to appropriate genomic regions, both in cis and in trans. ²⁸ Recently, the discovery and analysis of noncoding RNAs have been enhanced by RNA-seq technology. We discovered 46 DELs between control samples and irradiated samples, of which 32 were downregulated and 14 were upregulated. The target genes of cis-lncRNAs were Serpina9 and Rnps1. SART3 is involved in the regulation of mRNA splicing, probably via its complex formation with RNPS1 in cytotoxic T lymphocytes. ²⁹ Serpina9 expression could be used as a potential diagnostic marker in the differential diagnosis between follicular lymphoma and marginal zone lymphoma or mantle cell lymphoma. ³⁰ Furthermore, CCL28 is upregulated by tumor hypoxia, which promotes tumor tolerance and angiogenesis via cell migration. ^31,32 Radiotherapy for non-small cell lung cancer induces the DNA damage response and is associated with changes in the plasma levels of Mip cytokines. ³³

In conclusion, we identified a substantial number of radiation-induced–specific mRNAs and lncRNAs during radiation-induced esophageal injury, and we have imputed potential inflammatory and immunological functions for them in the pathogenesis of this disease. Moreover, our results provide interesting potential clues into the mechanisms of DNA damage response gene regulation. Our studies demonstrated that mRNA, miRNA, lncRNA, and circRNA were present differently expressed after radiation and may be served as a novel potential target for the protection of normal esophageal tissue while exposed to radiation therapy. As the roles of lncRNAs in other ESCC have not yet been fully identified and understood, this analysis should provide a valuable resource of information for future studies.

Supplemental Material

supplementary_materials - An RNA-seq-Based Expression Profiling of Radiation-Induced Esophageal Injury in a Rat Model

supplementary_materials for An RNA-seq-Based Expression Profiling of Radiation-Induced Esophageal Injury in a Rat Model by Zhiqiang Sun, Jinhui Li, Min Lin, Shuyu Zhang, Judong Luo and Yiting Tang in Dose-Response