Integrative analysis for the role of long non-coding RNAs in radiation-induced mouse thymocytes responses

Abstract

Long non-coding RNAs (lncRNAs) are a critical class of regulatory molecules involved in a variety of biological functions; however, their role in immune cells response to radiation is unknown. Therefore, in this study we used integrative analysis to determine the expression profile of lncRNAs in mouse thymocytes and the potential functions of lncRNAs in response to radiation. Microarray data profiling indicated that 53 lncRNAs (36 up-regulated and 17 down-regulated) and 74 coding genes (39 up-regulated and 35 down-regulated) were highly differentially expressed in the high dose radiation (HDR) group compared with the control group. In the low dose radiation (LDR) group, only one lncRNA was down-regulated. Moreover, as compared with the control group, 109 lncRNA pathways in the HDR group and 14 lncRNA pathways in the LDR group were differentially expressed. Our data revealed the expression pattern of lncRNAs in mouse thymocytes and predicted their potential functions in response to LDR and HDR. In the HDR group, GO analysis showed that the role of lncRNAs in damage responses of thymocytes to HDR mainly involved chromatin organization and cell death. These findings might improve our understanding of the role of lncRNAs in LDR- and HDR-induced immune cells and provide a new experimental basis for further investigation.

long non-coding RNAs, radiation effects, mouse thymocytes, response mechanism

Introduction

The human body is constantly exposed to various sources of radiation, including occupational exposure, environmental exposure, and exposure during medical diagnostic or therapeutic procedures. Medical procedures, especially medical imaging, are the most significant man-made sources of radiation exposure to the general population [1]. Important advances in medical imaging coupled with extraordinary growth in its utilization have been achieved in recent years [2]. The immune system is equipped with an arsenal of strategies to combat infectious threats and maintain normal health. It is well known that the immune system responds to ionizing radiation (IR) with distinct characteristics depending on the dosage and dose rate. Medium dose of IR reduces immune function and it is a fundamental cause of radiation carcinogenesis. Low dose radiation (LDR) may enhance immune function and it is an important part of the stimulatory effect of radiation. Thus, the stimulatory effect of LDR as shown in recent years has even greater significance in understanding the health effects of environmental low-level radiation and has received considerable attention in the radiation medical field [3–5]. As such, a more comprehensive understanding of the immune system cellular and molecular responses to IR is needed to avoid adverse radiation effects on human health.

It is known that alterations in gene expression are significant components in the response to IR. Therefore, with the advancement in high-throughput genomic technologies, global gene (mRNA) expression profiling is an effective approach to obtain information regarding the global cellular response to radiation [6,7]. Recently, long non-coding RNAs (lncRNAs) have emerged as important regulators of gene expression in diverse biological contexts. LncRNAs control gene expression in the nucleus by modulating transcription or via post-transcriptional mechanisms targeting the splicing, stability or translation of mRNAs. Our knowledge of lncRNAs’ biogenesis, their cell-type specific expressions, and their versatile molecular functions are rapidly progressing in all areas of biology [8–12].

In the immune cells, lncRNAs promote the distribution of Tregs in peripheral blood T cells, which causes an enhanced cell proliferation of gastric cancer cells by recruiting TGF-β and activating the TGF-β signal pathway [13]. Here we discuss these exciting new regulators and highlight an emerging paradigm of lncRNA-mediated control of gene expression in the mouse thymocytes in response to LDR and high dose radiation (HDR). In this study, the differential expression of lncRNAs and mRNAs in mouse thymocytes was determined upon exposure to different doses of X-ray whole-body irradiation (WBI), for the purpose of exploring the potential roles of lncRNAs in the regulation of coding gene expression for cellular response and biological pathways after LDR or HDR leading to thymocyte damage. Our results highlight the functions of lncRNAs in response to radiation-induced immune cell response mechanisms.

Materials and Methods

Mice and irradiation

