PTEN Hypomethylation Could Be A Biomarker For Early Detection of Silicosis

Weijiang Hu Chinese Center for Disease Control and Prevention, Beijing 100050, China Jingbo Zhang Tongji University School of Medicine Kai Liu Chinese Center for Disease Control and Prevention, Beijing 100050, China Jie Liu Department of Occupational Disease Yuxin Zheng Qingdao University Xin Sun Chinese Center for Disease Control and Prevention, Beijing 100050, China Liangying Mei Hubei Provincial Center for Disease Control and Prevention Zushu Qian Huangshi Center for Disease Control and Prevention Qiangguo Sun Huangshi Center for Disease Control and Prevention Qiang Liu Suzhou Center for Disease Control and Prevention Zhijun Wu Chinese Center for Disease Control and Prevention Hengdong Zhang Jiangsu Provincial Center for Disease Control and Prevention Yanping Li Honghe Prefecture Third People’s Hospital Daoyuan Sun (  dysun@163.com ) Tongji University School of Medicine Meng Ye Chinese Center for Disease Control and Prevention


Introduction
Globally, 39% of total pneumoconiosis cases were ascribed to silicosis in 2017 [1]. The silicosis is mainly caused by excessive inhalation of silica dust which, over time, leads to lung in ammation and brosis [2].
Besides, dysregulated mRNAs and miRNAs are implicated in immune response and in ammation [3], and differences in mRNAs and miRNAs expressions between silicosis and normal subjects have been reported in many studies [4][5][6][7][8][9][10]. However, most subjects of current silicosis related transcriptome studies concentrate mainly on human blood [6,7], bronchoalveolar lavage uid [4], macrophages [9], lung epithelial cells [8], or rats model [10][11][12] instead of silica dust major target organ-human lung tissue. In addition, despite the same silicosis rats model, the numbers of differential mRNA or miRNA expression were still controversial [5,10,13,14]. Hence, it is critical and essential to conduct a transcriptome study based on silicosis patients and normal people lung tissues. Our previous study revealed that PTEN promoter hypermethylation might be associated with decrease of PTEN protein in silicosis lung tissues [15], whereas its expression pro le in silicosis patients blood samples was still unclear.
Herein, we reported a transcriptome study based on fteen silicosis patients and eight normal people lung tissues, and made further quantitative real-time PCR (RT-qPCR) validation in other ten lung tissues and highlights the great signi cance of exploring silicosis underlying toxicological mechanisms in lung tissues.

