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 |log₂(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.

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 (Liu et al., 2018). 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) (Liu et al., 2018). 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₂), the concentration of carbon dioxide (CO₂) 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 (Liu et al., 2018). 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.

Materials and Methods

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 Enzyme-linked 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₂. 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¹ 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₂(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) (Gao et al., 2015), 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₂(Fold-change)| ≥ 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₂(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² 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.

TABLE 1

Table 1. Primer sequences used in quantitative reverse-transcription polymerase chain reaction (RT-qPCR)-based verification of RNA sequencing results.

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.

Results

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.

TABLE 2

Table 2. Summary of sequence after Illumina sequencing.

Using a cut-off of | log₂(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₂(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 up-regulated 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.

FIGURE 1

Figure 1. Summary of the differentially expressed (DE) RNAs between V79 cells cultured in either DUGL or AGL conditions. (A) Volcano plot of DE mRNAs. (B) Hierarchical cluster analysis of DE mRNAs. (C) Volcano plot of DE lncRNAs. (D) Hierarchical cluster analysis of DE lncRNAs. (E) Volcano plot of DE circRNAs. (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₂ fold change— ≥ 1.0 and q < 0.05 or (p < 0.05 for circRNAs)]. The X-axis shows the log₂ fold changes of RNAs expression, and the Y-axis shows the adjusted p-value (−log₁₀) for each gene. In the hierarchical cluster analyses, blue represents down-regulated, black represents no change, and yellow represents up-regulated.

TABLE 3

Table 3. The top ten differently expressed long non-coding RNAs (lncRNA) when comparing cells grown in either deep underground lab (DUGL) or above ground lab (AGL) conditions.

TABLE 4

Table 4. The top ten differently expressed circular RNAs (circRNA) when comparing cells grown in either deep underground lab (DUGL) or above ground lab (AGL) conditions.

TABLE 5

Table 5. The top ten differently expression messenger RNAs (mRNA) when comparing cells grown in either deep underground lab (DUGL) or above ground lab (AGL) conditions.

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 Supplementary Table 1). The up-regulated mRNAs in the DUGL group were enriched in 102 BP terms (Table 6, Figure 3A and Supplementary Table 2). Among those BP terms, 19 terms were negative regulation terms, which covered gene and metabolic processes. Moreover, seven negative regulation terms ranked in the top 10 enriched terms (Figure 3A and Supplementary Table 2). In the cellular component (CC) category, the down-regulated mRNAs were mainly enriched in terms of membrane part [e.g., plasma membrane (GO:0005886), extracellular region part (GO:0044421) and extracellular space(GO:0005615) (Figure 2B and Supplementary Table 1)]; whereas the up-regulated mRNAs were only enriched in terms of the nucleolus (GO:0005730) (Supplementary 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 (GO:0004888) and metalloendopeptidase activity(GO:0004222) (Figure 2C and Supplementary Table 1); whereas the up-regulated mRNAs were mainly enriched in terms of binding [e.g., transcription factor(GO:0008134),nucleic acid(GO:0003676) and regulatory region nucleic acid(GO:0001067) (Figure 3B and Supplementary Table 2)].

TABLE 6

Table 6. The overview of functional analysis of differentially expressed RNAs.

FIGURE 2

Figure 2. Functional annotation of the differentially expressed (DE) down-regulated RNAs messenger (mRNA) with the top 10 rich ration of biological processes (BP) (A), cellular component (CC) (B), and molecular function (MF) (C). The abscissa is the rich ration [(DE genes of the term/all the DE genes)/(all the annotated gene of the terms/all the annotated gene)]; the Y-axis shows the terms enriched. A higher rich ration correlates with lower p-values. Circle size represents the number of enriched genes, and the color indicates the degree of enrichment, with red representing the highest degree of enrichment.

FIGURE 3

Figure 3. Functional annotation of the differentially expressed (DE) up-regulated RNAs messenger (mRNA) with the top 10 rich ration of biological processes (BP) (A), molecular function (MF) (B), and KEGG pathways (C). The abscissa is the rich ration [(DE genes of the term or pathway/all the DE genes)/(all the annotated gene of the term or pathway/all the annotated gene)]; the Y-axis shows the terms enriched. A higher rich ration correlates with lower p-values. Circle size represents the number of enriched genes, and the color indicates the degree of enrichment, with red representing the highest degree of enrichment.

The KEGG analysis showed that the down-regulated mRNAs were enriched in six pathways including extracellular matrix-receptor 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 (145 up-regulated lncRNAs, 721 down-regulated lncRNAs; Figures 1C,D). Of those DE lncRNA, 497 lncRNAs were novel lncRNAs. The predicted target mRNAs of down-regulated lncRNAs were mainly enriched in 173 BP terms, 85 CC terms and 39 MF terms (Table 6 and Supplementary Table 4). 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).

FIGURE 4

Figure 4. Functional annotation of predicted targets of differentially expressed (DE) down-regulated long non-coding RNAs (lncRNA) with the top 10 enrichment scores of biological processes (BP) (A), cellular component (CC) (B), molecular function (MF) (C), and KEGG pathways (D). The abscissa is the rich ration [(DE genes of the term or pathway/all the DE genes)/(all the annotated gene of the term or pathway/all the annotated gene)]; the Y-axis shows the terms enriched. A higher rich ration correlates with lower p-values. Circle size represents the number of enriched genes, and the color indicates the degree of enrichment, with red representing the highest degree of enrichment.

FIGURE 5

Figure 5. Functional annotation of predicted targets of differentially expressed (DE) up-regulated long non-coding RNAs (lncRNA) with the top 10 enrichment scores of biological processes (BP) (A), cellular component (CC) (B), molecular function (MF) (C), and KEGG pathways (D). The abscissa is the rich ration [(DE genes of the term or pathway/all the DE genes)/(all the annotated gene of the term or pathway/all the annotated gene)]; the Y-axis shows the terms enriched. A higher rich ration correlates with lower p-values. Circle size represents the number of enriched genes, and the color indicates the degree of enrichment, with red representing the highest degree of enrichment.

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).

The Functional Analysis of DE circRNAs

circRNAs exert their functions through host genes, and 474 DE circRNAs (247 up-regulated, 227 down-regulated) 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 down-regulated circRNA were enriched in protein processing in the endoplasmic reticulum.

FIGURE 6

Figure 6. Functional annotation of host genes of differentially expressed (DE) down-regulated circular RNAs (circRNA) with the top 10 rich ration of biological processes (BP) (A), cellular component (CC) (B), and molecular function (MF) (C). The abscissa is the rich ration. The abscissa is the rich ration [(DE genes of the term/all the DE genes)/(all the annotated gene of the term / all the annotated gene)]; the Y-axis shows the terms enriched. A higher rich ration correlates with lower p-values. Circle size represents the number of enriched genes, and the color indicates the degree of enrichment, with red representing the highest degree of enrichment.

FIGURE 7

Figure 7. Functional annotation of host genes of differentially expressed (DE) up-regulated circular RNAs (circRNA) with the top 10 rich ration of cellular component (CC) (A) and molecular function (MF) (B). The abscissa is the rich ration [(DE genes of the term /all the DE genes)/(all the annotated gene of the term / all the annotated gene)]; the Y-axis shows the terms enriched. A higher rich ration correlates with lower p-values. Circle size represents the number of enriched genes, and the color indicates the degree of enrichment, with red representing the highest degree of enrichment.

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.

FIGURE 8

Figure 8. circRNA-miRNA co-expression network based on the RNA sequencing results. Circles represent circRNAs, and green triangles represent miRNAs. The label was the name of RNA.

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).

