Whole Transcriptome Analysis Revealed a Stress Response to Deep Underground Environment Conditions in Chinese Hamster V79 Lung Fibroblast Cells

Background: Prior studies have shown that the proliferation of V79 lung fibroblast cells could be inhibited by low background radiation (LBR) in deep underground laboratory (DUGL). In the current study, we revealed further molecular changes by performing whole transcriptome analysis on the expression profiles of long non-coding RNA (lncRNA), messenger RNA (mRNA), circular RNA (circRNA) and microRNA (miRNA) in V79 cells cultured for two days in a DUGL. Methods: Whole transcriptome analysis including lncRNA, mRNAs, circ RNA and miRNA was performed in V79 cells cultured for two days in DUGL and above ground laboratory (AGL), respectively. The differentially expressed (DE) lncRNA, mRNA, circRNA, and miRNA in V79 cells were identified by the comparison between DUGL and AGL groups. Quantitative real-time polymerase chain reaction(qRT-PCR)was conducted to verify the selected RNA sequencings. Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was analyzed for the DE mRNAs which enabled to predict target genes of lncRNA and host genes of circRNA. Results: With |log2(Fold-change)| ≥ 1.0 and p < 0.05, a total of 1257 mRNAs (353 mRNAs up-regulated, 904 mRNAs down-regulated), 866 lncRNAs (145 lncRNAs up-regulated, 721 lncRNAs down-regulated), and 474 circRNAs (247 circRNAs up-regulated, 227 circRNAs down-regulated) were significantly altered between the two groups. There was no significant difference in miRNA between the two groups. The altered RNA profiles were mainly discovered in lncRNAs, mRNAs and circRNAs. DE RNAs were involved in many pathways including ECM-RI, PI3K-Akt signaling, RNA transport and the cell cycle under the LBR stress of the deep underground environment. Conclusion: Taken together, these results suggest that the LBR in the DUGL could induce transcriptional repression, thus reducing metabolic process and reprogramming the overall gene expression profile in V79 cells.

Background: Prior studies have shown that the proliferation of V79 lung fibroblast cells could be inhibited by low background radiation (LBR) in deep underground laboratory (DUGL). In the current study, we revealed further molecular changes by performing whole transcriptome analysis on the expression profiles of long non-coding RNA (lncRNA), messenger RNA (mRNA), circular RNA (circRNA) and microRNA (miRNA) in V79 cells cultured for two days in a DUGL.
Methods: Whole transcriptome analysis including lncRNA, mRNAs, circ RNA and miRNA was performed in V79 cells cultured for two days in DUGL and above ground laboratory (AGL), respectively. The differentially expressed (DE) lncRNA, mRNA, circRNA, and miRNA in V79 cells were identified by the comparison between DUGL and AGL groups. Quantitative real-time polymerase chain reaction(qRT-PCR)was conducted to verify the selected RNA sequencings. Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was analyzed for the DE mRNAs which enabled to predict target genes of lncRNA and host genes of circRNA.