Methods
Lung tissues and blood samples of silicosis patients and normal people Silicosis patients and normal people lung tissues were obtained from National Institute for Occupational Health and Poison Control, Chinese Center for Disease Control and Prevention. We selected patients with silicosis who had undergone autopsying between 1967 and 1979, whereas cases combined with lung cancer or other pulmonary diseases were excluded. In addition, cases were divided into early stage group and advanced stage group according to diagnosis of occupational silicosis: early stage group included stage ; advanced stage group included stage and [15,34]. A total ofsix early stage silicosis cases, nine advanced stage silicosis cases, and eight normal lung tissues were included for microarray processing, and RT-qPCR validation.
Besides, blood samples of 414 male silicosis patients and 177 male normal people were collected for further RT-qPCR validation during 2015. The silicosis patients, diagnosed according to GBZ 70-2009, were selected from two silica dust-exposed enterprises: 353 of one copper mine, and 61 of one stone processing factory. 177 normal people were selected as control group from one local food processing company. The cases inclusion criteria were: 1) occupational silica dust exposure; 2) diagnosed silicosis according to GBZ 70-2009; and 3) no other occupational dust exposure except for silica dust. The control samples inclusion criteria were: 1) no occupational dust exposure; 2) no cancer; 3) no pulmonary brosis under imaging examination; and 4) no hereditary lung diseases. Finally, a total of 404 silicosis patients and 177 normal people blood samples were included.
All the lung tissues and blood samples were obtained through the medically prescribed test, including informing each participant, and written informed consent was obtained from all subjects. Meanwhile, this study was approved by the Ethics Committee Board in Occupational Health and Poison Control, Chinese Center for Disease Control and Prevention. Within ten minutes after sampling, lung tissue and blood samples were stored at -80℃ refrigerator. The main clinical characteristics of silicosis patients and normal people were presented in Supplementary Table 21 and Supplementary Table22: Supplementary  Table 21 of lung tissues, and Supplementary Table 22 of blood samples. Furthermore, there were no statistically signi cant differences in the age, gender, and smoking between silicosis patients and normal people in Supplementary Table 21 or Supplementary Table 22. Microarray Processing And Analyses RNA preparation Three early stage silicosis, ve advanced silicosis, and four normal lung tissues were randomly selected for microarray processing and analyse.The total RNA in lung tissues was extracted using Trizol Reagent (Invitrogen, Carlsbad, Canada), and puri ed with an RNeasy Mini Kit (Qiagen, Hilden, Germany) according to the manufacture's protocol, including a DNase digestion treatment. The amount and quality of RNA were determined by a UV-Vis Spectrophotometer (Thermo, NanoDrop 2000, USA) at the absorbance of 260 nm. miRNA microarray processiong The miRNA microarray pro ling was performed using GeneChip® miRNA 4.0 Array (Santa Clara, CA, USA) according to manufacturer's recommended protocol. Brie y, 1 µg of total RNA from samples was labeled by polyA polymerase using the Genisphere FlashTag HSR kit according to the manufacture's recommendations (Genisphere, Hat eld, PA). In addition, RNA was hybridized to the Affymetrix miRNA array as recommended by the vendor. Standard Affymetrix array cassette staining, washing and scanning was performed using the post-hybridization kit (#900720; Affymetrix) and GeneChip Scanner 3000. Feature extraction was performed using Affymetrix Command Console software. Furthermore, the raw data were treated by the following work ow: background detection, RMA global background correlation, quantile normalization, median polish and log2-transformed with miRNA QC tool software (Affymetrix). mRNA microarray processiong The mRNA expression pro ling was measured using GeneChip Human Transcriptome Array 2.0 (Affymetrix GeneChip, Santa Clara, CA, USA), which contained 59, 302 gene-level probe sets.The microarray analysis was performed by Affymetrix Expression Console Software (Version 1.2.1). Raw data (CEL les) were normalized at transcript level using robust multi-array average method (RMA work ow). Median summarization of transcript expression was calculated. Mrna expression data were then ltered to include only those sets that were in the 'core' metaprobe list, which represented RefSeq genes.
Differentially expressed RNAs (DERs)analyses and heatmap generating Given the random-variance model F-test can raise degrees of freedom effectively in the cases of small samples, it was applied to lter the differentially expressed RNA among early stage, advanced stage and control groups. The signi cance criteria of DERs was P-value < 0.05, and heatmaps of DERs were generated using the online tool Morpheus. Furthermore, fold change (FC) was calculated as geom mean of probeset intensities in case group divided by that of control group, and differential expressions of mRNA and lncRNA amongnormal, early stage, and advanced stagelung tissues were identi ed by fold change (FC) > 2.00 or < 0.50: FC > 2.00 indicated up-regulate, and FC < 0.50 indicated down-regulate.

Gene Ontology (Go) And Pathway Analyses
The above differentially expressed mRNAs were further used to conduct Gene ontology (GO) and pathway analyses. GO analysis was applied to analyze the main function of the differential expression mRNAs according to the Gene Ontology, which could organize mRNAs into hierarchical categories and uncover the mRNAs regulatory network on the basis of biological process and molecular function.
Pathway analysis was used to nd out the signi cant pathway of the differential genes according to Kyoto Encyclopedia of Genes and Genomes (KEGG). mRNA-mRNA-network mRNA-mRNA interaction network was constructed based on the data of differentially expressed mRNAs to identify key mRNAs (betweenness centrality > 0.05 and degree ≥ 6). In network, the degree measured how correlated a gene was with all other network mRNAs. For a mRNA in the network, the number of source mRNAs of a mRNA was called the indegree of the mRNA, and the number of target mRNAs of a mRNA was its outdegree. The character of mRNAs was described by betweenness centrality measures re ecting the importance of a node in a graph relative to other nodes.