FIGURE 9

Figure 9. Expression relationship of nine differentially expressed genes validated by qRT-PCR and RNA sequencing.

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 (Wang et al., 2015). 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) (Wang et al., 2018) 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-4-hydroxylase α 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 (Wang et al., 2019). 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 up-regulated 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.

Author Contributions

JL, LD, SL, JZ, JW, HX, and WL conceived and designed the experiments. JL, LD, YX, JC, and LW conducted the experiments. JL, HJ, YL, LD, TM, and MG analyzed the data. JL and LD wrote the manuscript. SL, HJ, and JZ revised the manuscript. All authors read and approved the final version of this manuscript.

Funding

The authors gratefully acknowledge grant support from the special fund for deep-underground medical research (Grant No. YB2018002) and 1.3.5 project for disciplines of excellence (Grant Nos. ZYJC18016 and ZYJC21048) provided by West China Hospital, Sichuan University, as well as the Sichuan International Technological Innovation Cooperation Project (Grant No. 2018HH0159), National Natural Science Foundation of China (Grant Nos. 51822403 and 81800892), and the research fund of Health Commission of Sichuan Province (Grant No. 20PJ029).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We sincerely thank Chairman Zhiliang Wu and his colleagues at Jiapigou Minerals Limited Corporation of China National Gold Group Corporation for their help providing an experimental site and their assistance in conducting research. We express our gratitude to Professor Aixiang Wu of the University of Science and Technology, Beijing for his guidance and help in searching for this experimental site.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.698046/full#supplementary-material

