RNA Sequencing Analyses Reveal the Potential Mechanism of Pulmonary Injury Induced by Gallium Arsenide Particles in Human Bronchial Epithelioid Cells

Extensive use of gallium arsenide (GaAs) has led to increased exposure to humans working in the semiconductor industry. This study employed physicochemical characterization of GaAs obtained from a workplace, cytotoxicity analysis of damage induced by GaAs in 16HBE cells, RNA-seq and related bioinformatic analysis, qRT-PCR verification and survival analysis to comprehensively understand the potential mechanism leading to lung toxicity induced by GaAs. We found that GaAs-induced abnormal gene expression was mainly related to the cellular response to chemical stimuli, the regulation of signalling, cell differentiation and the cell cycle, which are involved in transcriptional misregulation in cancer, the MAPK signalling pathway, the TGF-β signalling pathway and pulmonary disease-related pathways. Ten upregulated genes (FOS, JUN, HSP90AA1, CDKN1A, ESR1, MYC, RAC1, CTNNB1, MAPK8 and FOXO1) and 7 downregulated genes (TP53, AKT1, NFKB1, SMAD3, CDK1, E2F1 and PLK1) related to GaAs-induced pulmonary toxicity were identified. High expression of HSP90AA1, RAC1 and CDKN1A was significantly associated with a lower rate of overall survival in lung cancers. The results of this study indicate that GaAs-associated toxicities affected the misregulation of oncogenes and tumour suppressing genes, activation of the TGF-β/MAPK pathway, and regulation of cell differentiation and the cell cycle. These results help to elucidate the molecular mechanism underlying GaAs-induced pulmonary injury.

GO analysis of differential gene expression induced by GaAs. First, GO and KEGG analyses were performed on 7285 significantly dysregulated mRNAs in GaAs vs NC. We derived 85 highly enriched GO terms (adj. P val <0.05) and 27 significantly enriched pathways (adj. P val <0.05). The top 10 significant GO terms during biological process (BP), cellular component (CC) and molecular function (MF) induced by GaAs are displayed in Fig. 3a-c. Several GO terms, such as cellular response to chemical stimulus, regulation of signalling, cell differentiation, cell cycle, focal adhesion, transcription regulator activity and small molecule binding, were closely related to GaAs exposure. As shown in Fig. 3d, the significantly upregulated pathways induced by GaAs included transcriptional misregulation in cancer, cytokine-cytokine receptor interaction, MAPK signalling pathway, TGF-β signalling pathway and pathways in cancer. Significantly downregulated pathways, such as metabolic pathways, small cell lung cancer, cell cycle, and fatty acid metabolism, were detected (Fig. 3e). The summaries of genes involved in the significant up-and downregulation pathways are listed in Supplementary Table S2. DO analysis of differential gene expression induced by GaAs. DO (Disease Ontology) is a database describing human gene function and disease that can be used to consider the interactions and functions of differentially expressed genes with disease. The top 20 enrichment terms related to human disease were recognized based on the DO database. The P < 0.05 gene list was also analysed for DO annotation, enabling us to study gene-disease relationships. Significantly upregulated genes induced by GaAs in 16HBE using DO analysis may be related to lung disease, respiratory system disease, lower respiratory tract disease, pulmonary fibrosis, interstitial lung disease and chronic obstructive pulmonary disease (Fig. 4a). For example, hereditary breast ovarian cancer, retinal cancer and ocular cancer were relevant to the downregulation of DO-related diseases (Fig. 4b). Taken together, these results indicate that dysregulated DO of differentially expressed genes induced by GaAs in 16HBE may be related to lung carcinoma, non-small-cell lung carcinoma and thoracic cancer (Fig. 4c).
Protein-protein interaction (PPI) analysis revealed the key genes triggered by GaAs. PPI analysis of the STRING interactome was used to screen the key genes involved in GaAs-induced toxicity in 16HBE cells. According to the significant genes from the significant pathway and DO analysis, we next focused on the DEGs in the cell cycle pathway, pathways in cancer, transcriptional misregulation in cancer, small-cell lung cancer, the MAPK signalling pathway, the TGF-β signalling pathway and non-small-cell lung cancer, and the corresponding significant genes were mapped to the corresponding molecular interaction database (Fig. 5a). Seventeen key genes with the top degree (>20) involved in the above 6 pathway interaction network were screened (Table 2), which may play an important role in GaAs-induced pulmonary toxicity in 16HBE cells. These 17 key genes consisted of 10 upregulated genes (FOS, JUN, HSP90AA1, CDKN1A, ESR1, MYC, RAC1, CTNNB1, MAPK8 and FOXO1) and 7 downregulated genes (TP53, AKT1, NFKB1, SMAD3, CDK1, E2F1 and PLK1) (Fig. 5b).