Series test of cluster (STC) and STC-GO analyses
STC analysis was performed using normal, early stage and advanced stage lung tissues to explore possible changes in differentially expressed mRNAs and miRNAs expression patterns during the process of silicosis. Therefore, we set 16 distinct and representative signi cant model pro les for the above mRNAs and miRNAs, respectively. Given that there might be a large amount of signi cant trends, we selected mRNAs in patterns with opposite trends to those of miRNAs expression patterns to conduct further analyses. Furthermore, mRNAs and miRNAs in patterns with the completely opposite trends and signi cance were selected to conduct network analysis; mRNAs and miRNAs in patterns with opposite trends and more signi cant differences were selected to make RT-qPCR validation.

miRNAs target mRNAsand the intersecting mRNAs
The miRNAs in pattern with a completely opposite trend to mRNAs expression pattern were used to predict their target mRNAs by miRanda, and the intersecting mRNAs were extracted between the above target mRNAs and those mRNAs in pattern with an opposite trend to miRNAs expression pattern. miRNA-mRNA-network miRNA-mRNA interaction network was constructed based on the data of intersecting mRNAs and miRNAs, which were in pattern with an opposite trend to mRNAs expression pattern, according to their interactions in the Sanger miRNA database. Besides, the center of this network was represented by degree, which was the contribution of one miRNA to the mRNAs around or the contribution of one mRNA to the miRNAs around, and the key miRNA and mRNA in the network always had the biggest degrees.
RT-qPCR validation RNA preparation The blood samples and remaining lung tissues, including three early stage silicosis, four advanced stage silicosis, and three normal lung tissues, were used to make RT-qPCR. Moreover, the above seven silicosis lung tissues were considered as case group, and other three normal lung tissues as control group. Quantitative real-time PCR was performed using SYBR Green assays (ABI,7500Fast), including primers synthesized by CNKINGBIO (Beijing, China) (Supplementary Table 23 and Supplementary Table  24), and β-actin was used as an internal control.
RT-qPCR validation in lung tissues According to mRNA-mRNA-network, STC, and mRNA-miRNA-network analyses results, a representative set of ten mRNAs (Supplementary Table 23) and nine miRNAs (Supplementary Table24) were selected to make RT-qPCR validation. The above ten mRNAs were in pattern with more signi cant difference and opposite trend to that of miRNA expression pro le. In addition, among the above ten mRNAs and nine miRNAs, nine mRNAs had high degrees (degree > 6) in mRNA-mRNA-network, and PTEN was selected given our previous study 15 ; and the nine miRNAs were under selected mRNAs directed interactions.
RT-qPCR validation in blood samples Based on the above RT-qPCR validation results in lung tissues and previous studies [15,[35][36][37][38], PTEN and GNAI3, as well as hsa-miR-8063 and hsa-miR-181b-5p which regulated PTEN (Supplementary Table9to 12) were selected as markers for further validation in an independent cohort of 404 silicosis patients and 177 normal people blood samples.

Bisul te Sequencing PCR (BSP)
Accoding to the above RT-qPCR results, PTEN, an important tumor suppressor gene, had completely opposite expression levels between lung tissues and blood samples. So bisul te sequencing PCR was used to detect PTEN methylation state change in case group and control group blood samples. PTEN was extracted using the Genomic DNA Puri cation Kit (Promega, Madison, WI, USA) based on randomly selected control and case groups blood samples. In addition, 1 µg PTEN was treated with bisul te using EpiTect Fast Bisul te Conversion Kit (Cat No. 59824, Qiagen, Hilden, Germany) as follows: bisul te conversion of PTEN, puri cation of converted PTEN, primer design, and PCR ampli cation of bisul teconverted PTEN. Furthermore, the BSP speci c primers were designed according to the location of the target genes CpG islands (Supplementary Table 25 and Supplementary Text).

Statistical Analyses
Random-variance model F-test was applied to lter the differentially expressed RNAs, and differential expression RNAs were selected at a logical sequence according to RVM (Random variance model) corrective ANOVA. Besides, the signi cant GO and pathway categories were examined by Fisher's exact test and χ 2 test. The criteria of signi cance was de ned by P-value < 0.05. The above statistical analyses were performed using SPSS 23.0 (SPSS Inc., Chicago, IL, USA).

Results
Differentially expressed mRNAs and miRNAs in silicosis patients and normal people lung tissues According to quality-control tests results, there was one early stage silicosis lung tissue excluded due to no detection in mRNA microarray processing. Therefore, a total of eleven and twelve lung tissues included in mRNA and miRNA microarray analyses, respectively: seven and four of case and control groups, respectively, in mRNA microarray analysis; eight and four of case and control groups, respectively, in miRNA microarray analysis. As a result, we totally identi ed 1417 differentially expressed mRNAs ( Fig.  1a

GO and pathway analyses
On the basis of above results, we conducted gene ontology (Go) and pathway analyses to analyze main functions and nd out signi cant pathways of differential expression mRNAs. The GO results indicated that differentially expressed mRNAs were closely related to gene expression, viral reproduction, viral infectious cycle functions, and etc ( Supplementary Fig. 1a,Supplementary Table 3 and 4). In addition, these differential expression mRNAs were enriched in phagosome, ribosome, olfactory transduction, antigen processing and presentation, PI3K-Akt pathways, and etc ( Supplementary Fig. 1b, Supplementary  Table 5). Notably, an unexpected discovery was silica-exposure induced differentially expressed mRNAs involving in proteoglycans in cancer and pathways in cancer, and mRNAs included PIK3R3, HIF1A, FN1,  Table 6).

mRNA-mRNA-network
To identify key mRNAs in differentially expressed mRNAs, mRNA-mRNA functional interaction network was conducted, and the result revealed nine key mRNAs: SOCS3, KIT, STAT3, FIGF, CTNNA1, PIK3R3, HIF1A, NFKB1, and TJP1 ( Supplementary Fig. 2,Supplementary Table 7 and 8). Among the above nine mRNAs, SOCS3 had the highest centrality, and PIK3R3 had the highest degree value. SOCS3 expression was promoted by NFKB1, whereas inhibited by STAT3. Besides, FIGF could activate KIT, and KIT activated STAT3 and PIK3R3. Series test of cluster (STC) and STC-GO analyses STC analyses results showed ve and three patterns of mRNAs and miRNAs expressopm pro les, repectively,with signi cant trends (P< 0.05) (Supplementary Table 9 to 12). Among these patterns, mRNA expression pro le 5 had completely opposite trend to that of miRNA expression pro le 8: mRNAs expression in mRNA expression pro le 5 (Fig. 2a) increased in early stage, and reduced in advanced stage; whereas miRNAs expression in miRNA expression pro le 8 (Fig. 2b) reduced in early stage, and increased in advanced stage. STC-GO analysis suggested mRNAs in expression pro le 5 (Fig. 2a) were mainly involved in detection of chemical stimulus involved in sensory perception of smell, G-protein coupled receptor signaling pathway, and sensory perception of smell, etc. Besides, mRNA expression pro le 3 and miRNA expression pro le 13, mRNA expression pro le 13 and miRNA expression pro le 3 had opposite trends, respectively, with more signi cant differences.

miRNAs target mRNAs and the intersecting mRNAs
Based on STC analyses, miRNAs in pattern with completely opposite trend to mRNAs expression pattern were used to predict their target mRNAs by miRanda. Finally, the miRNAs in miRNA expression pro le 8 were used to make further prediction, and we found 12781 target mRNAs (Supplementary Table 13). Moreover, there were 39 intersecting mRNAs extracted between the above 12781 target mRNAs and mRNAs in mRNA expression pro le 5 ( Fig. 2a; Supplementary Table 9 and Supplementary Table 10).

miRNA-mRNA-network
The above intersecting mRNAs and miRNAs in miRNA expression pro le 8 were used to built miRNA-mRNA interaction network. Results revealed three and eight key mRNAs and miRNAs, respectively, in this network ( Fig. 3; Supplementary Table 14 and Supplementary Table 15). Both SLC7A6OS and FAM126B had the highest degree value among mRNAs, and hsa-miR-31-5p had the highest degree among miRNAs. Besides, other two miRNA-mRNA functional interaction networks were built based on mRNA expression pro le 3 and miRNA expression pro le 13, mRNA expression pro le 13 and miRNA expression pro le 3, respectively, and results suggested UBN2, QKI, hsa-miR-27a-3p, hsa-miR-27b-3p, and hsa-miR-3613-3p also played key roles (Supplementary Table 16 to 19; Supplementary Fig. 3 and Supplementary Fig. 4).

RT-qPCR validation in blood samples
Similarly, there were 16 blood samples of case group were excluded due to no detection: 11 of mRNAs RT-qPCR, and 5 of miRNAs RT-qPCR; and other 10 blood samples of control groupwere excluded in mRNAs RT-qPCR. As a result, a total of 560 samples were included in mRNAs validation: 393 of case group, and 167 of control group; and 576 samples were included in miRNAs validation: 399 of case group, and 177 of control group. The results indicated that PTEN and GNAI3 expressions were signi cantly up-regulated (P< 0.001) (Fig. 5a), whereas hsa-miR-8063 and hsa-miR-181b-5p were signi cantly down-regulated (P< 0.001) (Fig. 5b).

Bisul te Sequencing PCR (BSP)
On the basis of above results, we randomly selected 54 and 21 blood samples of case and control groups, respectively, to make further bisul te sequencing PCR (BSP) to detect PTEN methylation rate. The results indicated that PTEN in case group had signi cantly decreased methylation rate when compared with that of control group (control group: 3.54% ± 0.06, case group: 3.50% ± 0.05, P< 0.01) (Supplementary Table 20).

Discussion
Our study presents a comprehensive analysis of differential mRNA and miRNA expression between silicosis patients and normal people. To our knowledge, this is the rst study to report mRNA and miRNA expression on transcriptomic level in lung tissue. The microarray processing and analyses identi ed aberrant mRNA and miRNA expressions in global scale, and the number of mRNAs with signi cant expression differences was larger than that of miRNAs. In addition, gene expression function and phagosome pathway were signi cantly enriched in differential mRNAs expressions, and network analyses revealed nine key mRNAs in mRNA-mRNA-network. Interestingly, PTEN was identi ed as potential biomarker for silicosis.
Overdose and long-term exposure to silica dust could induce severe lung cells toxicities, including in ammation response [16], extracellular oxidative stress [17], cell apoptosis [18,19], DNA damage [18,19], brosis [20][21][22], etc. Dysregulated mRNAs and miRNAs associated with silica-induced pulmonary brosis, pneumoconiosis, or silicosis were reported in previous studies [14,[23][24][25][26][27], but the exact identify of differential mRNA and miRNA expression remains to be completely elucidated. Based on silicosis rats model, 39 altered miRNAs were identi ed involving in lung brosis [10]; based on silicosis patients bronchoalveolar lavage uid (BALF) samples, there were 110 dysregulated miRNAs, and different stages of silicosis were associated with distinct changes in miRNAs expression [4]; whereas our results showed 1417 and 241 differentially expressed mRNAs and miRNAs, respectively, and most mRNAs or miRNAs expressions had no signi cant difference between early stage and advanced stage silicosis lung tissues. The above differences might possibly be explained by miRNA high conservation, species diversity, and sample size expansion [10]. The only relevant report based on ve lung tissues of silicosis patients were conducted to validate decreased miR-486-5p expression [14], and that is consistent with our result. Furthermore, we identi ed other differentially expressed mRNAs and miRNAs which have not been reported previously.
Our enrichment analyses revealed that differential mRNAs expressions were involved in many functions and signal pathways which were rarely reported in previous studies. Most previous studies showed that in ammation and immune response were mainly involved in silica-exposure induced aberrantly expressed mRNAs and miRNAs [4-6, 8, 9, 13, 28], and these two functions were also observed in our study. Moreover, epithelial-mesenchymal transition (EMT) [23], cell autophagy [27], and alveolar structural damage [26] functions, as well as PI3K/Akt/mTOR/Snail [23], TGF-β1/Smad3 [24], and MAPK/ERK [27] signal pathways were associated with the dysregulated mRNAs and miRNAs. For instance, miR-503 was reported to modulate epithelial-mesenchymal transition in silica-induced pulmonary brosis by targeting PI3K/Akt/mTOR/Snail pathway [23]. In addition to those, functions like gene expression, viral reproduction, viral infectious cycle, translational initiation, and mRNA metabolic process, as well as signal pathways like phagosome, ribosome, olfactory transduction, and antigen processing and presentation found in this study suggested the need of further study. An unexpected discovery was silica-exposure induced differentially expressed mRNAs were involved in proteoglycans in cancer and pathways in cancer, and this result might reveal a critical link (direct or indirect) between silica exposure and lung cancer [29,30]. For instance, HIF1A was signi cantly associated with overall cancer risk [31], especially the lung cancer increased risk [32].
The following RT-qPCR results highlighted the importance of conducting silicosis related transcriptome study in human lung tissue, especially mRNA-related research. Based on lung tissues, RT-qPCR showed HIF1A, SOCS3, GNAI3, PTEN, hsa-miR-27a-3p, hsa-miR-27b-3p, hsa-miR-34b-3p, hsa-miR-575, hsa-miR-937-5p, and hsa-miR-181b-5p expressions were signi cantly down-regulated. In addition, PTEN had the most signi cantly down-regulated expression. This result was consistent with our previous study result: silicosis patients lung tissues had abnormal DNA methylation, and PTEN promoter hypermethylation might be associated with decrease of PTEN protein [15]. In addition, hsa-miR-8063 expression in both lung tissues and blood samples was down-regulated, which was not consistent with microarray analysis. This result might be explained by microarray detection potential limit such as limited representativeness [33]. However, PTEN and GNAI3 expressions were signi cantly up-regulated in blood samples, and the bisul te sequencing PCR result demonstrated PTEN had signi cantly decreased methylation rate in silicosis patients blood samples. Therefore, we speculated that there might be subsequent positive feedback regulation on PTEN decrease: lung tissue PTEN decrease promotes other tissues secrete PTEN into blood, then PTEN transferred from blood into lung tissues. Nevertheless, this study results provided critical reference for better clarifying the detailed biological mechanisms.
The major limitation of this study is the small number of lung tissues samples, precluding to some degree rm conclusions. However, Genome-wide pro ling based on silicosis patients and normal people lung tissues is inevitably limited by the small sample size. Besides, microarray-based analysis could not identify sequence variation, and the detection limit of the technology was also depended on the probe and array design [9].
In summary, we conducted a transcriptome study based on human lung tissues and blood samples, and broadened the subjects diversity in silicosis related transcriptome studies. Integrated analysis indicated PTEN hypomethylation could be a biomarker for early detection of silicosis.  a Hierarchical cluster of mRNA pro les. Every row represents one tissue sample: normal 1 to 4 represent normal lung tissues, early stage 1 and 2 represent early stage silicosis tissues, and advanced stage 1 to 5 represent advanced stage silicosis tissues. Every column represents one mRNA/miRNA probe. Red color indicates over-expression, and green indicates low-expression. The threshold of signi cance was P value < 0.05. b Hierarchical cluster of miRNA (b) pro les. Every row represents one tissue sample: normal 1 to 4 represent normal lung tissues, early stage 1 to 3 represent early stage silicosis tissues, and advanced stage 1 to 5 represent advanced stage silicosis tissues. Every column represents one mRNA/miRNA probe. Red color indicates over-expression, and green indicates low-expression. The threshold of signi cance was P value < 0.05. Figure 2 a Detailed trends of mRNA expression pro le 5.The horizontal axis represents different lung tissues: normal lung, early stage silicosis, advanced stage silicosis; and their vertical axis represents raw expression values that are converted into log2 ratio. b Detailed trends of miRNA expression pro le 8 (b).The horizontal axis represents different lung tissues: normal lung, early stage silicosis, advanced stage silicosis; and their vertical axis represents raw expression values that are converted into log2 ratio. Figure 3 miRNA-mRNA-network. Square represents microRNA, node represents mRNA, and straight line represents regulatory relationship between microRNA and mRNA. a RT-qPCR analyses of mRNAs relative expression levels based on lung tissues.Data show mean ± standard deviation, and the above asterisks between control and case groups represent statistical signi cance (unpaired two-tailed Student's t-test):*P < 0.05, **P < 0.01, ***P < 0.001. b RT-qPCR analyses of miRNAs (b) relative expression levels based on lung tissues. Data show mean ± standard deviation, and the above asterisks between control and case groups represent statistical signi cance (unpaired twotailed Student's t-test):*P < 0.05, **P < 0.01, ***P < 0.001. Figure 5 a RT-qPCR analyses of further selected mRNAs relative expression levels based on blood samples. Data show mean ± standard deviation, and the above asterisks between control and case groups represent statistical signi cance (unpaired two-tailed Student's t-test): *P < 0.05, **P < 0.01, ***P < 0.001. b RT-qPCR analyses of further selected miRNAs (b) relative expression levels based on blood samples. Data show