Footnotes

1. ^ http://ccb.jhu.edu/software/hisat2/index.shtml
2. ^ https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi

References

Antonelli, F., Belli, M., Sapora, O., Simone, G., Sorrentino, E., and Tabocchini, M. (2000). Radiation biophysics at the Gran Sasso laboratory: influence of a low background radiation environment on the adaptive response of living cells. Nucl. Phys. B Proc. Suppl. 87, 508–509. doi: 10.1016/S0920-5632(00)00735-0

CrossRef Full Text | Google Scholar

Beermann, J., Piccoli, M. T., Viereck, J., and Thum, T. (2016). Non-coding RNAs in Development and Disease: background, Mechanisms, and Therapeutic Approaches. Physiol. Rev. 96, 1297–1325. doi: 10.1152/physrev.00041.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Belli, P. (2017). Guest Editor’s Preface to the Special Issue on low background techniques. Int. J. Mod. Phys. A 32:1702001. doi: 10.1142/s0217751x17020018

CrossRef Full Text | Google Scholar

Bhutia, Y. D., Hung, S. W., Patel, B., Lovin, D., and Govindarajan, R. (2011). CNT1 Expression Influences Proliferation and Chemosensitivity in Drug-Resistant Pancreatic Cancer Cells. Cancer Res. 71, 1825–1835.

Google Scholar

Bigildeev, A. E., Zhironkina, O. A., Lubkova, O. N., and Drize, N. J. (2013). Interleukin-1 beta is an irradiation-induced stromal growth factor. Cytokine 64, 131–137. doi: 10.1016/j.cyto.2013.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Camon, E., Magrane, M., Barrell, D., Lee, V., Dimmer, E., Maslen, J., et al. (2004). The Gene Ontology Annotation (GOA) Database: sharing knowledge in Uniprot with Gene Ontology. Nucleic Acids Res. 32, D262–D266. doi: 10.1093/nar/gkh021

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, Y., Gao, X., Yang, Y., Ye, Z., Wang, E., and Dong, Z. (2018). Changing expression profiles of long non-coding RNAs, mRNAs and circular RNAs in ethylene glycol-induced kidney calculi rats. BMC Genomics 19:660. doi: 10.1186/s12864-018-5052-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Carbone, M. C., Pinto, M., Antonelli, F., Amicarelli, F., Balata, M., Belli, M., et al. (2009). The Cosmic Silence experiment: on the putative adaptive role of environmental ionizing radiation. Radiat. Environ. Biophys. 48, 189–196. doi: 10.1007/s00411-008-0208-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Castillo, H., Schoderbek, D., Dulal, S., Escobar, G., Wood, J., Nelson, R., et al. (2015). Stress induction in the bacteria Shewanella oneidensis and Deinococcus radioduransin response to below-background ionizing radiation. Int. J. Radiat. Biol. 91, 749–756. doi: 10.3109/09553002.2015.1062571

