Pai Zheng1
Yulin Kang2
Shuo Han1Huimin Feng1Feizai Ha2Changmao Long1Di Zhou1Guiping Hu3Zhangjian Chen1Zengmiao Wang4,5*Tiancheng Wang2Guang Jia1*- 1Department of Occupational and Environmental Health Science, School of Public Health, Peking University, Beijing, China
- 2Institute of Environmental Information, Chinese Research Academy of Environmental Sciences, Beijing, China
- 3School of Medical Science and Engineering, Beihang University, Beijing, China
- 4State Key Laboratory of Remote Sensing Science, College of Global Change and Earth System Science, Beijing Normal University, Beijing, China
- 5Department of Chemistry and Biochemistry, University of California, San Diego, La Jolla, CA, United States
Hexavalent chromium [Cr(VI)] is a well-known occupational carcinogen, but the mechanisms contributing to DNA damage and cell cycle alternation have not been fully characterized. To study the dose-response effects of Cr(VI) on transcription, we exposed BEAS-2B cells to Cr(VI) at concentrations of 0.2, 0.6, and 1.8 μmol/L for 24 h. Here, we identified 1,484 differentially expressed genes (DEGs) in our transcript profiling data, with the majority of differentially expressed transcripts being downregulated. Our results also showed that these DEGs were enriched in pathways associated with the cell cycle, including DNA replication, chromatin assembly, and DNA repair. Using the differential expressed genes related to cell cycle, a weighted gene co-expression network was constructed and a key mRNA-lncRNA regulation module was identified under a scale-free network with topological properties. Additionally, key driver analysis (KDA) was applied to the mRNA-lncRNA regulation module to identify the driver genes. The KDA revealed that ARD3 (FDR = 1.46 × 10–22), SND1 (FDR = 5.24 × 10–8), and lnc-DHX32-2:1 (FDR = 1.43 × 10–17) were particularly highlighted in the category of G2/M, G1/S, and M phases. Moreover, several genes we identified exhibited great connectivity in our causal gene network with every key driver gene, including CDK14, POLA1, lnc-NCS1-2:1, and lnc-FOXK1-4:1 (all FDR < 0.05 in those phases). Together, these results obtained using mathematical approaches and bioinformatics algorithmics might provide potential new mechanisms involved in the cytotoxicity induced by Cr.
Introduction
Inhalation exposure to hexavalent chromium [Cr(VI)] has been recognized as a significant occupational carcinogen according to the final document from the National Institute for Occupational Safety and Health (NIOSH) (NOISH, 2008). Several researchers have explored the underlying molecular mechanisms induced by Cr(VI) related to cellular transformation (Karaulov et al., 2019) and tumorigenesis (Holmes et al., 2008; Wang et al., 2019), whereas the transcriptomic responses remain elusive. Of all the possible mechanisms of carcinogenesis associated with exposure to Cr(VI), it is hard to ignore the effects of epigenetic modifications and cytogenetic damage (Rager et al., 2019), mainly the regulatory signaling pathways related to the two processes described above.
To better understand the history of transcriptome studies in the field, a brief bibliometric study was conducted to identify trends based on the frequently occurring keywords in published papers from different years. This study clearly demonstrated the changes in the main topics in this area during the last half-century. Since the first report of cancer among a population of Cr-exposed workers, Cr(VI) compounds have remained on the list of potential threats to human health. Before 1995, several studies laid a profound foundation in the toxicity field, especially regarding DNA expression under the condition of exposure to multiple metals (i.e., with nickel) (Alcedo et al., 1994) and the proposed “uptake-reduction” model, which suggests the hypothesis that molecular events of genes may induce gene expression changes in carcinogenicity (Wetterhahn and Hamilton, 1989; Dubrovskaya and Wetterhahn, 1998) and DNA damage (Standeven and Wetterhahn, 1989). From 1995 to 2000, researchers mainly focused on the toxicology of lung cells, including type II pneumocytes (Shumilla and Barchowsky, 1999) and fibroblasts (Carlisle et al., 2000), and the mechanisms of apoptosis induced by genetically programmed cell death or the effect of transcriptional inhibition (Singh et al., 1998; Röpke et al., 2000). From 2000 to 2005, Cr-induced apoptosis became the most debatable topic. Researchers studied the possible sensors or mediators involved in apoptosis, particularly the effect of ATM protein and p38 MAP kinase (Ha et al., 2003; Wakeman et al., 2005). In addition, the role of free radicals following Cr(VI)-induced DNA damage and carcinogenesis was speculated (Liu et al., 2001; Zhang et al., 2001). From 2006 to 2010, the potential mechanisms of Cr(VI) carcinogenesis in lung cells were extensively published, and the main hypothesis was related to genomic instability (Holmes et al., 2008), including microsatellite instability (Hirose et al., 2002; Takahashi et al., 2005), numerical chromosome instability (Xie et al., 2005), and consequences of the imbalance between cellular damage and repair systems (Yao et al., 2008). Then, by the year 2010, the Environmental Working Group had detected Cr-polluted drinking water in 42 states that affected 74 million Americans (Sutton, 2010), resulting in calling for a legal limit for Cr(VI) and studies investigating the environmental carcinogenicity targeted toward the digestive system (Stern, 2010; Kopec et al., 2012; Thompson et al., 2014). Additionally, the successful discovery of microRNAs (miRNA) (Hobert, 2008) provided a new perspective for gene regulation research. As a result, miRNA studies related to Cr were introduced into the field of environmental toxicology, and they aimed to elucidate the mechanisms of lung cancer induced by Cr (He et al., 2013; Li et al., 2014). Over the last 4 years, miRNA-related research has been a topic of high interest and has recently received increasing attention, especially topics related to DNA repair (Li et al., 2016) or glycolipid metabolism (Zhang et al., 2018).
Although all of these changes have occurred, studies investigating the genotoxic impact of Cr are still emerging, and Cr(VI) carcinogenicity is widely debated (Zhitkovich, 2011). In particular, the roles of non-coding RNAs in transcriptional responses during the exposure of physiological and toxicological levels are not well understood. Therefore, we conducted a toxicogenomics study using a data-driven analysis approach that aimed to outline integrated networks and to identify candidate key driver genes involved in the underlying mechanisms of cell cycle alterations after Cr exposure. Furthermore, considering the widespread involvement of long non-coding RNAs (lncRNAs) in multiple cellular functions (Mercer et al., 2009; Yao et al., 2019), including the cell cycle (Kitagawa et al., 2013) and other Cr-induced processes (Hu et al., 2019), we hypothesized that lncRNAs and coding genes might be key mediators of the responses to DNA damage by regulating the cell cycle in Cr-induced genotoxicity.
For the purpose of building a network to reveal the correlations between these RNAs, we integrated the results of several bioinformatic analysis approaches based on the expression data of RNAs. These included weighted correlation network analysis (WGCNA) to create modules according to highly correlated gene expression patterns (Langfelder and Horvath, 2008). Moreover, to pinpoint the key driver gene in these processes, we conducted weight key driver analysis (Shu et al., 2016) to extensively search for potential key elements in the regulation network (Sun et al., 2015; Bailey et al., 2018) and to detect the possible trigger genes in different phases of the cell cycle. Thus, our approach aimed to provide a candidate gene list for further research on explaining the underlying molecular mechanisms that regulate the cell cycle following exposure to relatively low concentrations of Cr(VI), and to identify the non-coding RNAs that might be novel candidate molecular targets for exposure biomarker studies.
Materials and Methods
Cell Culture and RNA Extraction
The human bronchial epithelial cell line (BEAS-2B) was purchased from American Type Culture Collection (ATCC, United States), maintaining at 37°C and 5% CO2 in a humidified incubator. Cells were cultured in Bronchial Epithelial Cell Growth Medium (BEGMTM, BulletKitTMLonza, Switzerland) supplemented with the necessary components and growth factors. The BEGM media was replaced every second day, and cells were passaged when they reached 70–80% confluency by incubation with 0.25% trypsin. In the presence of a diluted potassium dichromate stock solution (K2Cr2O7, Sigma, United States), BEAS-2B cells were seeded in six-well plates (105 cells/well) and exposed to low (0.2 μmol/L), medium (0.6 μmol/L), and high concentrations (1.8 μmol/L) of Cr(VI) for 24 h. A control group was established under the same conditions as the exposure groups. Each sample containing approximately 1 × 107 cells was disrupted in buffer RLT (Qiagen, United States) for RNA isolation. Total RNA was isolated using the miRNeasy Mini Kit (Cat#217004, QIAGEN, GmBH, Germany) according to the manufacturer’s recommended guidelines, and the RIN number was determined to analyze RNA integrity using an Agilent Bioanalyzer 2,100 spectrophotometer (Agilent Technologies, Santa Clara, CA, United States). To assess the purity of RNA, a NanoDrop 2000c (Thermo Fisher Scientific, United States) and UVP Imaging System were used to measure the 260/280 ratios. Samples with a RIN number > 7 and 260/280 ratio in the range of ∼2.0 were considered qualified samples.
RNA Isolation and RNA Microarray
Total RNA was amplified and labeled using a Low Input Quick Amp Labeling Kit, One-Color (Cat.# 5190-2305, Agilent Technologies). Each slide was hybridized with 1.65 μg of Cy3-labeled cRNA with a Gene Expression Hybridization Kit (Cat.# 5188-5242, Agilent Technologies), maintaining in a Hybridization Oven (Cat.# G2545A, Agilent technologies) for 17 h hybridization. Then, slides were washed with the Gene Expression Wash Buffer Kit (Cat.# 5188-5327, Agilent Technologies) and scanned by an Agilent Microarray Scanner (Cat#G2565CA, Agilent Technologies). Finally, data were extracted with Feature Extraction software 10.7 (Agilent Technologies). Raw data were normalized with the quantile algorithm and limma packages in R.
Differential Expression Analysis
The normalization of gene expression and differentially expressed (DE) gene analysis were performed by functions in R package limma (Ritchie et al., 2015). Briefly, the function lmFit was used to fit a linear model to estimate the variability in the data. Using the function eBayes in limma, the significance of differences in the variance of gene expression across biological replicates for each gene were calculated using empirical Bayes moderated t-statistics tests. For the multiple hypothesis testing correction, the false discovery rate (FDR) was applied. The fold change in logarithm form between the samples from treated group and the samples from control group was also calculated for each gene (Smyth, 2005). Finally, genes were defined as DE if FDRs were below 0.05 and |log2FC| ≥ 1.
Cell Cycle-Related Genes
After analyzing the pathways, significantly expressed pathways were identified with the cutoff of adjusted P < 0.05. From our results, we observed that significant pathways were related to biological processes involved in the regulation of the cell cycle, which had 11 related pathways. Furthermore, we extracted all genes with the GO annotations of the cell cycle (GO:0007049) from AmiGO 2 (Carbon et al., 2009) and withdrew the DE genes in our database. At last, a linear regression was performed with the expression of each gene. We used the 12 values in 4 dose groups (each dose had triplicated biological samples) as the dependent variable, and the dose 0–1.8 μmol/L corresponding to each sample as the independent variables. We took the linear regression slope to represent expression-dose-depend relationship for each gene and filtered the genes with a threshold of the median coefficient from the whole gene set. The g:Profiler (Reimand et al., 2016) (database built on 2020-03-07) was used to conduct GO enrichment analysis based on the DE gene list identified from the high dose group vs. with the control group comparison. We also visualized the g:GOSt (database built on 2020-03-07) enrichment results from different annotation resources, including GO, KEGG, REAC, and TF.
Co-expression Network Detection Using WGCNA
WGCNA (Langfelder and Horvath, 2008) was used to identify the co-expression module with the selected cell cycle-related coding RNAs. WGCNA is a guilt-by-association approach for constructing networks and module detection. We computed a correlation raised to a power as a soft thresholding between every pair of RNAs to amplify the disparity and transformed the result into an adjacency matrix. Then, the blockwiseModules function was used to compute topological overlap matrix (TOM) dissimilarity (Yip and Horvath, 2007) between genes, and a hierarchical clustering gene dendrogram was constructed. According to the standard of dynamicTreeCut, modules whose eigengenes are highly correlated are merged with the threshold of 0.2. At last, we applied the module preservation analysis to test whether the features of each module were preserved in an alternative set of samples.
The preservation analysis is to estimate the differences between the observation in our gene expression and random situations by permutation (nPermutations = 10). Furthermore, Z-summary statistical analysis was performed as a general summary of all of the different statistics used, and scores > 10 were considered preserved (Langfelder et al., 2011).
Key Driver Analysis
To conduct key driver analysis, we used the R package KDA (v1.14.0) (Shu et al., 2016). The purpose of KDA analysis was to analyze the detailed interactions between coding genes and lncRNAs to select genes that were over-represented in the regulation of the cell cycle. The package first required a sub-network file and gene list of interest as input files. Then, the enrichment in the k-step downstream neighborhood for the target gene list was assessed in the sub-network. In this study, we used a list of cell-cycle related genestarget list as the input file. We computed the co-expression pairs between targeted coding genes and all lncRNAs in our experimental study as the network file. We divided these coding genes into three group according to the key events in distinct phases, for example the checkpoint which regulated the transition of phases. Finally, we combined the lists with the name of cell cycle class as the input lists for KDA detection. The KDA results were visualized using Cytoscape (v3.7.2).
Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
For quantitative real-time PCR (RT-qPCR) analyses, total RNA was isolated from the control, middle, and high dose groups of the samples as described above, and those groups along with three independent replications of the biological experiments were analyzed. We selected and verified the key driver genes and genes in the network (Figure R1) for their biological potential in explaining the Cr(VI) cytotoxicity with the following criteria: (1) the genes participate in regulating at least two stages of cell cycle progress in the subnetwork; (2) the genes were reported to be involved in the biological functions, such as cell cycle arrest, carcinogenicity, and lung cancer, according to the literature. According to these criteria, we selected three groups of genes, in the (1) PARD3 group: ENST00000607815 and PARD3B, (2) SND1group: TEX14, (3) genes in both PARD3, and SND1 group:lnc-NCS1-2:1, FBXO6, CDK14, lnc-DHX32-2:1, lnc-FOXK1-4:1, POLA1, NR_002579, lnc-C3orf14-1:1 as the target genes for further validation.
Primer sequences were designed using Primer Express 3.0.1 designer software and then verified with NCBI Primer-BLAST software to confirm specific recognition of the target lncRNAs and mRNAs. qPCR was performed using a CFX96 Touch Real-Time PCR Detection System (Bio-Rad Laboratories, United States) to determine the expression levels of selected genes, and the results were presented as the mean value for the three wells. Data were calculated using the 2–△△Ct with glycerol-3-phosphate dehydrogenase (GAPDH) as the endogenous control.
Results
DE Genes Between Different Groups
In total, we collected 12 samples, 3 from the low-dosage group, 3 from the medium-dosage group, 3 from the high-dosage group and 3 from the control group. Genome wide analysis of RNA expression by microarray was used to examine the transcription changes between the Cr(VI) treated group and the control group. To fully characterize the regulation pattern between mRNA and lncRNA, 18833 mRNA and 68104 lncRNA were considered in this study. Differential expression analysis was conducted for each Cr(VI) dosage vs. control comparison. The largest differences in gene expression were evident in the high dose group (1.8 μmol/L) with log2 fold-change [log2(FC)] values ranging from −3.76 to 4.52, including 646 upregulated and 771 downregulated genes (Figures 1B,D). A relatively lower number of differentially expressed genes was observed in the middle (129) (Figures 1A,C) and low dose (17) groups compared with normal cultured cells. It was evident that there were more upregulated and downregulated lncRNAs in the 1.8 μmol/L group compared with the other groups (Figure 1C), and DE lncRNA from high-dose accounts for 89.3% of all DE lncRNAs across three groups. From the histograms of fold change profile for mRNA and lncRNA in the two groups, we found that the high dose group had more differentially expressed transcripts with log2FC values above 0.5 and below −0.5 than the middle dose group (Supplementary Figure 1A), whereas the differences in lncRNAs were not as clear (Supplementary Figure 1B). Taken together, the high level of hexavalent chromium induces the dramatic changes in transcription profiling. To get a unified DE gene set, we merged the DE genes from three different groups and 1,484 mRNAs and lncRNAs were left (Supplementary Table 1), while the detail information was listed in Supplementary Tables 2, 3.
Figure 1. Differentially expressed genes in 0.6 and 1.8 μmol/L Cr(VI) groups. (A,B) Differentially expressed mRNAs in 0.6 and 1.8 μmol/L Cr(VI) groups. (C,D) Differentially expressed lncRNAs in 0.6 and 1.8 μmol/L Cr(VI) groups.
Enrichment of Protein-Coding Genes
To further explore the potential biological processes induced by Cr, enrichment analyses using DEGs databases were conducted using g:Profiler. As shown in Figure 2, DNA replication (FDR = 5.97 × 10–18) and chromatin organization (FDR = 2.70 × 10–12) from GO biological process database, cell cycle-related pathways from the Reactome Pathway Database (REAC) were identified with almost the highest statistically significant confidence level (all FDR < 0.05), suggesting that these processes might substantially change after exposure to the Cr(VI) (Figure 2B). Moreover, after analyzing the most altered pathways in replicating cells in detail, we identified the pathway involving SIRT1 (FDR = 4.56 × 10–22), a component of the Energy-dependent Nucleolar Silencing Complex (eNoSC), which may serve as an important integration point in cell cycle regulation (Bouras et al., 2005). Other pathways, including DNA double-strand break repair (FDR = 0.017) and mitotic cycle alternation (FDR = 0.046), also suggested a link between the modification of the cell cycle and exposure to Cr(VI). Therefore, we selected genes with strong enrichment with specific processes as the target biological reactions for analysis.
Figure 2. Enrichment analysis of differentially expressed genes in the high-dose group. The total transcript level correlates to (A) cell cycle processes and (B) detailed information for the chosen pathways.
Selection of Cell Cycle-Related Genes
To determine the gene list for further analysis of transcript profiles and co-expression with lncRNAs, we select cycle-related genes based on three criterion: (1) genes annotated with the function of cell cycle, (2) DE genes, (3) existing linear relationship. We first downloaded the coding gene data and their information regarding GO class and evidence code from AmiGO 2 with the annotation of the cell cycle (GO:0007049) as the background list. Subsequently, 64 gene transcripts were matched with our differentially expressed genes, and the list in AmiGO 2 of 918 genes was related to 1,952 items with target processes. These genes were then screened according to the criterion for dose-expression relationships. To accurately identify the genes related to the effect of Cr(VI) treatment on cell proliferation, we filtered genes with absolute regression coefficients > 0.4 (close to the value of Q1 and Q3 in the cell cycle). Finally, 35 genes were selected as targeted mRNAs for the construction of co-expression modules with lncRNAs.
The expression levels in the four treatment groups from the mRNA profiles and evidence codes of these genes in the cell cycle are shown in Figure 3. As shown in Figure 3A, genes were clustered into four groups, even with similar pathways in the process of regulation of the cell cycle in all clusters, we still found the cluster-specific pathway, indicating that Cr(VI) may affect the regulatory functions in different aspects. In particular, cluster one was enriched in the function of cell division (P = 1.16 × 10–6), whereas cluster two focused on mitotic cell cycle phase transition (P = 6.11 × 10–6), especially the G2/M transition of the mitotic cell cycle (P = 2.30 × 10–4). The last cluster with 9 genes, including IL1A and IL1B, was associated with the regulation of cell cycle processes (P = 3.76 × 10–5). Additionally, as shown in Figure 3B, 17 of 35 gene annotations were derived from laboratory experiments, suggesting that the transcripts were closely related to the process of cell cycle regulation.
Figure 3. Expression levels and evidence codes of target cell cycle genes. (A) The expression level of the cell cycle-related genes in different groups and samples and (B) the gene ontology evidence codes of the selected genes in the cell cycle.
Co-expression Module With Cell Cycle-Associated Protein-Coding Genes and lncRNAs
To investigate lncRNAs that were highly connected to a given set of interesting protein-coding genes, WGCNA was employed. WGCNA is R software package to perform network construction, module detection, and calculation of topological properties based on the guilt-by-association strategy. The WGCNA network was generated using the selected 35 mRNA and all lncRNA profiles from the different dose groups. We selected 16 as the soft threshold based on the results of scale-free topology and mean connectivity for the construction of the network. All genes were clustered into a unique module and the cluster dendrogram shown in Supplementary Figures 2, 3. Altogether, we identified six modules that contained target cell cycle protein-coding genes from all 185 co-expression modules. Modules were color-coded as brown, plum1, black, gray60, magenta, and tan4 with 25, 5, 2, 1, 1, and 1 mRNA, respectively. A preservation test was applied to confirm the reliability and sensitivity of the results using a WGCNA function (modulePreservation). Briefly, to evaluate whether a module was conserved or not, the Zsummary (Z-score) was calculated. All six modules were highly preserved with Z-scores > 30, whereas the brown (Zsummary = 68) was chosen for obtaining the highest number of transcripts (Figure 4). Therefore, modules colored brown containing 25 mRNAs and 3,565 lncRNAs were regarded as highly representative modules according to the expression patterns of these genes, and served as a cell cycle-regulated network for further analysis.
Figure 4. Module preservation test for co-expression cell cycle mRNA data sets. Circles in both figure represent different co-expression modules. The horizontal and vertical axes represent gene number and Zsummary values for each module with random selection, respectively. The green dotted line indicates Zsummary = 60 (Zsummary > 10 indicates high preservation). (A) Preservation median rank of these modules. (B) Preservation Zsummary. Brown modules were indicated by arrows.
Key Driver Analysis
Previous analysis recommended that genes Within the brown module Were the key cell cycle regulatory gene set. We performed key driver analysis (KDA) for these target gene co-expression modules using our previously published method (Shu et al., 2016). Key driver genes or key regulatory components were defined as the subnetworks with respect to various biological contexts, when compare to the other random genes in the network. The KDA requires a gene list (generally a biological-associated gene list) and a gene network as input files to identify the key genes. Under different cell cycle phases, we first divided 25 mRNA in the brown module Into three categories (G1/S, G2/S, and M) according to their GO subcategory annotation (with the highest level of evidence code in each gene), which were connected to the biological regulatory events in cell cycle processes. Then, Pearson’s correlation coefficient was calculated to each selected protein-coding gene and all differentially expressed lncRNAs. Pairs with absolute values of Pearson’s correlation coefficients ≥ 0.90 were selected as the network file.
Finally, three top key driver genes, including PARD3 (FDR = 1.46 × 10–22), SND1 (FDR = 5.24 × 10–8), and lnc-DHX32-2:1 (FDR = 1.43 × 10–17), were particularly highlighted in the category of G2/M, G1/S, and M phases and mathematically identified as the causal modulators of the regulatory state of the functionally relevant gene group based on prior knowledge (Figure 5). Moreover, several genes showed great connectivity in our causal gene network and interrelated with every single key driver gene, including CDK14, POLA1, lnc-NCS1-2:1, and lnc-FOXK1-4:1 (all FDR < 0.05 in those phases). Furthermore, lnc-DHX32-2:1 and PARD3B existed in three categories in our list, illustrating their potential relationship within the probabilistic causal gene network, which might suggest their complex cellular context. Additionally, some significantly differentially expressed genes existed in the network, such as FBXO6 and ENST00000607815. It was also revealed that most genes might participated in multiple phases regulation, especially in the progression G1/S and G2/M phases.
Figure 5. Key driver gene subnetworks with the phases of the cell cycle. Red, blue, and yellow circles indicate G1/S, G2/M, and M phases, respectively. Diamonds represent the key driver genes.
Validation by Real-Time RT-PCR
To determine whether the expression patterns of these genes could be recapitulated, we selected ENST00000607815, lnc-NCS1-2:1, FBXO6, CDK14, lnc-DHX32-2:1, lnc-FOXK1-4:1, PARD3, PARD3B, POLA1, SND1, NR_002579, TEX14, lnc-C3orf14-1:1 for real-time RT-PCR analysis (Figure 6). The results found 11 different expressed genes in our network, including the key driver genes SND1 and PARD3, which had significantly change between control and 1.8 μmol/L group. Eight genes which both connected with PARD3 and SND1, including lnc-NCS1-2:1, FBXO6, CDK14, lnc-DHX32-2:1, lnc-FOXK1-4:1, POLA1, NR_002579 and ENST00000607815 could be divided into three groups, (1) six of them located in G1/S and G2/M phase sub-nodes, (2) lnc-DHX32-2:1 related with all cell cycle phases change, and (3) FBXO6 was in G2/M or M phase node. Among these genes, lnc-FOXK1-4:1, POLA1 and ENST00000607815 had the most significant change in Cr group, indicating the potential effect on G1/S and G2/M phases.
Figure 6. Results of RT-PCR analysis of 11 genes in the network. ** indicates significantly different at 0.05, *** indicates significantly different at 0.001.
Discussion
The present study provided a regulatory network containing protein-coding genes and lncRNAs, including PARD3, SND1, POLA1, FBXO6, lnc-NCS1-2:1, lnc-FOXK1-4:1, lnc-DHX32-2:1, and ENST00000607815, which revealed the key regulatory mechanisms in different cell cycle phases after exposure to Cr(VI). Because our in vitro experiment was based on relatively low doses of Cr, it can be referenced as a common environmental exposure level that mimics the damage on human respiratory epithelial cells. Our results regarding the cell cycle were consistent with previous findings that Cr exposure could induce cell cycle changes by regulating checkpoint pathways to control the order and timing of cell cycle phase transition (Alcedo et al., 1994; Chiu et al., 2010; Nickens et al., 2010; Proctor et al., 2014; Rager et al., 2019), even under very low concentrations (Ha et al., 2004; Wakeman et al., 2004). In particular, the bioinformatics approach we used provided us with potential candidate lncRNAs that might play important roles in Cr(VI) exposure via the regulation of different cell cycle phases.
In the past two decades, the number of papers related to the cell cycle has rapidly increased (Kastan and Bartek, 2004; Malumbres and Barbacid, 2009) because of the complex and diverse processes involved in the cellular responses to DNA damage. Both endogenous and exogenous sources that cause DNA damage are considered major contributors in the development progress of human cancers, and thus it is reasonable to speculate that defects in cell cycle genes have critical significance in increasing the respiratory cancer risk for occupational exposure to Cr(VI) (Seidler et al., 2013). Previous transcriptomic-based studies have indicated that genes relevant to cell proliferation and DNA repair showed differential expression patterns after Cr(VI) exposure. Among these pathways, the effects of inhibitors of ATM activation have received the most attention because of the observed toxicological importance in activation of the DNA damage-responsive kinase ATM after Cr(VI) exposure (Luczak et al., 2016).
However, this pathway did not appear in our experimental study, which may be because of the threshold of checkpoint activation or DNA-damage response (Liang et al., 2014). Our results offered novel candidate regulatory genes using mathematical algorithms and cells treated with relatively low doses (< 2 μmol/L) of Cr(VI), which altered the cell cycle phases in a concentration-dependent manner. The high level of complexity in cell cycle regulation after exposure to potassium dichromate provides a great opportunity to discover novel factors in several signaling and response pathways to address the specific nature of cell damage. To minimize the adverse effects of DNA damage situation, DNA repair is a mechanism that allows cells to properly repair these defects. This is exemplified in the results of a study involving PARD3 and PARD3B, which encode proteins belonging to the Par-3 family of cell polarity regulators, that play key roles in asymmetrical cell division and polarized cell growth. Studies showed that Par-3 complexes could regulate DNA-PK directly or via Ku70 which could affect the progress of DNA double-strand break repair (Fang et al., 2007), thus indicated that Par-3 enhanced NHEJ and HR pathway required for efficient DNA repair in G2 and M phases. Besides, the loss or attenuation of epithelial polarity is a hallmark of epithelial to mesenchymal transformation (EMT) (Thiery, 2003), indicating the possibility of regulating cell cycle and enhancing the genetic instability by PARD3 or PARD3B (Supplementary Figures 2, 3). Together, this suggests the possibility that Par-3 is involved in cancer development and progression in the G1/S or G2/M phase. Moreover, the gene SND1 (staphylococcal nuclease and tudor domain containing 1) in the center of the network mediated miRNA decay of both protein-free and AGO2-loaded miRNAs and also regulated mRNAs involved in G1-to-S phase transition (Elbarbary et al., 2017). Extensive research also supports the conclusion that SND1 is an oncoprotein in a variety of cancers involving multiple processes (Jariwala et al., 2017; Xin et al., 2019) and that it also acts as an essential effector in promoting EMT in cancer (Xin et al., 2019; He et al., 2020). In addition, FBXO6 was found to closely interact with Chk1 (Tu et al., 2017), which indirectly affected the ATR-Chk1 signal axis by various kinds of DNA damage insults, including replication stress, inter-strand cross-link (ICL), and DSBs (Jazayeri et al., 2006). Similarly, a study demonstrated that FBXO6 correlated with the survival of non-small cell lung cancer (NSCLC) patients (Cai et al., 2019). CDK14 plays fundamental role in regulating the G2/M phase of the cell cycle by mediating the phosphorylation of LRP5/6 in the Wnt signaling pathway (Wang et al., 2016). Correspondingly, CDK14 played an important role in lung tumorigenesis. For example, the cigarette smoke-induced downregulation of CDK14 in lung cells correlated with β-catenin levels, suggesting impaired Wnt signaling (Pollack et al., 2015). Another study also suggested that CDK14 is regulated by SNHG15 by competitively sponging miR-486, which contributed to NSCLC tumorigenesis (Jin et al., 2018).
As RNAs do not function alone, we conducted WGCNA and probabilistic causal gene analyses to model underlying molecular relationships and causal gene connections, and used mathematical approaches and key driver analysis (KDA) to prioritize key drivers of the genotoxicity of Cr(VI). WGNCA analysis was conducted to narrow down the extent of related lncRNAs with our selected genes with a dose-response relationship, and KDA was used to identify the key regulator within the co-expression network lncRNAs, including ENST00000607815, lnc-NCS1-2:1, lnc-DHX32-2:1, and lnc-FOXK1-4:1, were verified as the important lncRNAs in our network. One of the main regulators of lnc-DHX32-2:1, which targeted for gene ADAM12, has been implicated in a variety of biological processes, including lung cancer and the development of giant cell tumors, and positively involved in the regulation of the MAPK/ERK pathway. Moreover, the MAPK/ERK pathway plays an important role in integrating external signals to promote the transition of G1 to S phase, while ERK activation downstream of mitogen-induced Ras signaling is sufficient to alleviate cell cycle arrest and allow cells to progress to the S-phase (Chambard et al., 2007). This is consistent with a previous study in which Cr(VI) activated ERK signaling pathways (Chuang et al., 2000), although no cytotoxicity effects were observed. Indeed, our analysis identified lncRNAs as interesting candidates involved in cell-cycle gene regulation, but all guilt-by-association approaches must be treated with caution and recognized as hypothesis-generating research. Additionally, all exploratory data analysis techniques require extensive targeted studies to confirm suggested molecular networks and the potential mechanisms.
Finally, our study revealed that cells exhibited cytotoxicity even in the low-dose exposure group. The choice of relatively low concentrations based on the conclusion of previous studies showed that BEAS-2B cells exhibited increased proliferation in the presence of 1 μM (Nickens et al., 2010; Cerveira et al., 2014), whereas other studies showed slightly decreased proliferation following treatment with low micromolar range concentrations after stimulating human bronchial epithelial cells (Abreu et al., 2018; Hu et al., 2019). The 2 μM concentration was considered to cause significant cytotoxic effects and alterations of the cell cycle. Therefore, increasing the sample size and number of doses would allow for more convincing conclusions. Our research design follows the principle of Toxicity Testing in the 21st Century, which emphasized the usage of cell lines, computational modeling, and bioinformatics approaches to explain cell response progress (Krewski et al., 2010; Seidler et al., 2013). This research also implied that the DNA repair or immune defense mechanisms might involve greatly to cope with the toxicity in relatively low level, which is different from the proved results, and can help assess the dose-response characteristics in exploring the perturbation of the adverse pathways (Krewski et al., 2010).
The current study had several limitations. First, the use of single omics technology reveals the cellular transcription responses patter on Cr, which might have limitation in explaining the results in the aspect of functional outcomes. Subsequently, the transcriptome experiment conducted by microarrays with high background levels owing to cross-hybridization confines the accuracy of expression estimation, particularly for the transcripts in low abundance (Zecevic et al., 2009). Finally, this study based on a single cell line that cannot depict the whole picture of complicated toxic responses in human pulmonary bronchial. Thus, further research on different human cell as well as the in vivo approaches, such as cell cycle checkpoint and related protein analysis, that aiming to test the regulation of cyclin-dependent kinases secretion with knock-out amplification these genes can offer evidence on how they activate or inhibit cell cycle phases, which can be used to explain the carcinogenicity and Cr(VI) toxicity. Another potential limitation is that repeated experiments under different passages and freeze downs are needed for further analysis. In the future, we will apply this method to other cell lines, for example, HBECs and hTERT-immortalized lung Cells. With the development on 3D-cell-culture model in vitro (Gunness et al., 2013; Amann et al., 2014). Our further research will try to analysis the transcriptome character in different condition with the 3D Cultures method to construct lung model (Rayner et al., 2019). Nonetheless, further integrated proteomics, epigenomics or the RNA modifications analysis on multi-omics level and time series studies with RNA-seq methods are needed to confirm the mechanisms of toxicity involved in the cell cycle alterations induced by Cr(VI).
Conclusion
We designed a novel computational workflow that showed the toxicity effects of Cr(VI) on the cell cycle. We identified gene networks and candidate lncRNAs and mRNAs in relatively low dose range suggesting that the gene sets could provide a clue for toxic responses in cell cycle regulation to environmental and occupational Cr(VI) exposures.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/reefur/Chromium-transcriptome/blob/main/database_cr.xlsx.
Author Contributions
PZ and ZW participated in the study design, statistical and bioinformatics analysis, and manuscript preparation. SH, HF, and FH conducted the experiment. DZ and CL provided assistance with gene set network analysis. GH, ZC, ZW, TW, and GJ reviewed the figures and manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was funded by the National Natural Science Foundation of China (81872598, 81673118, and 82073616) and Beijing Natural Science Foundation (JQ18025).
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.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.597803/full#supplementary-material
References
Abreu, P. L., Cunha-Oliveira, T., Ferreira, L. M. R., and Urbano, A. M. (2018). Hexavalent chromium, a lung carcinogen, confers resistance to thermal stress and interferes with heat shock protein expression in human bronchial epithelial cells. Biometals 31, 477–487. doi: 10.1007/s10534-018-0093-7
Alcedo, J. A., Misra, M., Hamilton, J. W., and Wetterhahn, K. E. (1994). The genotoxic carcinogen chromium(VI) alters the metal-inducible expression but not the basal expression of the metallothionein gene in vivo. Carcinogenesis 15, 1089–1092. doi: 10.1093/carcin/15.5.1089
Amann, A., Zwierzina, M., Gamerith, G., Bitsche, M., Huber, J. M., Vogel, G. F., et al. (2014). Development of an innovative 3D cell culture system to study tumour-stroma interactions in non-small cell lung cancer cells. PLoS One 9:e92511. doi: 10.1371/journal.pone.0092511
Bailey, M. H., Tokheim, C., Porta-Pardo, E., Sengupta, S., Bertrand, D., Weerasinghe, A., et al. (2018). Comprehensive characterization of cancer driver genes and mutations. Cell 173, 371–385.
Bouras, T., Fu, M., Sauve, A. A., Wang, F., Quong, A. A., Perkins, N. D., et al. (2005). SIRT1 deacetylation and repression of p300 involves lysine residues 1020/1024 within the cell cycle regulatory domain 1. J. Biol. Chem. 280, 10264–10276. doi: 10.1074/jbc.m408748200
Cai, L., Li, J., Zhao, J., Guo, Y., Xie, M., Zhang, X., et al. (2019). Fbxo6 confers drug-sensitization to cisplatin via inhibiting the activation of Chk1 in non-small cell lung cancer. FEBS Lett. 593, 1827–1836. doi: 10.1002/1873-3468.13461
Carbon, S., Ireland, A., Mungall, C. J., Shu, S., Marshall, B., Lewis, S., et al. (2009). AmiGO: online access to ontology and annotation data. Bioinformatics 25, 288–289. doi: 10.1093/bioinformatics/btn615
Carlisle, D. L., Pritchard, D. E., Singh, J., Owens, B. M., Blankenship, L. J., Orenstein, J. M., et al. (2000). Apoptosis and P53 induction in human lung fibroblasts exposed to chromium (VI): effect of ascorbate and tocopherol. Toxicol. Sci. 55, 60–68. doi: 10.1093/toxsci/55.1.60
Cerveira, J. F., Sánchez-Aragó, M., Urbano, A. M., and Cuezva, J. M. (2014). Short-term exposure of nontumorigenic human bronchial epithelial cells to carcinogenic chromium(VI) compromises their respiratory capacity and alters their bioenergetic signature. FEBS Open Bio. 4, 594–601. doi: 10.1016/j.fob.2014.06.006
Chambard, J.-C., Lefloch, R., Pouysségur, J., and Lenormand, P. (2007). ERK implication in cell cycle regulation. Biochim Biophys. Acta 1773, 1299–1310. doi: 10.1016/j.bbamcr.2006.11.010
Chiu, A., Shi, X., Lee, W., Hill, R., Wakeman, T., Katz, A., et al. (2010). Review of chromium (VI) apoptosis, cell-cycle-arrest, and carcinogenesis. J. Environ. Sci. Health C Environ. Carcinog Ecotoxicol. Rev. 28, 188–230.
Chuang, S. M., Liou, G. Y., and Yang, J. L. (2000). Activation of JNK, p38 and ERK mitogen-activated protein kinases by chromium(VI) is mediated through oxidative stress but does not affect cytotoxicity. Carcinogenesis 21, 1491–1500. doi: 10.1093/carcin/21.5.491
Dubrovskaya, V. A., and Wetterhahn, K. E. (1998). Effects of Cr(VI) on the expression of the oxidative stress genes in human lung cells. Carcinogenesis 19, 1401–1407. doi: 10.1093/carcin/19.8.1401
Elbarbary, R. A., Miyoshi, K., Myers, J. R., Du, P., Ashton, J. M., Tian, B., et al. (2017). Tudor-SN-mediated endonucleolytic decay of human cell microRNAs promotes G/S phase transition. Science 356, 859–862. doi: 10.1126/science.aai9372
Fang, L., Wang, Y., Du, D., Yang, G., Tak Kwok, T., Kai Kong, S., et al. (2007). Cell polarity protein Par3 complexes with DNA-PK via Ku70 and regulates DNA double-strand break repair. Cell Res. 17, 100–116. doi: 10.1038/sj.cr.7310145
Gunness, P., Mueller, D., Shevchenko, V., Heinzle, E., Ingelman-Sundberg, M., and Noor, F. J. T. S. (2013). 3D organotypic cultures of human HepaRG cells: a tool for in vitro toxicity studies. Toxicol. Sci. 133, 67–78. doi: 10.1093/toxsci/kft021
Ha, L., Ceryak, S., and Patierno, S. R. (2003). Chromium (VI) activates ataxia telangiectasia mutated (ATM) protein. Requirement of ATM for both apoptosis and recovery from terminal growth arrest. J. Biol. Chem. 278, 17885–17894. doi: 10.1074/jbc.m210560200
Ha, L., Ceryak, S., and Patierno, S. R. (2004). Generation of S phase-dependent DNA double-strand breaks by Cr(VI) exposure: involvement of ATM in Cr(VI) induction of gamma-H2AX. Carcinogenesis 25, 2265–2274. doi: 10.1093/carcin/bgh242
He, A., He, S., Huang, C., Chen, Z., Wu, Y., Gong, Y., et al. (2020). MTDH promotes metastasis of clear cell renal cell carcinoma by activating SND1-mediated ERK signaling and epithelial-mesenchymal transition. Aging (Albany NY) 12, 1465–1487. doi: 10.18632/aging.102694
He, J., Qian, X., Carpenter, R., Xu, Q., Wang, L., Qi, Y., et al. (2013). Repression of miR-143 mediates Cr (VI)–Induced tumor angiogenesis via IGF-IR/IRS1/ERK/IL-8 Pathway. Toxicol. Sci. 134, 26–38. doi: 10.1093/toxsci/kft101
Hirose, T., Kondo, K., Takahashi, Y., Ishikura, H., Fujino, H., Tsuyuguchi, M., et al. (2002). Frequent microsatellite instability in lung cancer from chromate-exposed workers. Mol. Carcinogenesis 33, 172–180. doi: 10.1002/mc.10035
Hobert, O. J. S. (2008). Gene regulation by transcription factors and microRNAs. Science 319, 1785–1786. doi: 10.1126/science.1151651
Holmes, A. L., Wise, S. S., and Wise, J. P. (2008). Carcinogenicity of hexavalent chromium. Ind. J. Med. Res. 128, 353–372.
Hu, G., Feng, H., Long, C., Zhou, D., Li, P., Gao, X., et al. (2019). LncRNA expression profiling and its relationship with DNA damage in Cr(VI)-treated 16HBE cells. Sci. Total Environ. 655, 622–632. doi: 10.1016/j.scitotenv.2018.11.074
Jariwala, N., Rajasekaran, D., Mendoza, R. G., Shen, X.-N., Siddiq, A., Akiel, M. A., et al. (2017). Oncogenic role of SND1 in development and progression of hepatocellular carcinoma. Cancer Res. 77, 3306–3316. doi: 10.1158/0008-5472.can-17-0298
Jazayeri, A., Falck, J., Lukas, C., Bartek, J., Smith, G. C. M., Lukas, J., et al. (2006). ATM- and cell cycle-dependent regulation of ATR in response to DNA double-strand breaks. Nat. Cell Biol. 8, 37–45. doi: 10.1038/ncb1337
Jin, B., Jin, H., Wu, H.-B., Xu, J.-J., and Li, B. (2018). Long non-coding RNA SNHG15 promotes CDK14 expression via miR-486 to accelerate non-small cell lung cancer cells progression and metastasis. J. Cell Physiol. 233, 7164–7172. doi: 10.1002/jcp.26543
Karaulov, A. V., Renieri, E. A., Smolyagin, A. I., Mikhaylova, I. V., Stadnikov, A. A., Begun, D. N., et al. (2019). Long-term effects of chromium on morphological and immunological parameters of Wistar rats. Food Chem. Toxicol. 133:110748. doi: 10.1016/j.fct.2019.110748
Kitagawa, M., Kitagawa, K., Kotake, Y., Niida, H., and Ohhata, T. (2013). Cell cycle regulation by long non-coding RNAs. Cell Mol. Life Sci. 70, 4785–4794. doi: 10.1007/s00018-013-1423-0
Kopec, A. K., Kim, S., Forgacs, A. L., Zacharewski, T. R., Proctor, D. M., Harris, M. A., et al. (2012). Genome-wide gene expression effects in B6C3F1 mouse intestinal epithelia following 7 and 90 days of exposure to hexavalent chromium in drinking water. Toxicol. Appl. Pharmacol. 259, 13–26. doi: 10.1016/j.taap.2011.11.012
Krewski, D., Acosta, D., Andersen, M., Anderson, H., Bailar, J. C., Boekelheide, K., et al. (2010). Toxicity testing in the 21st century: a vision and a strategy. J. Toxicol. Environ. Health B Crit. Rev. 13, 51–138.
Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 9:559. doi: 10.1186/1471-2105-9-559
Langfelder, P., Luo, R., Oldham, M. C., and Horvath, S. (2011). Is my network module preserved and reproducible? PLoS Comput. Biol. 7:e1001057. doi: 10.1371/journal.pcbi.1001057
Li, Y., Hu, G., Li, P., Tang, S., Zhang, J., and Jia, G. J. T. (2016). miR-3940-5p enhances homologous recombination after DSB in Cr (VI) exposed 16HBE cell. Toxicology 344, 1–6. doi: 10.1016/j.tox.2016.02.003
Li, Y., Li, P., Yu, S., Zhang, J., Wang, T., and Jia, G. J. T. (2014). miR-3940-5p associated with genetic damage in workers exposed to hexavalent chromium. Toxicol. Lett. 229, 319–326. doi: 10.1016/j.toxlet.2014.06.033
Liang, H., Esposito, A., De, S., Ber, S., Collin, P., Surana, U., et al. (2014). Homeostatic control of polo-like kinase-1 engenders non-genetic heterogeneity in G2 checkpoint fidelity and timing. Nat. Commun. 5:4048.
Liu, K., Husler, J., Ye, J., Leonard, S. S., Cutler, D., Chen, F., et al. (2001). On the mechanism of Cr (VI)-induced carcinogenesis: dose dependence of uptake and cellular responses. Mol. Cell. Biochem. 222, 221–229. doi: 10.1007/978-1-4615-0793-2_26
Luczak, M. W., Green, S. E., and Zhitkovich, A. (2016). Different ATM Signaling in Response to Chromium(VI) Metabolism via ascorbate and nonascorbate reduction: implications for in vitro models and toxicogenomics. Environ. Health Perspect. 124, 61–66. doi: 10.1289/ehp.1409434
Malumbres, M., and Barbacid, M. (2009). Cell cycle, CDKs and cancer: a changing paradigm. Nat. Rev. Cancer 9, 153–166. doi: 10.1038/nrc2602
Mercer, T. R., Dinger, M. E., and Mattick, J. S. (2009). Long non-coding RNAs: insights into functions. Nat. Rev. Genet. 10, 155–159.
Nickens, K. P., Patierno, S. R., and Ceryak, S. (2010). Chromium genotoxicity: a double-edged sword. Chem. Biol. Interact. 188, 276–288. doi: 10.1016/j.cbi.2010.04.018
NOISH (2008). Occupational Exposure to Hexavalent Chromium: Department of health and human services, center for Disease Control and Prevention, National Institute for Occupational Safety and Health. Washington, D.C: NOISH.
Pollack, D., Xiao, Y., Shrivasatava, V., Levy, A., Andrusier, M., D’Armiento, J., et al. (2015). CDK14 expression is down-regulated by cigarette smoke in vivo and in vitro. Toxicol. Lett. 234, 120–130. doi: 10.1016/j.toxlet.2015.02.006
Proctor, D. M., Suh, M., Campleman, S. L., and Thompson, C. M. (2014). Assessment of the mode of action for hexavalent chromium-induced lung cancer following inhalation exposures. Toxicology 325, 160–179. doi: 10.1016/j.tox.2014.08.009
Rager, J. E., Suh, M., Chappell, G. A., Thompson, C. M., and Proctor, D. M. (2019). Review of transcriptomic responses to hexavalent chromium exposure in lung cells supports a role of epigenetic mediators in carcinogenesis. Toxicol. Lett. 305, 40–50. doi: 10.1016/j.toxlet.2019.01.011
Rayner, R. E., Makena, P., Prasad, G. L., and Cormet-Boyaka, E. J. (2019). Optimization of normal human bronchial epithelial (NHBE) cell 3D cultures for in vitro lung model studies. Sci. Rep. 9:500.
Reimand, J., Arak, T., Adler, P., Kolberg, L., Reisberg, S., Peterson, H., et al. (2016). g: profiler—a web server for functional interpretation of gene lists (2016 update). Nucleic Acids Res. 44, W83–W89.
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007
Röpke, G., Ebert, O., Märten, A., Lefterova, P., Micka, B., Buttgereit, P., et al. (2000). Increase in proliferation rate and normalization of TNF-alpha secretion by blockage of gene transfer-induced apoptosis in lymphocytes using low-dose cyclosporine A. Cancer Gene Therapy 7, 1411–1413. doi: 10.1038/sj.cgt.7700256
Seidler, A., Jähnichen, S., Hegewald, J., Fishta, A., Krug, O., Rüter, L., et al. (2013). Systematic review and quantification of respiratory cancer risk for occupational exposure to hexavalent chromium. Int. Arch. Occup. Environ. Health 86, 943–955. doi: 10.1007/s00420-012-0822-0
Shu, L., Zhao, Y., Kurt, Z., Byars, S. G., Tukiainen, T., Kettunen, J., et al. (2016). Mergeomics: multidimensional data integration to identify pathogenic perturbations to biological systems. BMC Genom. 17:874. doi: 10.1186/s12864-016-3198-9
Shumilla, J. A., and Barchowsky, A. (1999). Inhibition of protein synthesis by chromium(VI) differentially affects expression of urokinase and its receptor in human type II pneumocytes. Toxicol. Appl. Pharmacol. 158, 288–295. doi: 10.1006/taap.1999.8704
Singh, J., Carlisle, D. L., Pritchard, D. E., and Patierno, S. R. (1998). Chromium-induced genotoxicity and apoptosis: relationship to chromium carcinogenesis. Oncol. Rep. 5, 1307–1318.
Smyth, G. K. (2005). “Limma: linear models for microarray data,” in Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Statistics for Biology and Health, eds Gentleman R. Carey, W. Huber, R. A. Irizarry, and S. Dudoit (New York, NY: Springer), 397–420. doi: 10.1007/0-387-29362-0_23
Standeven, A. M., and Wetterhahn, K. E. (1989). Chromium (VI) toxicity: uptake, reduction, and DNA damage. J. Am. College Toxicol. 8, 1275–1283. doi: 10.3109/10915818909009118
Stern, A. H. J. E. R. (2010). A quantitative assessment of the carcinogenicity of hexavalent chromium by the oral route and its relevance to human exposure. Environ. Res. 110, 798–807. doi: 10.1016/j.envres.2010.08.002
Sun, M., Gadad, S. S., Kim, D.-S., and Kraus, W. L. (2015). Discovery, annotation, and functional analysis of long noncoding RNAs controlling cell-cycle gene expression and proliferation in breast cancer cells. Mol. Cell 59, 698–711. doi: 10.1016/j.molcel.2015.06.023
Takahashi, Y., Kondo, K., Hirose, T., Nakagawa, H., Tsuyuguchi, M., Hashimoto, M., et al. (2005). Microsatellite instability and protein expression of the DNA mismatch repair gene, hMLH1, of lung cancer in chromate-exposed workers. Mol. Carcinogenesis 42, 150–158. doi: 10.1002/mc.20073
Thiery, J. P. (2003). Epithelial-mesenchymal transitions in development and pathologies. Curr. Opin. Cell Biol. 15, 740–746. doi: 10.1016/j.ceb.2003.10.006
Thompson, C. M., Kirman, C. R., Proctor, D. M., Haws, L. C., Suh, M., Hays, S. M., et al. (2014). A chronic oral reference dose for hexavalent chromium-induced intestinal cancer. J. Appl. Toxicol. 34, 525–536. doi: 10.1002/jat.2907
Tu, Y., Liu, H., Zhu, X., Shen, H., Ma, X., Wang, F., et al. (2017). Ataxin-3 promotes genome integrity by stabilizing Chk1. Nucleic Acids Res. 45, 4532–4549. doi: 10.1093/nar/gkx095
Wakeman, T. P., Kim, W.-J., Callens, S., Chiu, A., Brown, K. D., and Xu, B. (2004). The ATM-SMC1 pathway is essential for activation of the chromium[VI]-induced S-phase checkpoint. Mutat. Res. 554, 241–251. doi: 10.1016/j.mrfmmm.2004.05.006
Wakeman, T. P., Wyczechowska, D., and Xu, B. (2005). Involvement of the p38 MAP kinase in Cr(VI)-induced growth arrest and apoptosis. Mol. Cell. Biochem. 279, 69–73. doi: 10.1007/s11010-005-8216-1
Wang, X., Jia, Y., Fei, C., Song, X., and Li, L. (2016). Activation/Proliferation-associated Protein 2 (Caprin-2) positively regulates CDK14/Cyclin Y-mediated Lipoprotein Receptor-related Protein 5 and 6 (LRP5/6) Constitutive Phosphorylation. J. Biol. Chem. 291, 26427–26434. doi: 10.1074/jbc.m116.744607
Wang, Z., Lin, H. P., Li, Y., Tao, H., Yang, P., Xie, J., et al. (2019). Chronic hexavalent chromium exposure induces cancer stem cell-like property and tumorigenesis by increasing c-Myc expression. Toxicol. Sci. 172, 252–264. doi: 10.1093/toxsci/kfz196
Wetterhahn, K. E., and Hamilton, J. W. (1989). Molecular basis of hexavalent chromium carcinogenicity: effect on gene expression. Sci. Total Environ. 86, 113–129. doi: 10.1016/0048-9697(89)90199-x
Xie, H., Wise, S. S., Holmes, A. L., Xu, B., Wakeman, T. P., Pelsue, S. C., et al. (2005). Carcinogenic lead chromate induces DNA double-strand breaks in human lung cells. Mutat. Res. 586, 160–172. doi: 10.1016/j.mrgentox.2005.06.002
Xin, L., Zhao, R., Lei, J., Song, J., Yu, L., Gao, R., et al. (2019). SND1 acts upstream of SLUG to regulate the epithelial-mesenchymal transition (EMT) in SKOV3 cells. FASEB J. 33, 3795–3806. doi: 10.1096/fj.201801728r
Yao, H., Guo, L., Jiang, B. H., Luo, J., and Shi, X. (2008). Oxidative stress and chromium(VI) carcinogenesis. J. Environ. Pathol. Toxicol. Oncol. 27, 77–88. doi: 10.1615/jenvironpatholtoxicoloncol.v27.i2.10
Yao, R.-W., Wang, Y., and Chen, L.-L. (2019). Cellular functions of long noncoding RNAs. Nat. Cell Biol. 21, 542–551.
Yip, A. M., and Horvath, S. (2007). Gene network interconnectedness and the generalized topological overlap measure. BMC Bioinform. 8:22. doi: 10.1186/1471-2105-8-22
Zecevic, A., Menard, H., Gurel, V., Hagan, E., DeCaro, R., and Zhitkovich, A. (2009). WRN helicase promotes repair of DNA double-strand breaks caused by aberrant mismatch repair of chromium-DNA adducts. Cell Cycle 8, 2769–2778. doi: 10.4161/cc.8.17.9410
Zhang, Q., Sun, X., Xiao, X., Zheng, J., Li, M., Yu, M., et al. (2018). Maternal chromium restriction induces insulin resistance in adult mice offspring through miRNA. Int. J. Mol. Med. 41, 1547–1559.
Zhang, Z., Leonard, S. S., Wang, S., Vallyathan, V., Castranova, V., and Shi, X. (2001). Cr (VI) induces cell growth arrest through hydrogen peroxide-mediated reactions. Mol. Cell. Biochem. 222, 77–83. doi: 10.1007/978-1-4615-0793-2_10
Keywords: hexavalent chromium, transcriptome, cell cycle, lncRNA, regulation network
Citation: Zheng P, Kang Y, Han S, Feng H, Ha F, Long C, Zhou D, Hu G, Chen Z, Wang Z, Wang T and Jia G (2021) A Novel Transcriptome Integrated Network Approach Identifies the Key Driver lncRNA Involved in Cell Cycle With Chromium (VI)-Treated BEAS-2B Cells. Front. Genet. 11:597803. doi: 10.3389/fgene.2020.597803
Received: 22 August 2020; Accepted: 26 November 2020;
Published: 13 January 2021.
Edited by:
Radu Corneliu Duca, Laboratoire National de Santé (LNS), LuxembourgReviewed by:
Tomokazu Tomo Fukuda, Iwate University, JapanHelmut Zarbl, The State University of New Jersey, United States
Copyright © 2021 Zheng, Kang, Han, Feng, Ha, Long, Zhou, Hu, Chen, Wang, Wang and Jia. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Guang Jia, jiaguangjia@bjmu.edu.cn; Zengmiao Wang, wangzengmiao@gmail.com