qRT-PCR verification and survival curve analysis.
To verify the RNA-seq analysis results, the expression of the above 17 significantly dysregulated genes selected from PPIs was verified using qRT-PCR (Fig. 6a). These selected genes were related to various pathways (Table 2), such as small-cell lung cancer (AKT1, TP53, NFKB1 and E2F1), MAPK signalling pathway (JUN, RAC1, MAPK8, MYC and FOS) and transcriptional misregulation in cancer (CDKN1A, MYC and FOXO1). The potential correlation between the expression of the above 17 significantly dysregulated genes from RNA-seq and qRT-PCR verification was analysed with the Spearman rank www.nature.com/scientificreports www.nature.com/scientificreports/ order correlation test. We found that the gene expression results of qPCR were positively correlated with the corresponding results of RNA-seq ( Fig. 6b; r = 0.887, P = 0.000). The results suggest that good consistency was obtained between RNA-seq analysis and qRT-PCR verification.
Then, we performed 10 upregulated genes for survival curve analysis using TCGA data by UALCAN. Further survival analyses on these key genes were employed to evaluate their effects on the survival of lung cancer. Figure 6c-e shows that HSP90AA1, RAC1 and CDKN1A expression levels were clearly related to the prognosis of lung adenocarcinoma (LUAD) or lung squamous cell carcinoma (LUSC) patients, which indicated that high expression of the three genes was significantly associated with a lower rate of overall survival.  www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Workers in the semiconductor manufacturing industry are potentially exposed to inhalable GaAs dusts. Because GaAs has a lower solubility than any other arsenic compound and is associated with lung toxicity, even linked to lung neoplasm risk, exploring the toxicogenomic effects of GaAs on human bronchial epithelial cells is critical to understanding the mechanisms of these adverse health effects 7,16 . In this study, we first examined the physicochemical characterization of GaAs from a workplace and cytotoxicity induced by GaAs and then conducted RNA-seq, related bioinformatic analysis, qRT-PCR verification and survival analysis to comprehensively understand potential genes and pathways that lead to lung toxicity induced by GaAs. Our aim was to elucidate the role of GaAs in 16HBE cells and to provide a mechanistic explanation for GaAs exposure increasing pulmonary injury risk in humans. This research was the first attempt to perform transcriptome sequencing analysis to fully understand GaAs-induced toxicity in 16HBE cells.
GaAs could release gallium and arsenic moieties in vitro though a relatively lower dissolution rate 17 . In our study, the dissolution rate of GaAs particles was approximately 2.3% after 24 h of exposure to 16HBE cells in the cell supernatant. In addition to dissolution of the above two metals, the toxic effect depends on particle size, exposure duration and exposure route 18 . We collected the GaAs from a workplace that reflects the real environmental exposure, in which the hydrodynamic median sizes of GaAs particles are 2.51 μm or 2.54 μm when suspended in dH 2 O or MEM, respectively. It was demonstrated that micron-sized GaAs particles can result in hazardous pulmonary effects in animal studies, which indicates that a smaller fraction of GaAs is a relatively more severe pneumotoxicant 19,20 . It was reported that a strong inflammatory response induced by solid GaAs in the lungs disappeared when the gallium and arsenic oxides dissolved, as shown in animal models, in which the unchanged GaAs particles were cleared from the lung 20 . Then, the toxic effects of micron-sized GaAs particles exposed to 16HBE were evaluated.
The lung is one of the target organs for the toxic effects of GaAs, as this organ has higher absorption after intratracheal administration than is observed after oral dosing 5 . In inhalation studies, lung retention of inhaled GaAs dust has been shown to be influenced by toxic effects from GaAs itself 19 . GaAs is a compound of gallium and arsenic that is responsible for the toxicity properties of these two metals. Gallium is considered mildly toxic because most of this metal is excreted through urine or faeces, which does not contribute significantly to the lung toxicity of GaAs 21 . Arsenic exposure has been shown to induce human tumourigenesis, and the lung is one of the main targets 22 . Toxicities of arsenic-induced carcinogenicity, including generation of oxidative stress, altered cell proliferation, changes in DNA methylation and co-carcinogenesis, though the exact molecular mechanism governing this toxicity is not well understood 18,23,24 . The toxicity of GaAs in pulmonary tissue includes inflammation, lung weight increase, fibrosis, seropurulent pneumonia and pneumocyte hyperplasia 18,20,25 .
Our results indicated that GaAs induced abnormal gene expression mainly related to cellular response to chemical stimulus, regulation of signalling, cell differentiation, the cell cycle and regulation of signalling. GaAs exposure in 16HBE cells induced the upregulation of cell differentiation-related genes and the downregulation of cell cycle-related genes, which indicated that GaAs is an inducer of cell differentiation. In 16HBE cells after 24 h of GaAs exposure, we also observed morphological changes, such as irregular shapes, smaller volumes and abnormal nuclei. GaAs-related dysregulated genes are involved in signal transduction pathways, such as transcriptional www.nature.com/scientificreports www.nature.com/scientificreports/ misregulation in cancer, cytokine-cytokine receptor interactions, the MAPK signalling pathway and the TGF-β signalling pathway. It has been summarized that arsenic-induced lung tumours may, through disruption of the PI3K/AKT signalling pathway, activate the EGFR signalling pathway and affect the NRF2 signalling pathway to play a role in carcinogenesis 22 . The TGF-β and MAPK pathways play critical roles in cell development and cell cycle regulation, even in tumour formation and metastasis. It has been shown that activation of MAPK and TGF-β could be induced by ROS accumulation, and the MAPK pathway and TGF-β pathways were closely related to lung fibrosis 26,27 . TGF-β can activate MAPK signalling, and inhibition of the TGF-β/MAPK pathway could protect against lung fibrosis 28,29 , which may explain the GaAs toxicity in pulmonary tissue. Then, a view of the global gene-disease relationships of GaAs exposure in 16HBE was obtained by DO analysis. Lung disease and respiratory system disease were found in significantly upregulated genes induced by GaAs in 16HBE. All dysregulated genes in the significant DO term were clustered into tumour-related diseases, including lung carcinoma and non-small-cell lung carcinoma. These results indicated that pulmonary disease-related pathways were affected when GaAs exposure occurred in 16HBE.
Then, we screened the 17 key genes with the top degree (>20) involved in significantly related pathways of GaAs-induced pulmonary toxicity in 16HBE cells, whose gene expression levels were verified by qPCR. Among these key genes, several oncogenes (FOS, JUN, MYC) and tumour suppressor genes (TP53) were associated with lung cancer development. ESR1/2 might directly or indirectly regulate oncogenic pathways in non-small-cell lung cancer, and FOXO1 expression is a favourable prognostic factor in non-small-cell lung cancer 30,31 . CTNNB1 has been found to be genetically mutated in various human cancers, including lung adenocarcinoma, whose gene expression was increased in lung tissue with pulmonary fibrosis 32 . MAPK8 was found to be potentially related to lung cancer 33 . Finally, the survival curve analysis of 10 upregulated genes using TCGA data indicated that high  www.nature.com/scientificreports www.nature.com/scientificreports/ expression of HSP90AA1, RAC1 and CDKN1A induced by GaAs was clearly related to the lower survival probability in LUAD or LUSC patients. HSP90AA1 was found to be directly associated with lung cancer 34 . The overexpression of Rac1 was linked to aggressive growth and other malignant characteristics of tumours, and a high level of Rac1 could predict a poor prognosis in different types of cancer 35 . CDKN1A expression is upregulated in long-term oxidative stress-induced experimental bronchopulmonary dysplasia and in a hyperoxia-induced lung injury rat model 36,37 . The results of this study indicate that GaAs-associated toxicities affect the misregulation of oncogenes and tumour suppressing genes, activation of the TGF-β/MAPK pathway, and regulation of cell differentiation and the cell cycle.
Methods particle sample collection. GaAs particles, obtained from a manufacturer of gallium arsenide crystals in Beijing, were collected using an IOM personal sampler (SKC USA) equipped with a cylindrical body, 37 mm cassette, and PVC filter with enough time. Then, the samples were rinsed with distilled water, and the weight change was recorded to calculate the collection quality of the particles. We chose the wafer manufacturing process as the sampling site in which workers are only exposed to GaAs particles. cell culture. The human bronchial epithelial cell line (16HBE14O-, abbreviated as 16HBE) was a gift from Dr D.C. Gruenert (University of California, San Francisco, USA). The cells were maintained in MEM culture medium (Gibco, USA) supplemented with 10% foetal bovine serum (Gibco, USA), 100 U/mL penicillin and 100 μg/mL streptomycin and cultured at 37 °C in a 5% CO 2 humidified environment.