PubMed Abstract | CrossRef Full Text | Google Scholar

Castillo, H., and Smith, G. B. (2017). Below-Background Ionizing Radiation as an Environmental Cue for Bacteria. Front. Microbiol. 8:177. doi: 10.3389/fmicb.2017.00177

CrossRef Full Text | Google Scholar

Chien, C. Y., Hung, Y. J., Shieh, Y. S., Hsieh, C. H., and Lee, C. H. (2017). A novel potential biomarker for metabolic syndrome in Chinese adults: circulating protein disulfide isomerase family A, member 4. PLoS One 12:e0179963. doi: 10.1371/journal.pone.0179963

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, N. J., Gao, J. L., Liu, Y., Song, W., Zhang, Z. Y., and Gao, C. F. (2015). Label-free quantitative mass spectrometry reveals a panel of differentially expressed proteins in colorectal cancer. BioMed Res. Int. 14:365068. doi: 10.1155/2015/365068

PubMed Abstract | CrossRef Full Text | Google Scholar

Forman, H. J., Maiorino, M., and Ursini, F. (2010). Signaling Functions of Reactive Oxygen Species. Biochemistry 49, 835–842.

Google Scholar

Fratini, E., Carbone, C., Capece, D., Esposito, G., Simone, G., Tabocchini, M. A., et al. (2015). Low-radiation environment affects the development of protection mechanisms in V79 cells. Radiat. Environ. Biophys. 54, 183–194. doi: 10.1007/s00411-015-0587-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedländer, M. R., Mackowiak, S. D., Na, L., Wei, C., and Nikolaus, R. (2011). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52.

Google Scholar

Gao, Y., Wang, J., and Zhao, F. (2015). CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 16:4.

Google Scholar

Ge, J. C., Wang, Y. X., Chen, Z. B., and Chen, D. F. (2020). Integrin alpha 7 correlates with poor clinical outcomes, and it regulates cell proliferation, apoptosis and stemness via PTK2-PI3K-Akt signaling pathway in hepatocellular carcinoma. Cell. Signal. 66:109465. doi: 10.1016/j.cellsig.2019.109465

PubMed Abstract | CrossRef Full Text | Google Scholar

Griffiths-Jones, S., Saini, H. K., van Dongen, S., and Enright, A. J. (2008). miRBase: tools for microRNA genomics. Nucleic Acids Res. 36:D154–D158.

Google Scholar

Han, Y., Hong, Y., Li, L., Li, T., Zhang, Z., Wang, J., et al. (2018). A Transcriptome-Level Study Identifies Changing Expression Profiles for Ossification of the Ligamentum Flavum of the Spine. Mol. Ther. Nucleic Acids 12, 872–883. doi: 10.1016/j.omtn.2018.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacob, A., Linklater, E., Bayless, B. A., Lyons, T., and Prekeris, R. (2016). The role and regulation of Rab40b-Tks5 complex during invadopodia formation and cancer cell invasion. J. Cell Sci. 129, 4341–4353. doi: 10.1242/jcs.193904

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, R., Shen, J., Zhang, T., Liu, Q., Liao, C., Ma, H., et al. (2017). The highly expressed COL4A1 genes contributes to the proliferation and migration of the invasive ductal carcinomas. Oncotarget 8, 58172–58183. doi: 10.18632/oncotarget.17345

PubMed Abstract | CrossRef Full Text | Google Scholar

Kukleva, L. M., Shavina, N. Y., Odinokov, G. N., Oglodin, E. G., Nosov, N. Y., Vinogradova, N. A., et al. (2015). Analysis of diversity and identification of the genovariants of plague agent strains from Mongolian foci. Russ. J. Genet. 51, 238–244.