A total of 48 ICR mice (8 weeks old, 24 male and 24 female, weighing 18–22 g) were purchased from the Experimental Animal Center of Jilin University (Changchun, China). The mice were maintained in a 12:12 light/dark cycle at a room temperature of 22 ± 1°C with free access to food and water. The male and female ICR mice were randomly divided into three groups, including control group (sham irradiation), LDR group (WBI with X-rays at a single dose of 0.075 Gy), and HDR group (WBI with X-rays at a single dose of 4.0 Gy).

An X-ray generator (Model X-RAD320iX; Precision X-ray, Inc., North Branford, USA) was used to deliver radiation at a dose rate of 0.0134 Gy/min (170 kV; 2 mA) in the LDR group and 1.020 Gy/min (180 kV; 20 mA) in the HDR group. Thymic tissues were harvested 12 h after treatment for the subsequent experiments. These experiments were approved by the Animal Experiment Ethics Committee of Jilin University [Permission Number: SCXK(J)2013-0001].

RNA labeling and array hybridization

Sample labeling and array hybridization were performed according to the manufacturer's protocols. Briefly, total RNA was extracted from each sample using TRIzol reagent (Invitrogen, Carlsbad, USA) following the manufacturer's protocol. RNA quantity and quality were measured with a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Inc., Wilmington, USA). RNA integrity was evaluated by standard denaturing agarose gel electrophoresis. LncRNAs/mRNAs were purified from total RNA after removal of rRNA by using the Arraystar rRNA Removal Kit (Arraystar, Rockville, USA). Then, each sample was amplified and transcribed into fluorescent cRNA along the entire length of the transcripts without 3′ bias utilizing a random priming method with the Arraystar Flash RNA Labeling Kit (Arraystar). The labeled cRNAs were purified by RNeasy Mini Kit (Qiagen, Valencia, California, USA). The concentration and specific activity of the labeled cRNAs (pmol Cy3/μg cRNA) were measured with the NanoDrop ND-1000 spectrophotometer. The labeled cRNAs were hybridized onto the LncPath^TM Mouse Epigenetic Pathway Array (8 × 15 K; Arraystar) according to the manufacturer's protocol. The slides were then washed and the arrays were scanned using the Agilent Scanner G2505C (Agilent Technologies, Santa Clara, USA).

Expression profiling data and analysis

Quantile normalization of raw data and subsequent data processing were performed using the R limma software package. After quantile normalization of the raw data, low-intensity filtering was performed, and the lncRNAs/coding genes that at least three out of nine samples had flags in ‘P’ or ‘M’ (All Targets Value) were retained for further analyses. All the lncRNAs were annotated with the corresponding target genes and potential regulatory mechanisms. The gene- and transcript-level expression data of the target genes were also provided to determine the expressional relationships between lncRNAs and their potential target genes.

When comparing two groups of profile differences (such as disease versus control), the ‘fold change’ (i.e. the ratio of the group averages) between the groups for each lncRNA or coding gene was computed. The statistical significance of the difference was estimated using the Student's t-test. LncRNAs or coding genes with fold change ≥2 and P < 0.05 were considered significantly differentially expressed.

Predicted regulation of lncRNAs and coding genes

LncPath™ arrays profile lncRNAs and their target genes simultaneously in a pathway-focused or disease-specific manner. By quantifying lncRNAs (at the transcript level) and their potential target genes (at both the gene and transcript level) in parallel, LncPath™ arrays provide comprehensive insights into the expressional relationship between lncRNAs and their target genes, which is helpful for establishing rapid connections between the lncRNA regulatory mechanisms and their biological functions.