INTRODUCTION
An increasing number of countries have begun to develop deep Earth in order to cope with the space and resources of surface Earth being gradually consumed in the future (Xie et al., 2018). As China, exploring deep underground space and resources has become a national priority since 2016 (Xie et al., 2017), which shown that an increasing number of people could live and/or work in the underground space in the near future (Ranjith et al., 2017;Liu et al., 2018). Currently in South Africa, gold miners are able to work 4000 m underground (Ranjith et al., 2017). However, little is still known about the biological effects of the deep underground environment . The limited knowledge may result from the shortage of deep underground laboratories (DUGLs),such as the Gran Sasso National Laboratory (LNGS) in Italy, and the Waste Isolation Pilot Plant (WIPP) in the United States. Researchers who have historically focused on the growth of cultures in DUGLs have observed some interesting changes(e.g., cell growth delay, enzyme activity change and sensitivity to genetic damage factors) in living organisms exposed to low levels of radiation (Belli, 2017;Castillo and Smith, 2017;Liu et al., 2018). However, with limited publications it is challenging to know the reliability of the results (Castillo and Smith, 2017).
To understand the biological effects of the deep underground environment on humans and service for the development of deep underground space and resources, a DUGL in a cave with a rocky cover of 1470 m was set up in our previous research in Erdaogou Mine, with Jiapigou Minerals Limited Corporation of China National Gold Group Corporation (CJEM) . Similar to the 4000 m water equivalent (WE) of LNGS, the flux of the cosmic rays in the DUGL of CJEM could be considered negligible compared to the cosmic rays at the above ground surface (Liu et al., 2020). The radon concentration in the DUGL of CJEM was 3.7-5.5 pCi/L, which is slightly higher than LNGS (Liu et al., 2020). However, the total gamma (γ) radiation dose rate of terrestrial radiation was 0.03-0.05 µSv/h,which was consistent with the levels at LNGS (Liu et al., 2020). Meanwhile, the relative humidity was approximately same as LNGS at 99%. Except for oxygen (O 2 ), the concentration of carbon dioxide (CO 2 ) and air pressure of the DUGL were slightly higher than an above ground laboratory (AGL) (Liu et al., 2020).
Prior series studies have investigated biological changes of Chinese hamster V79 cells in LNGS (Antonelli et al., 2000;Satta et al., 2002;Carbone et al., 2009;Fratini et al., 2015). Published research also demonstrated that V79 cells were able to be cultured in the DUGL of CJEM and the proliferation of these V79 cells could be inhibited by low background radiation (LBR) in this deep underground environment . Furthermore, our prior study found that V79 cells cultured in the DUGL of CJEM presented an altered protein profile related to the ribosome, RNA transport, translation, energy metabolism, and DNA repair (Liu et al., 2020). Recent evidence revealed that non-coding RNAs participate in modulating numerous biological functions by regulating gene expression at transcriptional and post-transcriptional levels (Beermann et al., 2016). In the current study, we revealed further molecular changes by performing whole transcriptome analysis on the expression profiles of long non-coding RNA (lncRNA), messenger RNA (mRNA), circular RNA (circRNA) and microRNA (miRNA) in V79 cells cultured for two days in a DUGL. These data will be helpful to further understand the biological effects of the deep underground environment.

Cell Culture
The detailed cell culture methods are described in our previous research (Liu et al., 2020). Briefly, frozen Chinese hamster V79 lung fibroblast cells were purchased from Shanghai Enzymelinked Biotechnology (China) and cultured in Dulbecco's modified eagle medium (DMEM; Gibco, United States) with 10% fetal calf serum (Gemini, United States), and 50 U dm −3 penicillin and streptomycin (Gibco, United States). At more than 80% confluency, cells were passaged and divided into four T 25 flasks, and two of them were assigned to either the DUGL or AGL of CJEM, respectively. All cells were maintained in incubators at 37 • C with 5% CO 2 . One flask of each location was used for growth curve and morphology analysis in previous researched (Liu et al., 2020), the another flask was cultured for passage. After passaging and two days culture, three flasks with enough cells for further analysis from each location were collected for the following experiments.

