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.

   Keywords: RNA-sequence, mRNA, lncRNAs, ionizing radiation, esophageal
   injury

Introduction

   Lung and esophageal cancers are the most common cancer.^[35]1,[36]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%.^[37]3,[38]4 Treatment includes surgery, adjuvant
   chemotherapy, and radiation therapy; however, chemoradiotherapy is
   burdened by many adverse effects.^[39]5 Radiation-induced acute injury
   is the main reason for the suspension of radiotherapy and the
   unsuccessful treatment of cancer.^[40]6,[41]7 Esophageal epithelial
   cells are extremely sensitive to ionizing radiation (IR). Different
   doses of IR may cause different degrees of esophageal epithelium
   damage.^[42]8 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.^[43]9-[44]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.^[45]12 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^2+ and mPTP), and components of the
   protein kinase and ubiquitination/proteasome system.^[46]13-[47]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.^[48]16 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.^[49]17 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.^[50]18 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
   ([51]http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/).
   Trimmed sequence files underwent quality control analysis using FastQC
   ([52]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
   ([53]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.^[54]19 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.^[55]20 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; [56]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 ([57]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 [58]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%)
   [59]Open in a new tab

   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
   ([60]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 ([61]Figure 1B) show the
   differential expression of the transcript in 2 samples.

Figure 1.

   [62]Figure 1.
   [63]Open in a new tab

   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 ([64]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 ([65]Figure 2A). The droplet plots ([66]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 ([67]Figures 2C and [68]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 ([69]Figures 2D and [70]3B
   and Supplementary Table S2).

Figure 2.

   [71]Figure 2.
   [72]Open in a new tab

   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.

   [73]Figure 3.
   [74]Open in a new tab

   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.^[75]21 As shown in [76]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.

   [77]Figure 4.
   [78]Open in a new tab

   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 [79]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 ([80]Figure 5).

Figure 5.

   [81]Figure 5.
   [82]Open in a new tab

   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 ([83]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
   [84]Open in a new tab

   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
   [85]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
   [86]Open in a new tab

   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.^[87]22
   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
   ([88]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.^[89]23
   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.^[90]24 Studies of epigenetic
   regulation to potentiate esophageal squamous cell carcinoma have
   recently emerged.^[91]25 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.^[92]26 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^[93]27 and targeting their action to appropriate genomic
   regions, both in cis and in trans.^[94]28 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.^[95]29 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.^[96]30 Furthermore, CCL28 is upregulated by tumor hypoxia,
   which promotes tumor tolerance and angiogenesis via cell
   migration.^[97]31,[98]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.^[99]33

   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
   [100]Click here for additional data file.^ (1.8MB, rar)

   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

Footnotes

   Authors’ Note: Zhiqiang Sun and Jinhui Li contributed equally to this
   work.

   Declaration of Conflicting Interests: The author(s) declared no
   potential conflicts of interest with respect to the research,
   authorship, and/or publication of this article.

   Funding: The author(s) disclosed receipt of the following financial
   support for the research, authorship, and/or publication of this
   article: This work is supported by the National Natural Science
   Foundation of China (81773224), Social Development Projects of Jiangsu
   Province (BE2018643), 333 project of Jiangsu scientific research
   program (BRA2018168), and Scientific Program of Changzhou (CJ20180029;
   CJ20179032; CE20175025).

   Supplemental Material: Supplemental material for this article is
   available online.

References