Google Scholar

Lan, S., Zheng, X., Hu, P., Xing, X., Ke, K., Wang, F., et al. (2020). Moesin facilitates metastasis of hepatocellular carcinoma cells by improving invadopodia formation and activating β-catenin/MMP9 axis. Biochem. Biophys. Res. Commun. 524, 861–868. doi: 10.1016/j.bbrc.2020.01.157

PubMed Abstract | CrossRef Full Text | Google Scholar

Landgren, H., and Curtis, M. A. (2011). Locating and labeling neural stem cells in the brain. J. Cell. Physio. 226, 1–7. doi: 10.1002/jcp.22319

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B. (2010). Aligning Short Sequencing Reads with Bowtie. Curr. Protoc. Bioinformatics 11:11.7.

Google Scholar

Li, H. (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. Cambridge: Broad Institute of Harvard and MIT.

Google Scholar

Lin, X., Zhihong, Y., Ya, Z., Dongqi, L., Fengdi, H., Yedan, L., et al. (2018). Deep RNA sequencing reveals the dynamic regulation of miRNA, lncRNAs, and mRNAs in osteosarcoma tumorigenesis and pulmonary metastasis. Cell Death Dis. 9:772. doi: 10.1038/s41419-018-0813-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Ma, T., Gao, M., Liu, Y., Liu, J., Wang, S., et al. (2020). Proteomics provides insights into the inhibition of Chinese hamster V79 cell proliferation in the deep underground environment. Sci. Rep. 10:14921. doi: 10.1038/s41598-020-71154-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Ma, T., Liu, Y., Zou, J., Gao, M., Zhang, R., et al. (2018). History, advancements, and perspective of biological research in deep-underground laboratories: a brief review. Environ. Int. 120, 207–214. doi: 10.1016/j.envint.2018.07.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550.

Google Scholar

Parsa, S., Sharifzadeh, S., Monabati, A., Seyyedi, N., Ranjbaran, R., Baghbani, M. R., et al. (2019). Overexpression of Semaphorin-3A and Semaphorin-4D in the Peripheral Blood from Newly Diagnosed Patients with Chronic Lymphocytic Leukemia. Int. J. Hematol. Oncol. Stem Cell Res. 13, 25–34.

Google Scholar

Planel, H., Soleilhavoup, J. P., Tixador, R., Richoilley, G., Conter, A., Croute, F., et al. (1987). Influence on cell proliferation of background radiation or exposure to very low, chronic gamma radiation. Health Phys. 52, 571–578. doi: 10.1097/00004032-198705000-00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, Z. Q., Tan, W. F., Zhang, B. H., and Jiang, X. Q. (2011). Chemokine CCL20 and malignant tumors. Tumor 31, 961–963.

Google Scholar

Rack, P. G., Ni, J., Payumo, A. Y., Nguyen, V., Crapster, J. A., Hovestadt, V., et al. (2014). Arhgap36-dependent activation of Gli transcription factors. Proc. Natl. Acad. Sci. U. S. A. 111, 11061–11066.

Google Scholar

Ranjith, P. G., Zhao, J., Ju, M., De Silva, R. V. S., Rathnaweera, T. D., and Bandara, A. K. M. S. (2017). Opportunities and Challenges in Deep Mining: a Brief Review. Engineering 3, 546–551. doi: 10.1016/J.ENG.2017.04.024

CrossRef Full Text | Google Scholar

Rausch, V., and Hansen, C. G. (2020). The Hippo Pathway, YAP/TAZ, and the Plasma Membrane. Trends Cell Biol. 30, 32–48. doi: 10.1016/j.tcb.2019.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Rupp, O., Macdonald, M. L., Li, S., Dhiman, H., and Lee, K. H. (2018). A reference genome of the Chinese hamster based on a hybrid assembly strategy. Biotechnol. Bioeng. 115, 2087–2100.

Google Scholar

Satta, L., Antonelli, F., Belli, M., Sapora, O., Simone, G., Sorrentino, E., et al. (2002). Influence of a low background radiation environment on biochemical and biological responses in V79 cells. Radiat. Environ. Biophys. 41, 217–224. doi: 10.1007/s00411-002-0159-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, G. B., Grof, Y., Navarrette, A., and Guilmette, R. A. (2011). Exploring biological effects of low level radiation from the other side of background. Health Phys. 100, 263–265. doi: 10.1097/HP.0b013e318208cd44

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, P., Ning, Z., Hao, X., Lin, X., Zheng, Q., and Li, S. (2019). Comparative Analysis of Whole-Transcriptome RNA Expression in MDCK Cells Infected With the H3N2 and H5N1 Canine Influenza Viruses. Front. Cell. Infect. Microbiol. 9:76. doi: 10.3389/fcimb.2019.00076

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, C., and Buolamwini, J. K. (2019). A novel RNA variant of human concentrative nucleoside transporter 1 (hCNT1) that is a potential cancer biomarker. Exp. Hematol. Oncol. 8:18. doi: 10.1186/s40164-019-0144-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Longchao, X., Hongliang, L., Ping, C., Yibin, F., Jingxuan, Z., et al. (2015). The Role of HMGB1 Signaling Pathway in the Development and Progression of Hepatocellular Carcinoma: a Review. Int. J. Mol. Sci. 16, 22527–22540. doi: 10.3390/ijms160922527

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Jin, Y., Bhandari, A., Yao, Z., Yang, F., Pan, Y., et al. (2018). Upregulated LAMB3 increases proliferation and metastasis in thyroid cancer. Onco Targets ther. 11, 37–46. doi: 10.2147/ott.S149613

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Feng, Y., Li, J., Zou, J., and Fan, L. (2019). Integrative microRNA and mRNA analysis reveals regulation of ER stress in the Pacific white shrimp Litopenaeus vannamei under acute cold stress. Comp. Biochem. Physiol. Part D Genomics Proteomics 33:100645. doi: 10.1016/j.cbd.2019.100645

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilusz, J. E., and Sharp, P. A. A. (2013). Circuitous Route to Noncoding RNA. Science 340, 440–441. doi: 10.1126/science.1238522

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, H., Gao, F., Yang, J. U., Zhang, R., Gao, M., Deng, J., et al. (2017). Novel Idea and Disruptive Technologies for the Exploration and Research of Deep Earth. Adv. Eng. Sci. 49, 1–8.

Google Scholar

Xie, H. P., Liu, J. F., Gao, M. Z., Wan, X. H., Liu, S. X., Zou, J., et al. (2018). The Research Advancement and Conception of the Deep-underground Medicine. Sichuan Da Xue Xue Bao Yi Xue ban 49, 163–168.

Google Scholar

Zeng, B., Zhou, M., Wu, H., and Xiong, Z. (2018). SPP1 promotes ovarian cancer progression via Integrin β1/FAK/AKT signaling pathway. Onco Targets Ther. 11, 1333–1343. doi: 10.2147/ott.S154215

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Bie, Q., Wu, P., Zhang, J., You, B., Shi, H., et al. (2018). PGD2/PTGDR2 Signaling Restricts the Self-Renewal and Tumorigenesis of Gastric Cancer. Stem Cells 36, 990–1003.

Google Scholar

Zhang, G., Bi, M., Li, S., Wang, Q., and Teng, D. (2018). Determination of core pathways for oral squamous cell carcinoma via the method of attract. J. Cancer Res. Ther. 14, S1029–S1034. doi: 10.4103/0973-1482.206868

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Li, B., Xu, P., and Yang, B. (2019). Integrated Whole Transcriptome Profiling and Bioinformatics Analysis for Revealing Regulatory Pathways Associated With Quercetin-Induced Apoptosis in HCT-116 Cells. Front. Pharmacol. 10:798. doi: 10.3389/fphar.2019.00798

PubMed Abstract | CrossRef Full Text | Google Scholar