RNA Preparation and Quality Control
After the cells were cultured for two days in either the AGL or DUGL, the total RNA per sample was extracted using Trizol reagent (Invitrogen, NY, United States) and was then used for RNA sequencing. The purity and integrity of RNA were examined by 1% agarose gel electrophoresis (Sigma-Aldrich,United States). Subsequently, further RNA integrity was verified using the Agilent 2100 Bioanalyzer (Agilent Technologies, CA, United States), The RNA integrity number (RIN) of all samples were more than 7.0, which were considered to meet the sequencing requirement. The RNA quantity was measured using the NanoDrop-2000 (NanoDrop Technologies, DE, United States). The ribosomal RNA (rRNA) was removed with the Ribo-Zero GoldKits (Epicentre, WI, United States).
lncRNA and mRNA Sequencing and Data Processing RNA (3 µg) from each sample was used for cDNA library construction (NEB Next Ultra Directional RNA LibraryPrep Kit for Illumina, Ispawich, United States). After the removal of rRNA (Ribo-Zero GoldKits), the rRNA-depleted RNAs were fragmented and used as templates to construct the cDNA library. First strand cDNA was synthesized using random hexamer primers. Next, second strand cDNA synthesis was performed using DNA polymerase I and RNase H. Libraries were amplified by polymerase chain reaction (PCR).
NovaSeq 6000 Illumina sequencing system (Illumina, San Diego, CA, United States) was used for RNA sequencing. The sequencing data was analyzed using CASAVA software for base calling. Raw data were transformed into FASTQ stored documents. To obtain clear reads for further analysis, reads containing adapters and ploy-N as well as low-quality reads were removed from the raw data. HiSAT2 software 1 was used to align sequencing data according to the reference hamster genome [GCA_003668045.1 (Rupp et al., 2018), GenBank assembly accession].
The Coding-Non-Coding Index (CNCI), Coding Potential Calculator (CPC), Pfamscan, and Coding Potential Assessment Tool (CPAT) were used to analyze the coding potential of transcripts. The assembled transcripts without coding potential of their overlap became the candidate set of lncRNAs.
The read count for each gene in each sample was counted by HTSeq, and Fragments Per Kilobase of transcript, per Million mapped reads (FPKM) were then calculated to represent the expression level of mRNA and lncRNA in each sample. The DE mRNAs and lncRNAs between the two groups were identified using DEseq. The DE cut-off criteria included q < 0.05 (adjusted p value) and |log 2 (Fold-change)| ≥ 1.0.

circRNA Sequencing and Data Processing
The other cDNA library construction method was similar to that for the lncRNAs, except RNase R was used to remove linear RNAs. Clean reads were obtained by removing the following fragments: (1) low quality data, (2) reads containing an N ratio greater than 5%, and (3) reads containing jointed-sequence and rRNAs. BWA-MEM (Li, 2013) was used for mapping clean reads to the reference genome. circRNA Identifier (CIRI) , a highly efficient and fast circRNA recognition tool, was used for circRNA recognition. The BWA-MEM algorithm was used for sequence splitting, and then the resulting SAM file was scanned to detect paired chiastic clipping (PCC) and paired-end mapping (PEM) sites, as well as GT-AG splicing signals. Next, the sequence of junction sites was re-aligned using a dynamic programming algorithm to ensure the reliability of the identified circRNA. Spliced Reads per Billion Mapping (SRPBM) were used to determine the expression level of circRNA. The DE circRNAs between the two groups were identified using DEseq (Love et al., 2014). The cut-off criteria included p < 0.05 and |log 2 (Foldchange)| ≥ 1.0. circRNA-miRNA interactions were predicted using the miRanda (3.3a) prediction algorithm.
miRNA Sequencing and Data Processing 18-30 nucleotide (nt) or 15-35 nt RNA fragments were obtained through agarose gel separation technology (Tao et al., 2019). Next, those RNAs were reverse transcribed to synthesize cDNA. Clean reads were obtained by removing the following fragments: (1) low quality data, (2) reads containing an N ratio greater than 10%, (3) reads without a 3' linker sequence, (4) reads with polyA/T, and (5) reads without sequences. Bowtie1 (Langmead, 2010) was used to map the clean reads to the reference genome. Reads were mapped to mature miRNA and hairpin RNA that were recorded in miRBase (release 21) (Griffiths-Jones et al., 2008) to identify known miRNAs. After excluding the reads mapped to known miRNA/ncRNA/repeat regions/mRNA regions, the remaining reads were used to predict novel miRNAs based on their hairpin structure and stability. The miRDeep2 (Friedländer et al., 2011) software was applied for the identification and prediction of miRNAs. Transcripts per Million (TPM) was adopted to determine the expression levels of miRNAs. With a cut-off of q < 0.05 and | log 2 (Fold-change)| ≥ 1.0, the DESeq2 package in R software was employed to identify the DE miRNAs.