Assessment of cell viability.
To evaluate the cytotoxicity and cell viability induced by GaAs particles, the Cell Counting Kit-8 (CCK-8) assay and Cytotoxicity LDH Assay Kit-WST (LDH) 38 assay were employed to assess mitochondrial dehydrogenase activity and loss of cell integrity, respectively. Three biological replicates were performed. Cytotoxicity EC 20 (effective concentration resulting in 20%) of 16HBE cells exposed to GaAs particles was determined. The cells were exposed to GaAs particles at 5.  physicochemical characterization of GaAs particles. Scanning electron microscopy (SEM HITACHI S4800) was used to examine the particle morphology. Characterization of GaAs particles was performed at a concentration identical to the cytotoxicity EC 20 . The hydrodynamic diameter of GaAs particles was determined by dynamic light scattering (DLS) using a MASTERSIZER 2000 (Malvern Instruments, Malvern, UK). GaAs particles were suspended in dH 2 O and in MEM at the concentration of cytotoxicity EC 20 and sonicated prior to measurement.
Dissolution of GaAs particles was assessed in full MEM. GaAs particles were incubated at the concentration of Cytotoxicity EC 20 for 24 h at 37 °C. Three biological replicates were performed. After incubation, the samples were centrifuged at 25,000 × g for 30 min, and the supernatant was used for Ga and As determination by inductively coupled plasma mass spectrometry (ICP-MS, Thermo).