Three types of lncRNAs are included in LncPath™ arrays. (i) The lncRNAs that overlap in part with critical pathway genes or locate within 3 kb from the nearest pathway genes were defined as ‘neighboring lncRNAs’, which might exert their functional roles by regulating the neighboring pathway genes at the transcriptional or post-transcriptional level. (ii) The lncRNAs with a significant possibility of functioning as competing endogenous RNAs of the critical pathway genes were defined as ‘ceRNAs’, which share miRNA response elements with the mRNA transcripts of the corresponding target genes and may prevent these mRNAs from being degraded as post-transcriptional gene regulators. (iii) The lncRNAs that locate within 300 kb of the pathway genes and may have functional properties of enhancers were defined as ‘enhancer-like lncRNAs’, which might positively regulate the expression of the distal pathway genes at the transcriptional level.

Quantitative real-time PCR

Total RNA from the same samples as those used in the microarray was extracted using the TRIzol reagent and then reverse transcribed using the PrimeScript RT reagent Kit with gDNA Eraser (Perfect Real Time; TaKaRa, Dalian, China) according to the manufacturer's recommendations. The quantitative real-time PCR (qRT-PCR) was performed using SYBR Green assay (TaKaRa), and the mouse housekeeping ’GAPDH’ was used as an internal control. The primers for each lncRNA were designed according to Primer 3 (http://sourceforge.net/projects/primer3/) online and checked with Basic Local Alignment Search Tool (BLAST) of National Center of Biotechnology Information (NCBI) to ensure unique amplification product as listed in Table 1. PCR was performed in a total reaction volume of 25 μl, including 12.5 μl of SYBR Premix Ex Taq (2×), 1 μl of PCR Forward Primer (10 μM), 1 μl of PCR Reverse Primer (10 μM), 2 μl of cDNA, and 8.5 μl of double-distilled water. The cycling conditions consisted of an initial single cycle of 30 s at 95°C, and followed by 40 cycles of 5 s at 95°C and 30 s at 60°C. All reactions were run in triplicate. The relative expression of genes was calculated based on the 2^−ΔΔCt method.

Gene Ontology and coding–non-coding gene co-expression network