Verification by Real-Time Quantitative Polymerase Chain Reaction
To verify the RNA-Seq results, 10 DE mRNAs were randomly selected for Real-time quantitative polymerase chain reaction (qRT-PCR) analysis (Table 1). Primer-BLAST of NCBI 2 was used for primer design. Total RNA was reverse transcribed into cDNA using a PrimeScript RT Reagent Kit with gDNA Eraser according to the manufacturer's instructions (RR047A, Takara, Japan). A 7500 Real-time PCR system (Applied Biosystems, CA, United States) was used to perform qRT-PCR. Actin was selected as the reference gene. The forward primer sequence was GATCTGGCACCACACCTTCT, and the reverse was GGGGTGTTGAAGGTCTCAAA. Three repeats were performed for each group. The relative gene expression levels were calculated by the 2 − Ct method (Livak and Schmittgen, 2001). qRT-PCR data were analyzed using a t-test with SPSS 13.0 software. A p value ≤0.05 was considered statistically significant.

Biological Information Analysis
Hierarchical clustering was conducted using R software (v3.5.1, 2018). Since the functional annotation of most lncRNAs has not been obtained, the functions of lncRNAs were predicted according to the annotations of the co-expressed mRNA function (Han et al., 2018). In this study, the function of DE  lncRNAs were predicted based on position relationship (within 50 kb of lncRNA) and the Pearson correlation coefficient (the value of correlation ≥0.9, and p < 0.01) between lncRNA and mRNA. The lncRNAs/circRNAs/miRNAs and mRNAs that shared expression levels with significant correlations were used to conduct co-expression analyses. To further reveal the biological functions of the DE RNAs, GO and KEGG pathway enrichment (performed using GeneCodis3 bioinformatics resources) were applied to the DE mRNAs, predicted targets for DE lncRNAs and miRNAs, and host genes of circRNAs. Both GO terms and KEGG pathways with q values < 0.05 were considered significantly enriched.

Overview of Transcriptomic Analyses
Sequencing was performed on the cDNA and sRNA libraries of cell samples from the groups of DUGL and AGL cells grown for two days. Furthermore, total read counts and the ratio of mapped reads of the sequencing data are shown in Table 2.
Using a cut-off of | log 2 (Fold-change)| ≥ 1.0 and q < 0.05, a total of 1257 mRNAs and 866 lncRNAs were identified as being differentially expressed (DE) between the two groups. Using a | log 2 (Fold-change)| ≥ 1.0 and p < 0.05, 474 DE circRNAs were identified. However, only nine novel miRNAs were found to be down-regulated in DUGL cells, but these changes were not statistically significant. Among the DE RNAs, there were 353 up-regulated mRNAs, 904 down-regulated mRNAs, 145 up-regulated lncRNAs, 721 down-regulated lncRNAs, 247 upregulated circRNAs, and 227 down-regulated circRNAs in the DUGL cells (Figure 1). The top 10 dysregulated lncRNAs, circRNAs and mRNAs are shown in Tables 3-5, respectively.
Hierarchical clustering of the lncRNA, mRNA and circRNA expression suggested obvious discrimination in V79 cells between the DUGL and AGL growth conditions (Figure 1).