GaAs particle exposure.
For experiments, the cells were seeded in culture plates at a density of 1 × 10 5 cells/mL, allowed to attach for 24 h, and treated with GaAs particles suspended in MEM culture medium of certain concentrations for another 24 h. A suspension of GaAs particles was dispersed by a sonicator (Bioruptor UDC-200, Belgium), diluted to EC 20 concentrations, and immediately added to 16HBE cells. In this study, 6-well culture plates were used, the bottom area of which was 9.6 cm 2 . To ensure that each unit area of cells was exposed to the same amount of GaAs particles, the exposure volume of the GaAs particle suspension was calculated according to the bottom area of the cell culture plate. Cells maintained in MEM culture medium without GaAs particles were used as the control group. Three biological replicates were performed.
Total RNA extraction and RNA-seq. Total RNA was extracted from the 16HBE cells exposed to EC20 of GaAs particles and negative control (NC) using TRIzol reagent (Invitrogen, San Diego, CA, USA) according to the manufacturer's instructions. Both groups were conducted in triplicate. A total amount of 3 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using the NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, USA), and index codes were added to attribute sequences to each sample. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina). After cluster generation, the library preparations were sequenced on an Illumina HiSeq platform to generate 150 bp paired-end reads (Novogene, Beijing).
Bioinformatic analysis. The raw sequence files generated from 6 files (bam) in this study has been deposited to NCBI's Sequence Read Archive (SRA) database with the accession number PRJNA623863. Clean data (clean reads) were obtained by removing reads containing adapters, reads containing poly-N and low-quality reads from raw data using in-house Perl scripts (Novogene, Beijing). FPKM (fragments per kilobase of exon per million fragments) of each gene was calculated based on the length of the gene and read count mapped to this gene. Differential expression analysis of two groups (three biological replicates per group) was performed using the DESeq. 2 R package (1.16.1). The resulting P-values were adjusted using Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value <0.05 found by DESeq. 2 were assigned as differentially expressed.
Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of differentially expressed genes was implemented using the clusterProfiler R package and iDEP program 39 , in which gene length bias was corrected. Identified GO terms and KEGG pathways with a corrected p < 0.05 were considered significantly enriched 40 . Disease ontology (DO) annotates human gene-disease relationships, which is important annotation in translating molecular findings from high-throughput data to clinical relevance 41 . DO disease terms and semantic associations were obtained through the DOSE R package, which provided DO terms with a corrected p < 0.05 were considered significantly enriched. Protein-protein interaction (PPI) analysis of differentially expressed genes was based on the STRING database, and the images were examined using the NetworkAnalyst platform 42 .
Quantitative real-time PCR (qRT-PCR) analysis. The RNeasy Mini Kit (Qiagen, Hilden, Germany) was used to isolate total RNA from 16HBE cells exposed to EC20 of GaAs particles and NC, as described above for the GaAs particle exposure method. Reverse transcription to synthesize first-strand cDNA was conducted using the SuperScript II First-stand Synthesis System for RT-PCR (Invitrogen, USA). The gene expression levels were assessed with the use of TaqMan qRT-PCR assays to validate the GaAs particle-related key genes acquired from RNA-seq. All TaqMan qRT-PCR reactions were carried out using TaqMan ® Fast Advanced Master Mix (Applied Biosystems, USA) on a ViiA 7 Real-Time PCR system (Applied Biosystems, USA) in 384-well plates. The relative mRNA expression levels of the target genes were normalized to GAPDH (housekeeping gene), and the fold change was calculated by using the 2 −ΔΔCt method. The primers and reaction conditions are listed in Supplementary Table S1.
Survival analysis. The survival analysis of key GaAs particle exposure-related genes was performed on UALCAN (http://ualcan.path.uab.edu/) to conduct the Kaplan-Meier estimator 43 . The survival curves of samples with high gene expression and low/medium gene expression were compared by the log rank test. P values <0.05 were considered to be significant. Statistical analysis. Statistical analysis was performed using SPSS 17.0 software (SPSS, Chicago, IL, USA).
The Mann-Whitney test was used to compare the qPCR results between groups. The potential correlation between