The Gene Ontology (GO) project provides a controlled vocabulary to describe gene and gene product attributes in any organism (http://www.geneontology.org). The ontology covers three domains: Biological Process, Cellular Component, and Molecular Function. Fisher's exact test is used to determine if there is more overlap between the differentially expressed (DE) list and the GO annotation list than that would be expected by chance. The P-value denotes the significance of GO term enrichment in the DE genes. The lower the P-value, the more significant the GO term (P < 0.05 is recommended). A coding–non-coding gene co-expression network (CNC network) was constructed with three pairs of the most significantly changed lncRNAs and their target genes in the HDR group, which were presented to identify interactions among genes. The method was performed according to Liao [14] and performed by KangChen Bio-tech (Shanghai, China).

Statistical analysis

The results were presented as the mean ± SD. Statistical comparisons of the experimental results between the treated and control groups were evaluated using the Student's t-test. All statistical tests were performed with SPSS version 22.0 (SPSS, Inc., Chicago, USA). P < 0.05 was considered as statistically significant difference.

Results

Profile of microarray data

For the microarray analysis, 499 lncRNAs (Fig. 1A,B) and 649 coding transcripts (Fig. 1C,D) were detected. They were carefully collected from the most authoritative databases such as RefSeq, UCSC Knowngenes, Ensembl, and many related reports. Scatter-Plot, a visualization method, was used to assess gene expression variation between the two groups and hierarchical clustering showed the gene expression patterns of the samples. Volcano plots are useful tools for visualizing differential expression between two different conditions. They are constructed using fold-change values and P-values, and thus allow visualization of the relationship between the magnitude of change and statistical significance. Of all the lncRNAs, 36 were up-regulated (fold change ≥2 and P < 0.05) and 17 were down-regulated (fold change ≥2 and P < 0.05) in the HDR group compared with the control group (Fig. 2A–C). Of all the mRNAs, 39 were up-regulated (fold change ≥2 and P < 0.05) and 35 were down-regulated (fold change ≥2 and P < 0.05) in the HDR group compared with the control group (Fig. 3A–C). Compared with the control group, these expression differences were not found in the LDR group, with the exception only one lncRNA which was down-regulated (fold change ≥2 and P < 0.05) (Table 2).

Figure 1.

Profile of microarray data (A,B) Differentially expressed lncRNAs (A) and differentially expressed mRNAs (B) were analyzed using hierarchical clustering. Hierarchical clustering analysis arranges samples into groups based on their expression level, which allows us to hypothesize the relationships between samples. ‘Red’ indicates high relative expression, and ‘green’ indicates low relative expression. (C,D) Box plot of lncRNA (C) and mRNA (D). They were used to observe and compare the distributions of gene expression values in the samples after normalization.

Figure 2.

LncRNA profile of the HDR group microarray data (A) Hierarchical Clustering shows a distinguishable lncRNA expression profile between the two groups and homogeneity within each group. (B) Volcano Plot of the differentially expressed lncRNAs. The red points in the plot represent the differentially expressed lncRNAs with statistical significance. (C) Scatter-Plot of lncRNA expression. The lncRNAs above the top green line and below the bottom green line indicated a more than 2.0-fold change in lncRNAs between the HDR group and control samples.

Figure 3.

mRNA profile of HDR group microarray data (A) Hierarchical Clustering shows a distinguishable mRNA expression profile between the two groups and homogeneity within each group. (B) Volcano Plot of the differentially expressed mRNAs. The red points in the plot represent the differentially expressed mRNAs with statistical significance. (C) Scatter-Plot of mRNA expression. The mRNAs above the top green line and below the bottom green line indicated a more than 2.0-fold change in mRNAs between the HDR group and control samples.

Table 2.

The pathways in the LDR group

Seq name	Gene symbol	P-value	Fold change	Regulation	Co-ordinates	Potential mechanism	Genomic relationship	Symbol	P-value	Fold change	Regulation	Co-ordinates
ENSMUST00000058706	Gm7816	0.0145	2.1366	Down	chr5	Enhancer	Upstream	Hs3st1	0.0226	1.3030	Up	chr5
AK047282		0.0011	1.4372	Up	chr5	Enhancer	Downstream	Hs3st1	0.0226	1.3030	Up	chr5
AK083162		0.0321	1.3419	Up	chr18	Neighboring	Upstream	Mbd2	0.0199	1.3159	Down	chr18
ENSMUST00000058881	Gm5341	0.0005	1.3116	Down	chr7	Enhancer	Downstream	Eed	0.0099	1.4367	Down	chr7

LncRNA target prediction

LncPath^TM arrays profile lncRNAs and their target genes simultaneously in a pathway-focused or disease-specific manner by quantifying lncRNAs (at the transcript level) and their potential target genes (at both the gene and transcript level) in parallel. In terms of the regulated target genes, up to 681 lncRNAs were co-expressed with target coding genes in the sample. The LDR group had four pathways that contained the predicted target genes (fold change >1.3, P < 0.05), while the HDR group had 17 pathways that were changed significantly (fold change ≥1.5 and P < 0.05) (Tables 2 and 3).

Table 3.

The pathways in the HDR group

Seq name	Gene symbol	P-value	Fold change	Regulation	Co-ordinates	Potential mechanism	Genomic relationship	Binding miRNA	Symbol	P-value	Fold change	Regulation	Co-ordinates
ENSMUST00000141797	Gm16194	0.0068	7.5347	Up	chr17	Enhancer	Upstream		Cdkn1a	0.0000	7.8613	Up	chr17
AK043271		0.0003	2.6043	Down	chr19	Enhancer	Upstream		Hells	0.0001	4.0311	Down	chr19
AK021390		0.0007	2.4458	Down	chr19	Enhancer	Upstream		Hells	0.0001	4.0311	Down	chr19
ENSMUST00000121087	Gm7308	0.0018	2.1328	Down	chr6	Enhancer	Downstream		Ccnd2	0.0000	3.0458	Up	chr6
uc008ref.2	F3	0.0002	18.9992	Up	chr3	ceRNA		mmu-miR-106a-5p, mmu-miR-106b-5p, mmu-miR-20a-5p, mmu-miR-98-3p, mmu-miR-17-5p, mmu-miR-467a-3p, mmu-miR-546, mmu-miR-669b-3p, mmu-miR-669f-3p	Asf1a	0.0003	2.9712	Down	chr10
AA089093		0.0002	5.2867	Up	chr6	Enhancer	Upstream		Ezh2	0.0008	2.6988	Down	chr6
ENSMUST00000118617	Gm12513	0.0010	3.2815	Down	chr4	Enhancer	Upstream		Klf4	0.0019	2.1494	Up	chr4
uc008ref.2	F3	0.0002	18.9992	Up	chr3	ceRNA		mmu-miR-29b-3p, mmu-let-7a-1-3p, mmu-let-7b-3p, mmu-let-7c-2-3p, mmu-let-7f-1-3p, mmu-miR-29a-3p, mmu-miR-29c-3p, mmu-miR-98-3p, mmu-miR-467a-3p, mmu-miR-770-5p, mmu-miR-669b-3p, mmu-miR-669f-3p	Hat1	0.0009	2.1308	Down	chr2
AK146532		0.0112	2.0950	Up	chr16	Enhancer	Downstream		Bcl6	0.0020	1.9466	Down	chr16
ENSMUST00000120026	Gm13588	0.0112	2.0131	Up	chr2	Enhancer	Downstream		Nek6	0.0051	1.9079	Up	chr2
ENSMUST00000181355	Gm26588	0.0070	3.1890	Up	chr6	Enhancer	Upstream		Mgll	0.0158	1.8993	Up	chr6
AK007040		0.0040	2.2962	Up	chr10	Enhancer	Upstream		Nuak1	0.0052	1.8268	Up	chr10
AK082634		0.0148	2.4210	Up	chr8	Enhancer	Upstream		Cdh11	0.0009	1.7708	Up	chr8
uc007xjz.1	2610037D02Rik	0.0026	2.0648	Up	chr15	Enhancer	Upstream		Arid2	0.0003	1.6935	Down	chr15
NR_045318	4930552P12Rik	0.0108	16.8170	Up	chr19	ceRNA		mmu-miR-27b-3p, mmu-let-7a-2-3p, mmu-miR-18a-3p, mmu-miR-22-3p, mmu-miR-671-5p, mmu-miR-666-3p, mmu-miR-873a-5p, mmu-miR-3096b-5p, mmu-miR-6370	Arid2	0.0028	1.5827	Down	chr15
NR_045162	A330048O09Rik	0.0083	2.3592	Up	chr13	Enhancer	Downstream		Id4	0.0063	1.5671	Up	chr13
AK047282		0.0036	5.1151	Up	chr5	Enhancer	Downstream		Hs3st1	0.0037	1.5441	Up	chr5

These co-expressed coding genes were identified as target genes regulated by differentially expressed lncRNAs. For instance, in the HDR group the down-regulated coding genes Asf1a (fold change: 2.97) were predicted to be the target genes regulated by the up-regulated lncRNA uc008ref.2 (fold change: 18.999), which has ‘ceRNA’-shared miRNA response elements (including mmu-miR-106a-5p, mmu-miR-106b-5p, and mmu-miR-20a-5p) with their corresponding target genes, and may prevent Asf1a from being degraded as a post-transcriptional gene regulator.

qRT-PCR validation

According to the lncRNA target prediction of lncRNA/mRNA expression profiling, and to confirm the reliability of our microarray data, we randomly selected two pathways (lncRNAs: ENSMUST00000058706, AK083162, and their target coding genes: Hs3st1, Mbd2) from the LDR group (fold change >1.3) and three pathways (lncRNAs: AA089093, AK043271, and uc008ref.2; and their target coding genes: Ezh2, Hells, and Asf1a) from the HDR group (fold change ≥2). These lncRNAs were selected based upon their regulation directions and significance. The expression levels in nine samples were analyzed by qRT-PCR. The qRT-PCR results were consistent with the microarray data and showed the same trends in up- or down-regulation for each gene (Fig. 4A–D). The primers designed by Sangon Biotech (Shanghai, China) are shown in Table 1.

Figure 4.

The expression profiles of the selected lncRNAs and their co-expressed mRNAs in the two groups (A,B) The heatmap of selected lncRNAs and their co-expressed mRNAs in the HDR group and LDR group. (C, D) The expression patterns of the genes as monitored by qRT-PCR.

Table 1.

Sequence of primers for qRT-PCR

Gene	Forward primer (5′→3′)	Reverse primer (5′→3′)
GAPDH	CGGAGTCAACGGATTTGGTCGTAT	AGCCTTCTCCATGGTGGTGAAGAC
AA089093	GCCTCCTTCAGCAAATCCTA	GTGTGGGCAGTGTTCAGGTA
AK043271	GGTAGAGGCTACAGTGAGCG	TGTTTTAGGCTGGGATGACG
uc008ref.2	CAACCCAAACCCACCAACTA	CCAGGTCACATCCTTCACG
ENSMUST00000058706	TGTTCTTTGGAGGAGGAGGA	ATTCTTTTGCCGAGCCAGT
AK083162	GAACCCTGCTGTTTGGCTTA	GCACAGCGGAGGAATTTCAA
Ezh2	CGGGACTAGGGAGTGTTCAG	GTTGTAAGGGCGACCAAGAG
Hells	CTTCCTAACTGGATGGCTGAA	TGTCCCTTGTCTTTTGTGGA
Asf1a	AAGGGAAAACCCACCAGTAAA	CTTTGACGCTGACGGTAATG
Hs3st1	TCTTTGACTGGGAGGAGCAT	TGAAATAGGCGGGTGTCTTC
Mbd2	CATTCACAGGGGGAGATACG	GTGCCTCCTCCAGTTTCTTG

GO analysis

To investigate the roles of lncRNAs following radiation, GO analysis, a functional analysis which associates differentially expressed mRNAs with GO categories, was performed. The results for the HDR group were obvious, thus the GO analysis was performed to determine the gene and gene product enrichment in biological processes, cellular components, and molecular functions in mouse thymocytes. Fisher's exact test was used to determine whether there was more overlap between the DE and the GO annotation list than that would be expected by chance (P < 0.05 is recommended). The results in Fig. 5 showed that the down-regulated transcripts were related to chromatin organization (GO: biological processes), chromatin binding (GO: cellular components), and the ESC/E(Z) complex (GO: molecular functions); while the highest enriched GOs targeted by up-regulated transcripts were related to cell death (GO: biological processes), nucleus (GO: cellular components), and kinase regulator activity (GO: molecular functions).

Figure 5.

The GO analysis (A) GO down-regulated and (B) GO up-regulated. The P-value denotes the significance of GO term enrichment in the DE genes. The lower the P-value, the more significant the GO term (P ≤ 0.05 was recommended as the cut-off value).

CNC network

In order to identify the possible modulating mechanism of lncRNAs in the HDR group, a CNC network was constructed based on the strength of the correlation between the differentially expressed mRNAs and lncRNAs (Fig. 6). LncRNAs and mRNAs with Pearson correlation coefficients that exceeded 0.995 were selected to construct the network. The three pairs of lncRNA pathways selected from the HDR group were shown in Fig. 6. Interestingly, these network pairs were all positive. Totally, 161 lncRNAs and 156 mRNAs formed the CNC network nodes and 228 network pairs were obtained. In this co-expression network, a single lncRNA may correlate with many mRNAs and vice versa. In addition, among these network pairs, there were also correlations between several top lncRNAs and top mRNAs. For instance, lncRNA ENST00000141797 (fold change: 7.5347226) and mRNA Cdkn1a (fold change: 7.861313) were predicted to co-express with a Pearson correlation coefficient of 0.997. These were one of the most significantly expressed lncRNAs and mRNAs, respectively. This selected pair pathway might be a center in this CNC network. All the targets mentioned here have been reported to be linked to chromatin binding and transcriptional regulation.

Figure 6.

Coding–non-coding gene co-expression network The blue circular nodes represent the mRNA, the orange box nodes represent the lncRNA. The solid lines represent the correlation. The node size indicates the number associated with these genes, where nodes with more gene co-expression have a more extensive relationship with the gene.

Discussion

Recently, lncRNAs have also been recognized as important regulators with biological functions. However, the expression patterns of lncRNAs and the classification of lncRNAs subgroups in the immune cells affected by radiation have not been fully clarified. In terms of immunization, Petri [15] provided an important resource for future studies on the functions of lncRNAs in the development of adaptive immune response. Subsequent studies have revealed that lncRNAs may impact the transcriptional programs of immune cells required to fight against pathogens and maintain normal health and homeostasis [16]. These studies suggest that non-coding RNAs play important roles in the immune cells.

In terms of radiation, photon radiation (X-rays and gamma-rays) is still a predominant type of IR used in radiotherapy with low linear energy transfer radiation compared with proton and alpha particles [17]. The mechanism of cell injury induced by low linear energy transfer radiation is mainly via free radicals, which induce DNA and cell membrane damage [18–20]. Various studies have confirmed that exposure to IR increases the risk of cancer [21–23]. The establishment of dose–response relationship for biological parameters is one of the topics in radiobiology research. Due to the importance of immune surveillance against cancer and infection, radiation effects on the immune system are one of the chief research fields in radiation biology and radiation protection. Our results revealed that some parameters were influenced in response to HDR and LDR radiation in mouse thymocytes. These parameters include: cell activities such as cell cycle control and cell survival (apoptosis), surface molecules, and signal molecules [24]. Therefore, it is essential to clarify the role of lncRNAs in the damage of immune cells when exposed to radiation.

In this study, we used 0.075 and 4.0 Gy of X-rays in the LDR and HDR, respectively, because these doses were used in most immunological studies to evaluate the response of lymphocytes to WBI [25]. Then RNAs were extracted from the thymus for the LncPath^TM microarray, and a small number of lncRNAs and their predicted target genes were selected for validation by quantitative RT-PCR. When reanalyzing the microarray data, the genes differentially expressed in the LDR group were not obvious, while those in the HDR group were obvious and therefore were analyzed further. By quantifying lncRNAs (at the transcript level) and their potential target genes (at both the gene and transcript level) in parallel, LncPath^TM Pathway Focus LncRNA Microarrays provided comprehensive insights into the expressional relationship between lncRNAs and their target genes, which was helpful in establishing rapid connections between the lncRNA regulatory mechanisms and their biological functions in immune cells, which were subject to WBI.

It is well known that HDR could induce DNA damage and severe DNA damage might lead to apoptosis or cell death [26–28]. In this study, GO analysis revealed that the genes biological processes involving in ‘cell death’, ‘programmed cell death’, and ‘proliferation’ were significantly enriched among the up-regulated genes in response to HDR. When reanalyzing microarray data in the HDR group, we found significant expression of coding genes, which are involved in ‘cell death’. For example, the lncRNA NR_028591 expression was elevated after HDR, which up-regulated the potential target genes BCL2L11 expression on the transcript level. The BCL2L11 is a member of the BCL-2 family and locates in the outer membrane of mitochondria, which acts as an important regulator that mediates excitotoxic apoptosis, apoptosis-inducing factor translocation and mitochondrial depolarization [29,30]. The GO analysis of down-regulated genes showed that these genes corresponded to ‘chromatin organization’ or ‘modification’, and ‘regulation of cellular macromolecule biosynthesis’ (Fig. 5B). For example, the lncRNA AA089093 expression was elevated, which down-regulated the potential target gene Ezh2 expression. Previous studies revealed that, in response to IL-1β and TGFβ2, TAGLN expression in endothelial cells is co-regulated at the chromatin level through histone methylation (H3K27me3), most likely through regulation of the abundance of the Polycomb methyltrasferase EZH2. EZH2 appears to integrate and mediate the IL-1β and TGFβ signaling at the epigenetic level [31]. Our previous research confirmed that HDR could induce the production of CD4⁺CD25⁺Nrp1⁺ cells, which inhibited immune responses in general through the release of TGF-β1 and IL-10. This process may be related to the PLC-PIP2 and cAMP-PKA signal pathways [32]. This provided an example of how cellular signaling could be resolved at the level of epigenetic regulation.

Therefore, the lncRNAs of these regulated target genes exert their biological functions at both the transcription and translation level. They regulate target genes by interacting with the transcriptional machinery or collaborating with chromatin modifiers directly. For instance, several lncRNAs, such as CCND1, gadd7, ANRIL, and PANDA, were induced by DNA damage. The expression of the lncRNA RoR, lncRNA p21, p53-induced eRNA, and loc285194 were modulated upon DNA damage in a p53-dependent manner [33–36]. In order to identify the possible modulating mechanism of lncRNAs in the HDR group, we constructed the CNC network based on the Pearson correlation analysis data of 161 lncRNAs and 156 mRNAs (protein-coding), which were associated with 228 network pairs of co-expressing lncRNAs and mRNAs. Most of these pairs showed a positive correlation. In Fig. 6, it was demonstrated that the selected genes such as those in the center in several circles and a single lncRNA could be correlated with many mRNAs or lncRNAs. For example, uc008ref.2 and ENSMUST00000141797 were closely related to coding RNAs or lncRNAs and had some in common, such as the coding gene Cdkn1a (also called p21Cip1), which was previously proposed as a biomarker gene after radiation exposure. CDKN1A is a conserved protein belonging to the CDK inhibitor Cip/Kip family, and has been implicated in the regulation of many cellular processes, such as proliferation, differentiation, apoptosis, metastasis, cell survival, and stem cell renewal. Expression of CDKN1A could be regulated at the transcriptional level by oncogenes and tumor suppressor proteins that bind various transcription factors to specific elements in response to a variety of intracellular and extracellular signals [37,38]. Therefore, we speculate that uc008ref.2 and ENSMUST00000141797 may respond to HDR by interacting with cellular processes through CDKN1A. Other mRNAs were also co-expressed with their corresponding lncRNAs. All the targets mentioned here have also been reported to be linked to ‘cellular mechanisms’, ‘transcriptional regulation’, and ‘chromatin binding’, consistent with the GO analysis results. Thus, lncRNAs might interact with these targets. These findings suggest that lncRNAs and their target gene network may play important roles in the regulation of cellular response and biological response to IR damage. However, further studies are required to determine how these lncRNAs are involved in the transcriptional and post-transcriptional regulation of their related genes.

In summary, this study demonstrates the profile of differentially expressed lncRNAs and their co-expressed coding genes involved in the damage responses of thymocytes to LDR and HDR. It is suggested that lncRNAs may also act as co-regulators of transcription factors to regulate gene transcription in immune cells in response to IR, or to affect the expression of adjacent coding genes. Considering that only a small percentage of the total lncRNA population has been studied up to date, we are just at the beginning of this research. Further work is needed to clarify the significance of lncRNAs in response to IR damage, which may have a profound impact on our understanding of the biological and molecular mechanisms involved in IR-induced DNA damage response, and may result in novel therapeutic approaches.

Funding

This work was supported by the grants from the National Natural Science Foundation of China (Nos. 30870584, 81371890, and 81573085) and the Doctoral Program Foundation of Institutions of Higher Education of China (No. 20120061110063).

Author notes

†

These authors contributed equally to this work.