The Functional Analysis of DE mRNAs
According to GO analyses, in the biological process (BP) category, the down-regulated mRNAs in the DUGL group were enriched in 28 terms ( Table 6). The top three enriched terms were response to external stimulus(GO:0009605), defense response(GO:0006952) and response to stimulus(GO:0050896) (Figure 2A and Table 2). As to the molecular function (MF) category, the down-regulated mRNAs were mainly enriched in terms of oxidoreductase activity (GO:0016614), transmembrane signaling receptor activity (F) Hierarchical cluster analysis of DE circRNAs. DUGL, deep underground laboratory; AGL, above ground laboratory. In the volcano plots, red and green dots correspond to RNAs that are significantly up-regulated or down-regulated between the two groups [|log 2 fold change| ≥ 1.0 and q < 0.05 or (p < 0.05 for circRNAs)]. The X-axis shows the log 2 fold changes of RNAs expression, and the Y -axis shows the adjusted p-value (−log 10 ) for each gene. In the hierarchical cluster analyses, blue represents down-regulated, black represents no change, and yellow represents up-regulated.
Frontiers in Genetics | www.frontiersin.org The KEGG analysis showed that the down-regulated mRNAs were enriched in six pathways including extracellular matrixreceptor interaction (ECM-RI), arachidonic acid metabolism, linoleic acid metabolism, NOD-like receptor signaling pathway, glutamatergic synapse and PI3 kinase-Akt signaling pathway ( Figure 3C and Supplementary Table 3). There was no significant enrichment for up-regulated mRNAs in any pathway.

