Whole Transcriptome-Based ceRNA Regulatory Network Analysis of Radiation-Induced Esophageal Epithelial Cell Injury

Introduction

Radiation-induced esophageal injury (RIEI) is a significant complication observed during radiation therapy for esophageal tumors or tumors located in close proximity to the esophagus, such as lung cancer.¹ Clinically, RIEI is characterized by esophageal inflammation, mucosal congestion, edema, difficulties in swallowing, and has been associated with the development of esophageal cancer.^2–4 Concurrent chemoradiation therapy (CCRT) enhances survival rates in patients with locally advanced non-small cell lung cancer, but it increases toxicity, particularly radiation esophagitis (RE). Among 1082 patients who underwent CCRT, the development of RE was prevalent, with grade 2 observed in 348 patients (32.2%), grade 3 in 185 (17.1%), and grade 4 in 10 (0.9%). There were no deaths related to RE.⁵ The damage from the ionizing radiation to biological targets may occur due to direct, indirect, and bystander effects. Direct effects occurs due to targeted assault on intracellular biomolecules such as DNA bases, as a result of damage, double-strand breaks, and alterations in protein structure.⁶ Indirect effects occur due to the production of reactive oxygen species (ROS) and free radicals, leading to oxidative stress and subsequent cellular damage at molecular-level.⁷ Bystander effects occur when neighboring normal cells respond to signals from damaged cells, inducing toxic effects on adjacent healthy tissues.⁸ Following radiation exposure, these three above-mentioned mechanisms collectively contribute to the onset of RIEI; however, the precise molecular mechanisms underlying them remain to be fully elucidated. During the acute inflammatory response of RIEI, there is an upregulation of free radicals (ROS), interferon-alpha (IFN-α), interleukin-6 (IL-6), tumor necrosis factor-alpha (TNF-α), interferon gamma-inducing factor (IFNGR), and interferon gamma (IFN-γ), contributing to pathogenic damage.^9–13 Additionally, our previous research has revealed the roles of cyclooxygenase-2 (COX-2) – mediated arachidonic acid metabolism and IL-17 signaling-mediated inflammation in RIEI.^14,15 Nevertheless, additional research is required to elucidate the exact molecular mechanisms responsible for RIEI.

Competing endogenous RNA (ceRNA) is a relatively newly discovered RNA regulatory mechanism, where non-coding RNAs (such as lncRNA and circRNA) compete with mRNA for binding to miRNA, thereby modulating mRNA function at the post-transcriptional level.^16,17 Studies conducted over recent years have implicated ceRNA networks in the pathogenesis of diverse diseases, such as Parkinson’s disease,¹⁸ Alzheimer’s disease,¹⁹ intervertebral disc degeneration,²⁰ depression,²¹ and esophageal cancer.^22,23 Several studies in the past few years, have unraveled abnormal RNA expression profiles in the RIEI animal model. Sun et al investigated the relationship between differentially expressed 36 lncRNAs and 853 mRNAs, using RNA-seq data obtained from both normal and radiation-treated rat esophageal tissues.²⁴ Luo et al examined the association between miRNAs and circRNAs in a rat RIEI model.²⁵ Likewise, in our previous study, we investigated an immune-related ceRNA network involving lncRNAs, miRNAs, and mRNAs in the esophageal tissues of radiation-exposed rats.²⁶ Their findings unveiled potential regulatory associations and provided insights into the pathogenesis of RIEI. Nonetheless, these aforementioned studies focused primarily on investigating tissue-level mechanisms underlying disease progression in a rat model, irrespective of interspecies differences between rats and humans. Additionally, they did not explore potential relevant cellular-level mechanisms specific to esophageal tissues. Liu et al found that enhanced expression of miR-132-3p intensifies radiation-induced damage in esophageal epithelial cells.²⁷ Sun et al examined the role of terminal uridylation enzyme 4 (TUT4) in modifying miRNA sequences and demonstrated the impact of TUT4 deficiency on downstream miRNA uridylation, thereby influencing damage to radiation-induced esophageal epithelial cells.²⁸ Their cellular-level studies focused solely on the impacts of miRNAs and their modifications on injury caused by radiation-induced esophageal epithelial cells; however, they did not incorporate lncRNAs, circRNAs, miRNAs, and mRNAs to investigate potential underlying mechanisms.

In this study, we performed phenotypic experiments on human esophageal epithelial cells following different radiation exposure and subsequently extracted RNA for whole transcriptome sequencing (lncRNA, miRNA, mRNA, and circRNA). Through whole transcriptome analysis, we identified differentially expressed RNAs depending on radiation dosage. Furthermore, we constructed a ceRNA network involving lncRNAs, miRNAs, and mRNAs in radiation-induced injury of esophageal epithelial cells which may contribute to the development of RIEI in humans. These findings may provide and information for understanding the pathogenesis and novel insights for developing treatment strategies for RIEI.

Material and Methods

Cell Culture and Irradiation

Human esophageal epithelial cells (Het-1A) were purchased from ATCC (Catalog number: CRL-2692), the cells were cultured in Dulbecco’s modified Eagle medium supplemented with 10% fetal bovine serum (VivaCell, Shanghai, China) and 1% penicillin-streptomycin (VivaCell, Shanghai, China) under at 37 °C and 5% CO₂. Based on the evaluation of the morphology and density of Het-1A cells at various doses and time points (Supplementary Figure 1), the Het-1A cells were ultimately irradiated with 0, 2, and 4 Gy using an X-ray linear accelerator (KUBTEC XCELL 320, Milford, CT) at a fixed dose rate of 1.7 Gy/min.²⁹

Cell Apoptosis Assay

Het-1A cells were seeded in 6-well plates in triplicate, cultured overnight, and subsequently irradiated. After 6, 24 and 48 h of radiation, the cells were stained with fluorescein isothiocyanate-conjugated Annexin V and propidium iodide (AV/PI; Yeason, Shanghai, China). Flow cytometry analysis was performed using the FACS Celesta flow cytometer (Becton Dickinson, Franklin Lakes, NJ), and data were analyzed using FlowJoTM software (Version 10.7). In the AV/PI co-staining experiment, viable cells showed low fluorescence intensity. The apoptosis rate was estimated as the sum of the early apoptosis rate (lower-right quadrant) and late apoptosis rate (upper-right quadrant).

Lactate Dehydrogenase (LDH) Release Assay

The toxicity in Het-1A cells following exposure to various radiation doses was evaluated by quantifying the release of LDH into the supernatant using an LDH cytotoxicity assay kit (Beyotime, Nantong, China) following the manufacturer’s instructions. The results were reported as the mean ± standard error of the mean of absorbance measured at 490 nm using a microplate reader (Biotek, Winooski, VT).

Cell Viability Assay

Het-1A cells were seeded in a 96-well plate at a density of 5×10³ cells/ well, followed by treatment with predetermined radiation doses. After 6, 24 and 48 h of radiation, cell viability was evaluated using the Cell Counting Kit-8 assay kit (APExBIO). The optical density was quantified at 450 nm using a microplate reader (Biotek, Winooski, VT).

Cellular Reactive Oxygen Species Assay

Het-1A cells were seeded in triplicate within 6-well plates and cultured overnight. Subsequently, the cells were irradiated with radiation doses of 0, 2, or 4 Gy. Following 6, 24 and 48 h of irradiation, the levels of ROS in Het-1A cells were evaluated using the ROS-sensitive dye, 2’,7’-dichlorofluorescein diacetate (DCF-DA). The cells were incubated with 10 mM DCF-DA at 37°C for 20 min and washed with phosphate-buffered saline. Subsequently, ROS levels were measured using the FACS Celesta flow cytometer (Beyotime, Nantong, China), and data were analyzed using FlowJoTM software (Version 10.7).

Colony Formation Assay

For the colony formation assay, 600 cells were plated into 6-well plates and incubated for 11 days. Colonies were then fixed with methanol and stained with crystal violet solution when they reached an adequate size for visualization. Stained colonies were photographed and counted by imageJ software. Each group was tested independently in triplicate.

Cellular DNA Damage Detection

Het-1A cells were seeded into confocal dishes and allowed to adhere overnight. Subsequently, the cells were exposed to ionizing radiation at doses of 0, 2, and 4 Gy. To evaluate DNA damage, γ-H2AX expression levels were quantified using a γ-H2AX Immunofluorescence Kit (Beyotime, Nantong, China) at three time points (6, 24, and 48 hours) post-irradiation, following the manufacturer’s protocol. Fluorescence imaging and data acquisition were performed using a Zeiss LSM710 confocal microscope (ZEISS, Germany).

Whole Transcriptome Sequencing (miRNA, lncRNA, mRNA, and circRNA)

Total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific, USA) according to the manufacturer’s instructions. The quantity and purity of total RNA were estimated using a Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA). To remove rRNA, approximately 5 μg of total RNA was treated with the Ribo-Zero Gold rRNA Removal Kit (Illumina, San Diego, CA, USA). The remaining RNA was fragmented of shorter lengths using the NEBNext^® Magnesium RNA Fragmentation Module (New England BioLabs). The steps for library preparation and sequencing of the four types of RNA are as follows:

Sequencing of mRNA and lncRNA

RNA fragments were reverse transcribed into cDNA employing SuperScript™ II Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA). Subsequently, the second strand of DNA with U-label was synthesized. Next, the U-labeled second strand DNA underwent removal treatment, followed by PCR amplification of the resulting products. Finally, paired-end sequencing (2×150 bp) was performed on the Illumina Novaseq™ 6000 platform. The generated reads were filtered and trimmed using Cutadapt software (version: cutadapt-1.9). The filtered reads were assessed for quality using the FastQC tool (version 0.11.9).