The Functional Analysis of DE lncRNAs
To investigate the DE lncRNAs under the LBR stress of the deep underground environment, a functional analysis was performed for the predicted target mRNA of the DE lncRNAs  The BP terms were mainly related to several steps of the metabolic process (Figure 4A), the CC terms were related to organelle ( Figure 4B) and the MF terms were related to binding ( Figure 4C). In contrast, the predicted target mRNAs of the up-regulated lncRNAs were mainly enriched in 214 BP terms, 101 CC terms and 53 MF terms ( Table 6 and  Supplementary Table 5). In the BP category, the main enriched terms were related to metabolic processes (GO:0008152), such as primary (GO:0044238), organic substance(GO:0071704) and macromolecule metabolic (GO:0043170) ( Figure 5A). Part of those terms were enriched in the cellular response to stress (GO:0033554), and regulation of cellular response to stress (GO:0080135) (Supplementary Table 5). In the CC category, the target mRNAs of DE lncRNAs mainly related to organelle part (GO:0044422), organelle (GO:0043226), and nuclear part (GO:0044428) ( Figure 5B). As for the MF category, most of the identified terms related to binding (GO:0005488), catalytic (GO:0003824), and protein binding (GO:0005515) functions ( Figure 5C).
In the KEGG analysis, the target mRNAs of down-regulated lncRNAs were significantly enriched in four pathways (spliceosome, RNA degradation, proteasome and protein processing in the endoplasmic reticulum) (Figure 4D). Interestingly, the target mRNAs of up-regulated lncRNAs were enriched in these same three pathways and one other pathway (spliceosome, RNA degradation, proteasome and cell cycle) (Figure 5D).  were detected in this study. Function analyses were then performed to identify the host genes of these DE circRNAs. GO analyses showed that host genes of the down-regulated circRNAs were enriched in four BP terms [(cellular macromolecule metabolic process (GO:0044260), regulation of macromolecule metabolic process (GO:0060255), macromolecule metabolic process (GO:0043170) and regulation of metabolic process (GO:0019222)];eight CC terms including intracellular part (GO:0044424) and organelle (GO:0043229) and nine MF terms including Rab GTPase binding(GO:0017137), protein kinase activity(GO:0004672) and phosphotransferase activity(GO:0016773) ( Table 6 and Figures 6A-C, and Supplementary Table 6). The host genes of the up-regulated circRNA were enriched in 31 CC terms including intracellular part(GO:0044424) and organelle(GO:0043229);and three MF terms [GTPase activator(GO:0005096) and regulator activity(GO:0030695), transferase activity(GO:0016740)] ( Table 6 and Figures 7A,B). Several terms surrounding metabolic processes were enriched, which was similar to the target mRNAs of the DE lncRNAs. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses showed that the host genes of downregulated circRNA were enriched in protein processing in the endoplasmic reticulum.

Construction of the circRNA-miRNA Co-expression Network
A circRNA-miRNA co-expression network was constructed based on the RNA-Seq results, and when comparing DUGL to AGL cells, 286 miRNA-circRNA interaction pairs were obtained. For interest, the miRNA-circRNA interaction network was shown in Figure 8 and Supplementary Table 8.

Verification of DE RNA by qRT-PCR
To verify the RNA-Seq results, 10 DE mRNAs were selected for qRT-PCR analysis. Among them, nine mRNAs had comparable expression patterns between the RNA-Seq and qRT-PCR results (Figure 9).

DISCUSSION
Several researchers have found that the deep underground environment, where cosmic radiation is shielded, can reduce the growth rates of paramecia, bacteria and some mammalian cells (Planel et al., 1987;Smith et al., 2011;Castillo et al., 2015). Indeed, cell translational repression and gene expression profiles induced by stress can inhibit proliferation . However, as an environmental stress, little is known about the genetic profile changes that occur under the stress of LBR in a deep underground environment. In this study, whole transcriptomic analyses were conducted in V79 cells grown for two days in a DUGL with LBR and compared to cells grown in an AGL. The results showed a distinct genetic profile change, which covered circRNAs, mRNAs, and lncRNAs (most of them down-regulated). Although there was no significant difference in miRNA between the two groups, nine DE novel miRNAs were identified. Ten DE mRNAs found by RNA-Seq were selected for qRT-PCR validation, and a similar expression level in nine of the ten mRNAs confirmed the accuracy of the RNA-seq findings to some extent. Taken together, the changes in the gene profile suggested that LBR stress could cause a delay in growth through the inhibition gene transcription. Therefore, the genetic down-regulation induced by LBR stress might be the main molecular basis of the inhibition of proliferation in V79 cells cultured in DUGL.
The GO term annotation is helpful to reveal physiological and functional changes related to genes and protein expression in cells (Camon et al., 2004). To further explore the effect of the genetic changes in V79 cells cultured in DUGL, GO analyses were performed for the DE mRNAs. Under LBR environmental stress, the top three BP terms of the down-regulated mRNAs were related to responses to stimulus and defense, and the down-regulated mRNA mainly enriched in CC terms of plasma membrane. The cell senses and adapts to changes in the extracellular environment by plasma membrane,which allows the external singal inputing via intracellular signaling networks (Rausch and Hansen, 2020). The down-regulated mRNAs of plasma membrane related to stress function might help to explain the hypothesis that normal environmental radiation contributes to maintaining the defense systems of living organisms (Fratini et al., 2015). Contrast to γ irradiation could induce the expression of the interleukin-1(IL-1) gene (Bigildeev et al., 2013), our study showed that the LBR could decrease the IL-1 gene production, which presented with the down-regulate mRNAs enriched in ten GO terms involving IL-1 production. On the other hand, the up-regulated DE mRNAs were significantly enriched in BP categories involved in negative regulation terms, such as gene expression, metabolic process, and biosynthetic process. These up-regulated mRNAs related to negative metabolic and biosynthetic functions also could be the main causative factors of proliferative inhibition.   The KEGG pathway analysis is useful for the systematic understanding of large-scale gene functions. In this study, KEGG pathway analysis of DE mRNAs showed significant enrichment in ECM-RI, arachidonic acid metabolism, linoleic acid metabolism, NOD-like receptor signaling pathway, glutamatergic synapse and PI3K-Akt signaling pathway. As a biological regulation network, both the ECM-RI and PI3K-Akt signaling pathways were comprehensive net and played an essential role in cell proliferation and survival (Landgren and Curtis, 2011;Zhang G. et al., 2018). Twelve down-regulated were mRNAs shared in the two pathways. Most of the down-regulated genes had the function of promoting proliferation, such as collagen alpha 1 (Col1a1) (Jin et al., 2017), integrin alpha 7 (Itga7) (Ge et al., 2020), laminin beta 3 (Lamb3)  and secreted phosphoprotein 1 (Spp1) (Zeng et al., 2018). Furthermore, other down-regulated genes were detected with similar functions [e.g., Rab40b (Jacob et al., 2016), S100 calcium-binding protein A4 (S100A4) (Forman et al., 2010) and collagen prolyl-4hydroxylase α subunit 2 (P4HA2)]. These key down-regulated genes might play crucial roles in the inhibition of growth found in DUGL cultures. Moreover, arachidonic and linoleic acid metabolism, as well as NOD-like receptor signaling and the glutamatergic synapse might also be involved in the stress response of LBR. However, these pathways' role in LBR stress need further research.
Besides the DE mRNAs enriched in GO terms and KEGG pathways, several top DE mRNAs function in the regulation of cell proliferation. Matrix metallopeptidase 9 (MMP9) has been shown to be involved in promoting proliferation (Lan et al., 2020). Human concentrative nucleoside transporter-1 (hCNT1, SLC28A1) (Wang and Buolamwini, 2019), Ccl20 (Qiu et al., 2011), Sema4d (Bhutia et al., 2011), and Arhgap36 (Rack et al., 2014) have been shown to influence cellular growth and proliferation, and were the top down-regulated mRNAs in V79 cells cultured in a DUGL. Nuclear receptor 4A3 (NR4A3) has been shown to have the ability to suppress cell growth (Zhang B. et al., 2018;Parsa et al., 2019). Importantly, NR4A3 was significantly up-regulated in cells cultured in a DUGL. Therefore, it can be inferred that the down-regulation of genes that promote growth and the up-regulation of genes that function to suppress growth also contribute to the inhibition of proliferation in cells cultured in a DUGL.
Non-coding RNAs, by definition, do not code for protein. However, they have crucial roles in various cellular activities (Kukleva et al., 2015). Therefore, the lncRNA, circRNA, and miRNA in V79 cells with altered expression profiles between those cultured in a DUGL and AGL were comprehensively investigated. lncRNA is a class of non-coding transcripts longer than 200 nucleotides (Cao et al., 2018). In the present study, 866 lncRNAs were found to be differentially expressed between the two groups. The number of down-regulated lncRNAs was much higher than those that were up-regulated, suggesting that the expression of lncRNAs transcripts is repressed by LBR stress. Owing to a few studies in this field, more than half of the DE lncRNAs were newly identified in this study.
To further reveal the functions of the identified DE lncRNAs, GO and KEGG analyses were performed. In the BP category, both up-and down-regulated target mRNAs of lncRNAs were mainly enriched in many terms related to metabolic processes. This finding strongly suggested that lncRNAs, similar to mRNAs, might also play a crucial role by altering metabolic processes during LBR stress in a deep underground environment. In the KEGG analyses, the target mRNAs of dysregulated lncRNAs shared three pathways including spliceosome, RNA degradation and proteasome. Spliceosomes are known to precisely and efficiently perform mRNA processing, which is a critical step in organ development. Consistent with proteomic result of our previous research (Liu et al., 2020), the target mRNAs of lncRNAs enriched in these pathways indicated that the spliceosome played an important role in the LBR stress response of V79 cells. Additionally, the involvement of the RNA degradation and proteasome pathways might suggest that these pathways also function in the stress response in a LBR environment.
Regarding another class of non-coding RNA, accumulating evidence has highlighted that circRNAs can affect mRNA splicing and transcription (Zhang et al., 2019). Therefore, we analyzed the DE circRNAs between DUGL and AGL groups of V79 cells and identified 474 DE circRNAs. GO analyses showed that the host genes of down-regulated DE circRNAs were enriched in metabolic processes, which was similar to the results found in both DE mRNAs and target genes of DE lncRNAs. This finding indicated that circRNAs are also involved in the LBR stress of V79 cells by interacting with lncRNAs and mRNAs. The endoplasmic reticulum (ER) is a vital organelle that can perceive environmental changes . ER stress could induce changes in key mediators for cell survival (Chien et al., 2017). KEGG analyses revealed that the host genes of down-regulated DE circRNAs were enriched in the pathway of protein processing in the ER. This result was consistent with our previous studies, which have revealed that several proteins of the ER were down-regulated in cells cultured in a DUGL. However, further confirmation is required to verify the functional role of the ER in LBR stress. miRNAs, as a class of small non-coding RNAs (approximately 22 nucleotides), are essential elements to regulate gene expressions through partial base-pairing with target mRNAs (Lin et al., 2018). circRNAs may function similarly to regulate the activity of other miRNAs (Wilusz and Sharp, 2013). Due to little research in this area, and since the minority of the miRNA functions were annotated, we failed to detect significance difference in the DE miRNAs between the two groups. However, we identified several circRNAs that contained one or more miRNA binding sites and obtained 286 miRNA-circRNA interaction pairs between the DUGL and AGL groups. These interactions and their functions involved in LBR stress are worthy of being investigated further.
Nucleic acid binding plays an important role in translation regulation (Liu et al., 2020). In our present study, the upregulated mRNA and target mRNAs of dysregulated lncRNAs enriched in many GO terms which were related to nucleic acid binding. The change of gene expression was consistent with the proteomic result of V 79 cells conducted in our previous research (Liu et al., 2020). Those molecular change from RNAs to proteins further indicated that these nucleic acid binding was affected when V 79 cells under the stress of reduced background radiation.
Environmental stress can trigger an increase in reactive oxygen species (Fan et al., 2015). Castillo and colleagues (Castillo et al., 2015) have shown that oxidative stress and the SOS response (katB and recA) as well as metal efflux activity (SOA0154) were elevated in cells grow under LBR conditions. Although we did not detect the differential expression of these specific genes between the two groups, down-regulated DE mRNAs in DUGL cells were found to be enriched in two terms involving oxidoreductase activity. Similar to DE RNAs, the predicted target genes of DE lncRNAs were observed to be enriched in cellular responses to oxygen-containing compounds and mitochondrial terms. These findings were also consistent with our previous research, which showed that cells grown in a DUGL presented a change in energy metabolism, morphologic changes of mitochondrion and oxidative phosphorylation (Liu et al., 2020). Taken together, the results of this study suggest that the oxidative response could be involved in the LBR stress response in the deep underground environment.
The limitations of our study included the short growth time of V79 cultured cells (two days) in the DUGL for the analysis of lncRNAs, circRNAs, and miRNA. Also, we did not verify the sequencing results for DE lncRNAs and circRNAs by RT-qPCR. Additionally, we failed to construct a competing endogenous network and verify an interactive relationship. Therefore, additional in-depth research is required to reveal the biological specifics of the LBR stress response in a deep underground environment.

CONCLUSION
In conclusion, our study investigated the transcription patterns of lncRNAs, mRNAs, circRNAs, and miRNAs of V79 cells cultured in a DUGL and an AGL by whole-transcriptome sequencing and integrated analysis. We confirmed that the LBR of a deep underground environment could induce V79 cell transcription repression, metabolic process delaying and overall gene expression profile reprogramming. The altered RNA profiles were mainly discovered in lncRNAs, mRNAs and circRNAs. DE RNAs were involved in many pathways including ECM-RI, PI3K-Akt signaling, RNA transport and the cell cycle under the LBR stress of the deep underground environment. These profile changes might be the molecular basis of the inhibition of cell proliferation. This study provided a systematic perspective on the potential effects of the deep underground environment on V79 cells.

DATA AVAILABILITY STATEMENT
The datasets generated and/or analyzed during the current study can be found with this article and Supplementary Information online. The raw transcriptomic data have been deposited to the NCBI Short Read Archive (SRA) under the BioProject accession number PRJNA623056: https://www.ncbi.nlm.nih.gov/ bioproject/PRJNA623056.