Reads from all samples were aligned to the reference genome using HISAT2 software (version: hisat2-2.2.1). Subsequently, the transcript was reconstructed using a combination of HISAT2 and StringTie software (version: stringtie-2.1.6). The resulting transcripts were categorized using gffcompare software (version: gffcompare-0.9.8). Transcripts overlapping with known mRNAs, lncRNAs, and those shorter than 200 bp were filtered out. Novel transcripts with coding potential were predicted using CPC0.9-r2 and CNCI2.0 tools. Transcripts with Coding Potential Calculator (CPC) scores < 0.5 and CNCI scores < 0 were retained and considered as novel lncRNAs. The levels of all transcripts were estimated using StringTie and ballgown software, and the expressions were analyzed using Fragments Per Kilobase of transcript per Million mapped reads (FP) values.

Sequencing of circRNA

First, the Cutadapt tool was used to remove reads containing adaptor contamination, low-quality bases, and undetermined bases. Circular RNA was then assembled using a combination of CIRCExplorer2 (version 2.2.6) and CIRI (version 2.0.2) with the aligned reads. Back-spliced reads among the unmapped reads were identified using Tophat-fusion, CIRCExplorer2, or CIRI. We then normalized the back-spliced reads supporting circular RNAs based on their read length and the number of mapped reads (splice-aligned reads per billion mapped, SRPBM) to facilitate quantitative comparisons between different RNA-seq experiments. The SRPBM value for each circular RNA in a sample was calculated by dividing the number of back-spliced reads by the total number of mapped reads (in billions) and the read length. SRPBM = number of circular reads / number of mapped reads (in billions) / read length.

Sequencing of miRNAs

The miRNA library was prepared after isolating total RNA from the samples and purifying using the provided protocol. The amplified products were then purified and enriched using polyacrylamide gel electrophoresis for miRNA library preparation. Single-end sequencing (SE50) was performed on the Illumina HiSeq 2500 platform through standard procedures. The raw reads were initially processed using an in-house developed tool, the ACGT101-miR (v4.2) program, to remove low-quality reads and eliminate any biases introduced during library preparation or sequencing. Subsequently, unique sequences with lengths in the range of 18 to 26 nucleotides were aligned to specific species precursors in miRBase 22.1 using BLAST search; this allowed for the identification of known miRNAs and novel 3p- and 5p-derived miRNAs. Unique sequences aligning to mature miRNA precursors of specific species were classified as known miRNAs.

Principal Component Analysis (PCA) and t-Distributed Stochastic Neighbor Embedding (t-SNE)

The “stats” and “Rtsne” R packages were employed to conduct Principal Component Analysis (PCA) and t-Distributed Stochastic Neighbor Embedding (t-SNE), respectively. These analyses were performed to assess the expression differences of long non-coding RNAs (lncRNAs), circular RNAs (circRNAs), microRNAs (miRNAs), and messenger RNAs (mRNAs) across distinct irradiation groups (0 Gy, 2 Gy, and 4 Gy).

Gene Set Enrichment Analysis (GSEA)

Enrichment analysis involved gene conversion using the “org.Hs.eg.db” and “clusterProfiler” R packages, and the GSEA method was applied. According to NES values, pathways were ranked in descending order, and the top five activated and inhibited pathways were selected for presentation.

Identification of Dose-Dependent Differentially Expressed RNAs

Differential expression was analyzed using the “DESeq2” R package to examine the variations in RNA expression between different groups for each RNA type: 2 Gy vs 0 Gy, 4 Gy vs 0 Gy, and 4 Gy vs 2 Gy. The threshold for differential conditions was set at |log2Fold Change| ≥ 1 and P-value < 0.05. The differentially expressed RNAs (DERs) included differentially expressed circRNAs (DECs), lncRNAs (DELs), miRNAs (DEMs), and mRNAs (DEGs). Following the differential analysis, we identified differentially expressed RNAs (referred to as dd-DERs) exhibiting dose-dependent upregulation or downregulation. The filtering criteria for differentially expressed RNAs (dd-DERs) were based on our previous study,²⁶ upregulated dd-DERs were defined as 4 Gy ≥ 2 Gy > 0 Gy, and downregulated dd-DERs were defined as 4 Gy ≤ 2 Gy < 0 Gy.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis

GO and KEGG pathway analyses were conducted on the differentially expressed genes employing the R package “clusterProfiler”. The GO analysis focused solely on the Biological Process (BP) category. A P-value < 0.05 indicated statistical significance.

Construction of Competing Endogenous RNA (ceRNA) Network

In section 2.9 of the research methodology, we could identify a limited number of differentially expressed circRNAs. Subsequent screening revealed no circRNAs with dose-dependent regulation. Hence, for the construction of, we focused exclusively on utilizing dose-dependent differentially expressed lncRNAs (dd-DELs), dd-DEMs, and dose-dependent differentially expressed genes (dd-DEGs) to construct the ceRNA network. The criteria for correlation significance were defined as a P-value < 0.05 and a correlation coefficient < −0.7. For this assessment, we conducted a co-expression analysis. The sequences of the co-expressed RNAs were compiled, and target gene was predicted using the OmicStudio Target Gene Prediction Cloud Analysis platform (https://www.omicstudio.cn). Predictions were generated by integrating outcomes from TargetScan (version 5.0) and miRanda (version 3.3a) softwares, employing the threshold for target selection as TargetScan_score ≥ 50 and miranda_Energy < −20, and the overlapping results from both prediction tools were considered. This approach enabled the identification of dd-DELs and dd-DEGs, potentially interacting with dd-DEMs. Subsequently, a ceRNA network was constructed by integrating co-expression analysis and predicted binding sites using shared dd-DELs, dd-DEMs, and dd-DEGs.²⁶ The ceRNA network was visualized using Cytoscape software (version 3.10.1).

Reverse Transcription Quantitative Polymerase Chain Reaction (RT-qPCR)

Total RNAs were extracted from rat esophagus using TRIzol reagent (Thermo Fisher Scientific, USA), and three biological replicates were used for each sample group. Total RNA was reverse transcribed into cDNA using the miRNAs 1st Strand cDNA Synthesis Kit (Vazyme, China) and FastKing RT Kit (TIANGEN, China). Subsequently, RT-qPCR was performed on a BioRad real-time fluorescent quantitative PCR instrument employing the miRNA Universal SYBR qPCR Master Mix (Vazyme, China) and SuperReal PreMix Plus (TIANGEN, China). The relative expression levels of RNAs were assessed using the 2^−ΔΔCt method. All primer sequences for the experiment are provided in Supplementary material 1.

Statistical Analysis

Statistical analysis and data visualization were performed using R version 4.3.1 and GraphPad version 9.5.1. The threshold for statistical significance was set at P-value < 0.05. For comparisons between two groups, independent samples t-tests were employed, while for comparisons involving multiple groups, either one-way ANOVA or Kruskal–Wallis tests were utilized.

Results

Cell Phenotype Analysis of Radiation-Induced Injury in Het-1A Cells

To systematically assess radiation-induced damage in Het-1A cells, we conducted a comprehensive series of phenotypic analyses following exposure to 0, 2, and 4 Gy radiation doses. Our observations revealed significant radiation-induced morphological alterations accompanied by a dose-dependent reduction in cell density (Figure 1A). Quantitative analyses demonstrated consistent dose-response relationships: apoptosis rates, intracellular LDH release, and ROS levels showed progressive increases, while cell viability exhibited a corresponding decline (Figure 1B–E). The clonogenic survival assay further confirmed the radiation sensitivity of Het-1A cells, with colony-forming capacity decreasing proportionally with escalating radiation doses (Figure 1F). Parallel evaluation of DNA damage through γ-H2AX foci quantification at multiple time points revealed a clear dose-dependent accumulation of DNA damage markers (Figure 1G). Collectively, these findings demonstrate that radiation doses of 2 and 4 Gy induce substantial cellular damage in Het-1A cells, with the severity of damage exhibiting a positive correlation with radiation dose.

Figure 1 Continued.

Figure 1 Phenotypic outcomes of radiation damage in Het-1A cells (A) Microscopic observation of cell morphology. (B) The apoptosis rate was detected at 6, 24, and 48 h after irradiation. (C) The release rate of lactate dehydrogenase (LDH) was detected at 6, 24, and 48 h after irradiation. (D) Cell viability was measured at 6, 24, and 48 h after irradiation. (E) Intracellular ROS levels were measured at 6, 24, and 48 h after irradiation. (F) Assay of colony formation after irradiation. (G) Images were taken by confocal microscopy with γ-H2AX. * P-value < 0.05; ** P-value < 0.01; *** P-value < 0.001; **** P-value < 0.0001. Data are presented as mean ± standard error (n = 3).

Whole-Transcriptome Expression Profiling of Radiation-Induced Esophageal Epithelial Cell Injury

Whole transcriptome sequencing was performed on cell samples at 48 h after 0, 2, and 4 Gy irradiation to evaluate the differential expression of circRNA, lncRNA, miRNA, and mRNA. PCA and t-SNE evaluated the overall expression patterns of the four RNA types in each group and showed well-defined clusters of lncRNA (Figure 2A and Supplementary Figure 2A), miRNA (Figure 2B and Supplementary Figure 2B), mRNA (Figure 2C and Supplementary Figure 2C), and circRNA (Figure 2D and Supplementary Figure 2D) across the four groups. It was shown that the overall expression of these RNA types differed significantly between the groups.

Figure 2 Principl Component Analysis plots illustrate the expression levels of lncRNA (A), miRNA (B), mRNA (C) and circRNA (D) in Het-1A cells exposed to 0, 2 and 4 Gy.

Gene Set Enrichment Analysis of mRNA Expression Profiling in Radiation-Induced Esophageal Epithelial Cell Injury

GSEA was employed to investigate functional changes at the cellular-level by analyzing the expression profiles of all mRNAs. GSEA analyses were carried out on the three pairwise comparisons, including 2 Gy vs 0 Gy, 4 Gy vs 0 Gy, and 4 Gy vs 2 Gy. The GSEA analysis utilizing the GO gene set revealed significant enrichment in the activation of channel activity and mitotic inhibition following ionizing radiation exposure. Specifically, molecular functions (MFs) such as voltage-gated potassium channel activity, potassium channel activity, and monoatomic ion channel activity were upregulated (NES > 0), whereas pathways related to BP and cellular components (CCs), including sister chromatid cohesion, kinetochore, and condensed chromosome, centromeric region were downregulated (NES < 0) across all three comparative analyses (Figure 3A). The enrichment analysis using KEGG revealed radiation activated immune-related pathways involving viral protein interaction with cytokine and cytokine receptor, neuroactive ligand-receptor interaction, cytokine-cytokine receptor interaction, complement and coagulation cascades, and calcium signaling pathway (Figure 3B). Conversely, radiation-inhibited repair-related pathways, such as nucleotide excision repair, base excision repair, and nucleocytoplasmic transport were also noted (Figure 3B). The enrichment analysis of both gene sets mentioned above revealed significant disruptions in channel activity, cell replication, cellular repair, and immune response in Het-1A cells upon radiation exposure.

Figure 3 The Gene Set Enrichment Analysis of mRNA expression profiles based on the Gene Ontology (A) and Kyoto Encyclopedia of Genes and Genomes (B) gene sets.

Identification of Differentially Expressed RNAs

DERs, including lncRNA, miRNA, mRNA, and circRNA, were analyzed in the three pairwise comparisons of radiation exposures: 2 Gy vs 0 Gy, 4 Gy vs 0 Gy, and 4 Gy vs 2 Gy. Pairwise comparisons of the three groups of lncRNAs revealed, respectively, 79 downregulated and 61 upregulated (Figure 4A and D), 91 downregulated and 86 upregulated (Figure 4B and D), and 55 downregulated and 73 upregulated DELs (Figure 4C and D). Considering miRNAs, we identified 19 downregulated and 20 upregulated DEMs under 2 Gy vs 0 Gy radiation doses (Figure 4E and H), 18 downregulated and 30 upregulated DEMs under 4 Gy vs 0 Gy (Figure 4F and H), and 9 downregulated and 6 upregulated DEMs in 4 Gy vs 2 Gy radiation doses (Figure 4G and H). Among mRNAs, we identified 182 downregulated and 161 upregulated DEGs in 2 Gy vs 0 Gy radiation doses (Figure 4I and L), 309 downregulated and 172 upregulated DEGs in 4 Gy vs 0 Gy (Figure 4J and L), and 62 downregulated and 69 upregulated DEGs under 4 Gy vs 2 Gy radiation doses (Figure 4K and L). Among circRNAs, we identified 25 downregulated and 35 upregulated DECs in 2 Gy vs 0 Gy radiation doses (Figure 4M and P), 45 downregulated and 21 upregulated DECs in 4 Gy vs 0 Gy (Figure 4N and P), and 22 downregulated and 23 upregulated DECs in 4 Gy vs 2 Gy radiation doses (Figure 4O and P). The human esophageal epithelial cells in the 4 Gy group exhibited a significantly higher number of DERs compared to those in the 2 Gy group, indicating a more pronounced disturbance at the transcriptome level. These cellular findings were in concordance with the previously documented pathological changes in the esophageal structure. The differentially expressed RNAs are shown in Supplementary Material 2.

Figure 4 Differentially expressed RNAs (DERs) were identified among the three pairwise comparisons. (A–C). Volcano plots illustrating the differential expression of lncRNAs (DELs): (A) (2 Gy vs 0 Gy), (B) (4 Gy vs 0 Gy), and (C) (4 Gy vs 2 Gy). The bar graph shows the number of DELs in each comparison (D). Volcano plots (E–G) depict the differential expression of microRNAs (DEMs): (E) (2 Gy vs 0 Gy), (F) (4 Gy vs 0 Gy), and (G) (4 Gy vs 2 Gy), while the corresponding bar graph represents the number of DEMs in each comparison (H). Volcano plots (I–K) display the differential expression of genes (DEGs) of (I) (2 Gy vs 0 Gy), (J) (4 Gy vs 0 Gy), and (K) (4 Gy vs 2 Gy). The bar graph exhibits the number of DEGs in each group (L). Volcano plots (M–O) display the differential expression of circRNAs (DECs): (M) (2 Gy vs 0 Gy), (N) (4 Gy vs 0 Gy), and (O) (4 Gy vs 2 Gy). The bar graph presents the number of DECs in each group (P).

Identification of Dose-Dependent Differentially Expressed RNAs

Utilizing the analysis performed in Section 3.4 as a foundation, We defined the criteria for radiation dose dependence as follows: either 4 Gy ≥ 2 Gy > 0 Gy or 4 Gy ≤ 2 Gy < 0 Gy. Subsequently, using these criteria, we screened and identified dd-DERs that were dose-dependent and exhibited either upregulation or downregulation. We noted 17 upregulated and 24 downregulated dd-DELs (Figure 5A), 11 upregulated and 7 downregulated dd-DEMs (Figure 5B), and 84 upregulated and 108 downregulated dd-DEGs (Figure 5C). A limited number of DECs were identified in the aforementioned differential RNA expression analysis; therefore, we did not screen out any dd-DECs. Consequently, dd-DECs were excluded from the subsequent construction of the ceRNA network. In summary, we identified a total of 41 dd-DELs, 18 dd-DEMs, and 192 dd-DEGs that potentially contribute to radiation-induced injury of esophageal epithelial cells. The Dose-Dependent Differentially Expressed RNAs are shown in Supplementary Material 3.

Figure 5 Identification of dose-dependent differentially expressed RNAs (dd-DERs). (A) Venn diagrams and heatmap of the dose-dependent differentially expressed lncRNAs (dd-DELs). (B) The Venn diagrams and heatmap of the dose-dependent differentially expressed miRNAs (dd-DEMs). (C) The Venn diagrams and heatmap of the dose-dependent differentially expressed mRNAs (dd-DEGs).

Investigation of the Potential Functions and Pathways

We performed Gene Ontology-Biological Pathway (GO-BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses to identify the potential biological functions and pathways associated with dd-DEGs in irradiated cells. According to the GO-BP analysis, the top five enriched BP terms included various metabolic processes, including hormone metabolism, cellular glucuronidation, uronic acid metabolism, glucuronate metabolism, and xenobiotic metabolism (Figure 6). The KEGG pathway analysis revealed the top five terms, including pentose and glucuronate interconversions, retinol metabolism, drug metabolism - cytochrome P450, bile secretion, and porphyrin metabolism (Figure 6). The inner circle in this figure represents the terms, while the outer circle represents the genes enriched in each term. Among the genes presented in the outer circle, nine members of the UGT1A family, including UGT1A1, UGT1A3, UGT1A4, UGT1A5, UGT1A6, UGT1A7, UGT1A8, UGT1A9, and UGT1A10, have roles in a variety of metabolic pathways. The above results suggest that dd-DEGs, especially the members of the UGT1A family, primarily exert their functions by regulating metabolic pathways, thereby providing insights into the effects of radiation on protein functions.

Figure 6 Gene Ontology-Biological Process and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses based on dose-dependent differentially expressed genes (dd-DEGs).

Construction of Competing Endogenous RNA (ceRNA) Network

The relationships among dd-DELs, dd-DEMs and dd-DEGs were investigated by constructing a ceRNA network to explore the reciprocal regulatory interactions. We able to identify a total of 113 co-expressed dd-DELs/dd-DEMs exhibiting a strong negative correlation (correlation coefficient < −0.7, P-value < 0.05) (Figure 7A). Further online predictions indicated interactions for 141 dd-DELs/dd-DEMs (Figure 7A). Additionally, correlation among 770 co-expressed dd-DEMs/dd-DEGs were observed, and online predictions suggested binding for 31,639 dd-DEMs/dd-DEGs (Figure 7B). We constructed a ceRNA network through the intersection of the co-expression and online prediction sets (Figure 7C). This network comprised 10 lncRNAs, 5 miRNAs, and 55 mRNAs. The heatmap was generated to visualize the expression levels of all RNAs extracted from the network (Figure 7D). These RNA molecules are believed to be involved in RIEI through the competitive mechanism of ceRNA, where they compete to regulate the expression of target genes. The ceRNA network offers valuable insights into the regulatory interactions between diverse RNA types and their potential roles in RIEI. The RNA relationships involved in constructing the ceRNA network are shown in Supplementary Material 4.

Figure 7 Competing endogenous RNA (ceRNA) network composed of lncRNA, miRNA, and mRNA. (A) Venn diagram illustrating the co-expression dose-dependent differentially expressed lncRNAs/ dose-dependent differentially expressed miRNAs/ (dd-DELs/dd-DEMs) and binding dd-DELs/dd-DEMs. (B) Venn diagram depicting the co-expression dd-DEMs/ dose-dependent differentially expressed genes (dd-DEMs/dd-DEGs) and binding dd-DEMs/dd-DEGs. (C) CeRNA network diagram constructed based on filtered dd-DELs, dd-DEMs, and dd-DEGs. lncRNA, miRNA, and mRNA are represented by “V” shape, circle, and triangle, respectively. (D) Heatmap of RNA expression in the ceRNA network.

Validation of Competing Endogenous RNA Network Through RT-qPCR

The accuracy of whole transcriptome sequencing was validated by selecting several lncRNAs, miRNAs, and mRNAs exhibiting high connectivity in the network for validation through reverse transcription-quantitative polymerase chain reaction (RT-qPCR). The RT-qPCR findings demonstrated an upregulation of IL11 mRNA expression with increasing radiation dose, while NRP2, SLC16A7, and KANK4 were downregulated. Moreover, the expression levels of lncRNAs (AC048341, AL353151) decreased in response to escalating radiation dose. Conversely, miRNA (hsa-mir-9902-1-p3_1ss13GC) exhibited enhanced expression, whereas the miRNA (hsa-miR-2116-3p_L+1R-1) level decreased (Figure 8). Significantly, the experimental RT-qPCR results aligned with the sequencing data, partially confirming the accuracy of our sequencing outcomes. These findings provide valuable insights into the potential mechanisms of radiation-induced damage in esophageal epithelial cells.

Figure 8 Validation of multiple dose-dependent differentially expressed RNAs (dd-DERs). (A) Summary of the validation of whole transcriptome results. (B) Validation of expression levels of multiple RNAs in the network by real time-quantitative polymerase chain reaction. *P-value < 0.05; ** P-value < 0.01; *** P-value < 0.001; **** P-value < 0.0001. Data are presented as mean ± standard error (n = 3).

Discussion

In our study, we employed whole transcriptome sequencing technology to investigate the differential expression of mRNA, lncRNA, circRNAs, and miRNA in human esophageal epithelial cells following exposure to different radiation doses. Using dose-dependent genes as an entry point, we identified key biological processes and pathways closely associated with radiation-induced esophageal epithelial cell injury. Through correlation analysis and sequence-based prediction, we successfully identified a potential ceRNA network of dose-dependent differentially expressed mRNAs, lncRNAs, and miRNAs.

We identified a ceRNA regulatory network comprising 10 lncRNAs, 5 miRNAs, and 55 mRNAs. In irradiated cells, the lncRNAs LINC01006 and SATB1-AS1 were significantly downregulated. Although some identified lncRNAs have not been previously reported to act through endogenous competition in RIEI, their roles in other diseases have been implicated. LINC01006, in particular, has been shown to influence cancer progression through endogenous competition across various cancers. In hepatocellular carcinoma (HCC), LINC01006 suppresses the development of tumors via the LINC01006/miR-194-5p/CADM1 axis.³⁰ Additionally, it regulates HCC cell proliferation, migration, and invasion through the LINC01006/miR-433-3p/CBX3 axis.³¹ In another example, the downregulation of SATB1-AS1 has been linked to enhanced migration and invasion of laryngeal squamous cell carcinoma cells through the SATB1-AS1/miR-1299/SAV1/CCNG2/SH3KBP1/JADE1/HIPK2 axis.³² The most significantly upregulated and downregulated mRNAs in the network are, respectively, RHBG and PLCXD3. Research shows that RHBG can initiate intracellular signaling cascades that ultimately activate NFκB,³³ with nuclear factor-κB activation contributing to inflammation, such as radiotherapy-induced mucositis.³⁴ Downregulated PLCXD3 is associated with impaired insulin signaling and biosynthesis, and reduced glucose sensing in pancreatic β cells.³⁵ The miRNAs involved in this endogenous competition network are variants of mature miRNAs available in the database, and their specific effects in the regulatory functions of cells or tissues remain unclear. Nonetheless, certain mechanisms associated with some miRNAs in the database have been confirmed. For example, Wu et al described a positive feedback loop in breast cancer cells (LINC01433/miR-2116-3p/MYC), whereas LINC01433 promotes cell progression by regulating the miR-2116-3p/MYC axis.³⁶ While these RNAs play a role in disease progression through endogenous competition, their roles in RIEI require further investigation.

The GSEA analysis identified the channel activity activation, including the enrichment of potassium ion and monovalent ion channel activities. It also revealed the activation of immune responses and the suppression of cell mitosis. Earlier studies have documented that ion channels are potentially promising targets for modulating diverse cellular functions. Furthermore, specific ion channels on the cell membrane are accompanied by intracellular counterparts, each serving distinct roles within various metabolic pathways.³⁷ Our phenotypic results revealed a gradual increase in the rate of cell apoptosis with escalating radiation dose in agreement with elevated levels of ROS. Notably, mitochondria serve as the main source of ROS production and enhanced ROS levels can trigger the release of cytochrome c from mitochondria through mitochondrial pathways, ultimately culminating in cell apoptosis.³⁸ Experimental evidence supports the hypothesis that ion channels play a potential regulatory role at the mitochondrial level during apoptosis. Among these channels, potassium (K+) channels are widely present in organisms and facilitate the transport of K+ across biological membranes. Several events can be triggered through the regulation of K+ channels, including K+ flux, such as mitochondrial membrane hyperpolarization, significantly impacting the generation of mitochondrial ROS.^39–41 Voltage-gated potassium channels (Kv) constitute the largest and most diverse family of potassium channels, divided into 12 families (Kv1-12).⁴² The Kv1.3 channel, a member of the Kv channels, belongs to the Shaker family.⁴³ It was initially discovered in the plasma membrane (PM),⁴⁴ where it is involved in proliferation and cell volume regulation. Studies have shown that MitoKv1.3 participates in apoptotic cell death.⁴⁵ In cells lacking Kv1.3, the expression of mitochondria-targeted Kv1.3 is sufficient to restore these cells’ response to several apoptotic stimuli.⁴⁵ MitoKv1.3 has been identified as a target of the pro-apoptotic protein Bax, a member of the Bcl-2 family.^45,46 The physical interaction between Bax and mitoKv1.3 occurs only upon the induction of apoptosis. Similar to Kv1.3, Kv1.1 and Kv1.5 also exhibit dual localization in the inner mitochondrial membrane (IMM) and PM, and can interact with Bax under apoptotic stimuli, sensitizing cells to death signals.^47,48 Radiation can induce DNA damage in cells, leading to several cellular responses. These include the removing and repairing DNA damage, restoring DNA double-strand continuity, activating DNA damage checkpoints to halt cell cycle progression for repair, and preventing the propagation of damaged or incompletely replicated chromosomes. Additionally, cells may elicit transcriptional responses that could induce beneficial changes or apoptosis to eliminate severely damaged or dysfunctional cells.⁴⁹ Our GSEA results indicate the activation of mitotic inhibition and DNA damage checkpoints within the cell. Furthermore, our enrichment analysis revealed a suppression of nucleotide excision repair (NER) and base excision repair (BER). Nucleotide excision repair removes large DNA lesions caused by radiation, chemicals, or adding proteins to DNA. The damaged bases are excised by “excision nucleases”, the multi-subunit enzyme system that creates double incisions in the damaged strand to remove the lesion.^50–53 In humans, excision repair involves six repair factors (RPA, XPA, XPC, TFIIH, XPG, and XPF•ERCC1) composed of 15 polypeptides.^54–56 Although our transcriptome data could detect expression changes in several repair factors, these changes were not statistically significant according to the defined threshold.

We identified a total of 55 dd-DEGs, with dose changes categorized as 4 Gy ≥ 2 Gy > 0 Gy or 4 Gy ≤ 2 Gy < 0 Gy. Enrichment analysis of these dd-DEGs revealed that several pathways related to glucuronic acid metabolism were enriched, and the UDP-glycosyltransferases (UGT) superfamily genes were found to be involved in multiple metabolic pathways. UGTs are crucial enzymes involved in biotransformation, regulating the levels and distribution of many endogenous signaling molecules and the metabolism of a broad range of endogenous as well as exogenous compounds. In mammals, the UGT superfamily includes the UGT1, UGT2, UGT3, and UGT8. UGTs serve as rate-limiting enzymes in the glucuronic acid pathway, and alterations in their activity can impact gluconeogenesis and lipid metabolism. UGTs conjugate various lipophilic molecules to sugars, which include galactose, xylose, acetylglucosamine, glucuronic acid, and glucose, thereby converting these molecules into water-soluble substrates, including carcinogens, drugs, steroids, lipids, fatty acids, and bile acids.⁵⁷ UGTs expression and their enzymatic activity are influenced by various factors and regulated through multiple mechanisms of epigenetic modifications (eg, DNA methylation and histone modifications), transcriptional regulation, post-transcriptional regulation (eg, microRNAs), and post-translational modifications (eg, structural and functional changes and protein-protein interactions). Methylation of the transcription factor HNF1A in colon cells negatively regulates UGT1A1.^58,59 UGT expression is also influenced by histone modifications, which together with DNA methylation, regulate UGT1A1 expression.⁶⁰ MicroRNAs can directly target UGT mRNAs to modulate their expression and indirectly affect UGT expression through the inhibition of transcription factors.^61–63 Additionally, the interaction between UGT1A9 and bisphenol leads to intracellular calcium overload through the activation of the calcium signaling pathway, a finding also supported by our analysis. This overload induces mitochondrial stress, disrupts mitochondrial homeostasis, and promotes bisphenol-induced cell death.⁶⁴ We hypothesize that members of the UGT superfamily may significantly contribute to radiation-induced esophageal epithelial cell injury, warranting further investigational studies.

Conclusion

In conclusion, our study provides significant insights into the effects of ionizing radiation on the entire transcriptome of human esophageal epithelial cells. We identified several differentially expressed RNAs (DD-RNAs) at the 48-hour time point following ionizing radiation and used these RNAs to construct an associated ceRNA network. However, the ceRNA network is a dynamic regulatory process. Our study focused on a single time point, which limits the ability to fully capture the temporal changes within the ceRNA network. Additionally, the internal regulatory relationships require further validation. Therefore, future experimental studies will aim to identify key factors and regulatory relationships that contribute to early damage and explore potential therapeutic strategies to mitigate radiation-induced damage to esophageal epithelial cells through the competitive regulatory interactions between RNAs in the ceRNA network. These findings offer potential avenues for radiation protection of the entire esophagus.