Integrated Analysis of lncRNA and mRNA Transcriptomes Reveals New Regulators of Ubiquitination and the Immune Response in Silica-Induced Pulmonary Fibrosis

Objectives As an epigenetic player, long noncoding RNAs (LncRNAs) have been reported to participate in multiple biological processes; however, their biological functions in silica-induced pulmonary fibrosis (SIPF) occurrence and development remain incompletely understood. Methods Five case/control pairs were used to perform integrated transcriptomes analysis of lncRNA and mRNA. Prediction of lncRNA and mRNA functions was aided by the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Additionally, we constructed a coexpression network of lncRNAs and mRNAs to identify targets of regulation. Results In total, 1069 differentially expressed mRNAs and 366 lncRNAs were identified with the changes more than 2 times (p<0.05), of which 351 downregulated mRNA and 31 downregulated lncRNA were <0.5 (p<0.05) and those of 718 upregulated mRNAs and 335 upregulated lncRNA were >2 (p<0.05). The levels of 10 lncRNAs were measured via qRT-PCR; the results were consistent with the microarray data. Four genes named of FEM1B, TRIM39, TRIM32, and KLHL15 were enriched significantly with ubiquitination and immune response. Cytokine-cytokine receptor interaction was the most significantly enriched KEGG pathway in both mRNAs and lncRNAs. The coexpression network revealed that a single lncRNA can interact with multiple mRNAs, and vice versa. Conclusions lncRNA and mRNA expression were aberrant in patients with SIPF compared to controls, indicating that differentially expressed lncRNAs and mRNAs may play critical roles in SIPF development. Our study affords new insights into the molecular mechanisms of SIPF and identifies potential biomarkers and targets for SIPF diagnosis and treatment.


Introduction
Evidence has been presented that silica dust acute and chronic exposures will trigger an inflammatory cascade, such as macrophages and alveolar epithelial cells proliferated successively, followed by inflammatory cascade; the fibrotic development will be ignited, which eventually will proceed to determine the silicosis [1,2]. Silicosis is an irreversible and incurable pulmonary disorder; even when patients are no longer exposed to silica, the fibrosis remains progressive [3,4]. Silica dust is widely distributed in workplaces in which drilling, grinding, and hammering activities occur [5,6]. China has the highest silicosis burden worldwide, with more than 600,000 cases recorded over the past 30 years. Furthermore, the number of new cases is increasing, with over 24,000 deaths annually [7]. In South African gold mines, 2 BioMed Research International the incidence of silicosis is 13-25% in long-term miners [8]. Silicosis is also prevalent in developing countries. The CDC states that approximately 2 million US workers are currently exposed to silica; from 2001 to 2010, 1,437 deaths were attributable, in whole or in part, to silicosis; the youngest patient was aged only 19 years [9]. In the UK, over 3.2 million persons are occupationally exposed to silica [10]. Many preventative efforts have been made, including dust control and the development of personal protective equipment; in some developed countries, the incidence has thus steadily declined over the past few decades, but new cases or outbreaks are sporadically reported [11]. In developing countries, silicainduced silicosis associated with progressive lung fibrosis still remains a major concern [6,12]. No effective therapy is available, and the molecular details of fibrotic progression remain unclear [13]. Early biomarkers are urgently required; these would aid in implementation of preventative measures and allow for early diagnosis, intervention, and treatment.
Some biomarkers aiding in early diagnosis have been identified and can be divided into exposure, effect, and susceptibility biomarkers. Serum KL-6, MMP-2, and SP-D are potential biomarkers of silicosis, with levels being significantly higher in cases than controls (p<0.05) [14]. HO-1, the enzyme catabolizing heme to bilirubin, is a recognized biomarker [15]. Silica increases the ceruloplasmin level, which thus serves as a diagnostic biomarker [16]. The serum CC16 level indicates the risk of silica exposure toxicity [17]. However, current biomarkers studies regarding lung fibrosis induced by silica mostly focus on the genes level, seldom associated with transcriptome changes, which is an emerging concern which have been reported to be potential biomarkers in the development of lung fibrosis. Recently, epigenetic research has shown that the changes in microRNA (miRNA)/long noncoding RNA (lncRNA) expression may affect the stability and translation of genes involved in silicainduced lung fibrosis. MiRNAs are short noncoding RNAs that bind to targeting mRNAs in a sequence-specific manner [18]. For example, miR-489 was increased significantly in an animal model of lung fibrosis, suppressing TGF-1 synthesis (and where TGF-1 is a transcription factor precipitating the inflammation associated with lung fibrosis [19]). MiR-449a inhibits fibrosis by targeting Bcl2, in turn regulating autophagy [20]. miR-19a is downregulated in the early stages of fibrosis and was proposed as a biomarker for early diagnosis [21]. The miR-146a level in bronchoalveolar lavage fluid aids in early diagnosis [22]. lncRNA research only commenced in the last decade; lncRNAs are long noncoding RNAs [23] modulating a series of important biological processes, e.g., metabolism and apoptosis, and serve also as useful biomarkers of various diseases, including silica-induced lung fibrosis. Sun et al. showed that the lncRNAs suc.77 and 2700086A05Rik regulated the progression of lung fibrosis [24]. lncRNA-ATB, first identified in 2014, was activated by TGF-and induced lung fibrosis [25]. Thus, the potential utility of lncRNAs as fibrosis biomarkers requires attention, especially the investigation on human cases of silica associated fibrosis. However, such work in the context of silicainduced lung fibrosis is limited. Here, we collected peripheral blood samples from patients with silica-induced lung fibrosis followed by transcriptome analysis using a microarray; we then employed bioinformatics tools to identify changes in lncRNA levels in patients with progressive lung fibrosis; we identified the top 10 most-affected lncRNAs and further analyzed the interaction networks of lncRNAs and mRNAs. Our study may improve our mechanistic understanding of the silicosis and identifies novel diagnostic and therapeutic biomarkers.

Methods and Materials
. . Participants. We enrolled 10 subjects: 5 with phase I lung fibrosis (cases) and 5 healthy subjects who worked in the same industry but were not exposed to silica (controls). The inclusion criteria were (1) an occupational history of silica exposure and a fibrosis diagnosis based on the Chinese Pneumoconiosis Diagnosis Standards (GBZ70-2013), as agreed by at least three occupational physicians; (2) lack of complications such as tuberculosis, lung cancer, or an infection; and (3) voluntary agreement to participate. The mean age of the cases was 55±6.2 years and that of controls was 51±7.9 years; all participants were male. The silica dust concentration in the workplace was 6.43±0.77 mg/m 3 , exceeding the exposure limit of 0.2-1 mg/m 3 (GBZ2-2002). At least 5 mL of peripheral blood was collected from each subject into a PAX gene Blood RNA tube (BD Biosciences, San Jose, CA, USA), held at room temperature (22-25 ∘ C) for at least 2 h to completely lyse all cells, placed at -70 ∘ C and immediately delivered to the Beijing Kangpusen Biotech Company for microarray analysis. The study was performed in accordance with all relevant tenets of the Declaration of Helsinki and was approved by the Ethics Committee of Xiangya School of Public Health (Approval no. XYGW-2018-11). All participants provided written informed consent.
. . Total RNA Isolation. Total RNA was isolated using the mirVanaRNA Isolation Kit following the manufacturer's instructions. Briefly, a 1.25 volume of absolute ethanol was added to each blood sample followed by mixing; the mixture was passed through a filter and the filter was washed three times with 700 L of miRNA wash solution and 500 L general wash solution; then, RNA was eluted into 100 L of elution solution at 95 ∘ C and stored at -70 ∘ C until analyzed. RNA purity and concentration were assessed using a QIA-GEN RNeasy Mini Kit and by deriving a spectrophotometric ratio with the aid of a NanoDrop-2000 instrument, respectively (Thermo-Scientific, Waltham, MA, USA).
. . Microarray Analysis. In brief, microarray analysis was performed as follows: RNA/-poly-A-RNA-control mixtures were prepared by mixing test RNA samples (50-500 ng), a diluted poly-A-RNA control, and nuclease-free water; 4 L of first-strand buffer mix and 1 L of first-strand enzyme solution were then mixed and added, and first-strand was cDNA synthesized; this was followed by second-strand synthesis. Purified cDNA and second-cycle single-strand cDNA were subjected to WT cartridge hybridization, washed, and scanned using a DNA Microarray Scanner (Agilent, Santa Clara, CA, USA). Expression Console software (ver. 1.4.1) was used to adjust the raw data background and standardize microarray data.
. . Differential mRNA and lncRNA Expression and Cluster Analysis. Transcriptome Analysis Console (version 3.1, Affymetrix, Santa Clara, CA, USA) software was employed to explore the differential expression of mRNAs and lncRNAs between cases and controls; we selected RNAs exhibiting the greatest fold differences with reference to the p values. Cluster analysis [26] was used to identify genes exhibiting similar biological functions.
. . Gene Ontology and Pathway Enrichment Analysis. Gene Ontology (GO) and pathway enrichment analysis were performed to predict biological processes, cellular components, and molecular functions affected in disease. We used a controlled gene vocabulary (http://www.geneontology.org). The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was employed for pathway enrichment analysis.
. . lncRNA-mRNA Coexpression Network. To explore interactions among differentially expressed lncRNAs and mRNAs in cases and controls, we prepared coexpression networks using only differentially expressed mRNAs with Pearson correlation coefficients ≥0.95 (p<0.01) and generated visual data. We included five essential lncRNAs and all essential mRNAs in a network assessment. Only differentially expressed mRNAs with Pearson correlation coefficients ≥ 0.90 (p<0.01) were used to construct the lncRNA-mRNA network and generate visual representations.
. . Quantitative Real-Time PCR (RT-PCR) Validation. Total RNA was isolated and reverse-transcribed into cDNA using a PEXBIO RNA kit according to the manufacturer's protocol. qPCR was performed with the aid of Super Real PreMix Plus (SYBR Green, Tiangen, Beijing, China) on a CFX96TM Real-Time System (Bio-Rad, Hercules, CA, USA). GAPDH served as the internal control and fold change calculations were made using the 2-Δ Δ CT method. The primer sequences are shown in Table S1. Differences between mRNA and lncRNA expression levels were compared using Student's two-tailed t-test. A p value <0.01 was considered to reflect significance after false discovery rate correction for multiple testing.
. . Statistical Analysis. R software (R Development Core Team, Vienna, Austria) was used to compare the expression levels of mRNAs and lncRNAs, for cluster analysis, and to draw volcano maps. One-way ANOVA was employed to explore the significance of differences in mRNA and lncRNA levels after log 2 transformation. The fold changes between cases and controls were calculated; changes >2 or <0.5 were considered to indicate differential expression; p<0.05 was taken to indicate significance.

Results
. . Comprehensive Transcriptome Analysis Reveals Differentially Expressed mRNAs and lncRNAs. Through microarray analysis, 1069 differentially expressed mRNAs were identified with the changes more than 2 times (p<0.05), of which oneway ANOVA showed that the fold changes of 351 downregulated mRNA were < 0.5 (p<0.05) and those of 718 upregulated mRNAs were > 2 (p<0.05). The most prominently upregulated mRNA was TC0700007850.hg.1 (fold change, 8.00; p<0.05). The top 25 up-and downregulated mRNAs are listed in Table 1. Meanwhile, 366 differentially expressed lncRNAs were identified with the changes more than 2 times (p<0.05), of which one-way ANOVA showed that the fold changes of 31 downregulated lncRNAs were <0.5 (p<0.05) and those of 335 upregulated lncRNAs were >2 (p<0.05). The most prominently downregulated lncRNA was TC1300009322.hg.1, with a fold change of 0.395 (p<0.05). The most prominently upregulated lncRNA was TC1700010761.hg.1, with a fold change of 6.70 (p<0.05). The top 25 up-and downregulated lncRNAs are listed in Table 2.
Unsupervised cluster analysis of the mRNAs of the 10 samples, in terms of differentially expressed genes, yielded the results shown in Figure 1(a). A volcano plot of these mRNAs based on p values and fold changes is shown in Figure 1(b). Figure 1(b) shows that more mRNAs were upthan downregulated. Unsupervised cluster analysis of the lncRNAs of the 10 samples is shown in Figure 1(c). Highly expressed lncRNAs are shown in red and those expressed at low levels are in green. A volcano plot based on p values and fold changes is shown in Figure 1

. . Comprehensive Transcriptome Functional Analysis of Differentially Expressed mRNAs and lncRNAs Reveals New
Regulators of Ubiquitination and the Immune Response. To highlight the biological functions and signaling pathways change in differentially expressed mRNAs and lncRNAs, we performed GO analysis and KEGG pathway assessment on mRNAs and lncRNAs, respectively ( Figure 2). For mRNAs, in the main GO biological process affected was sensory perception of smell and interleukin-12 "---" presents that no symbols are available in the microarray databases, but they were actually presenting dysregulation in the case group.
secretion, while in the affected GO cellular components were principally the immune response including secretory dimeric IgA immunoglobulin complex and monomeric IgA immunoglobulin complex, while the GO molecular functions affected were olfactory receptor activity, cytokine binding, and, in particular, thiol-dependent ubiquitinyl hydrolase activity, ubiquitinyl hydrolase activity, and ubiquitin-like protein-specific protease activity (Figure 2(a)). Four genes named FEM1B, TRIM39, TRIM32, and KLHL15 were enriched significantly with ubiquitination and immune response. KEGG pathway analysis emphasized olfactory transduction, chemokine signaling, cytokine-cytokine receptor interaction, maturity-onset diabetes of the young, and bile secretion (Figure 2(b)). GO and KEGG analyses of differentially expressed lncRNAs are shown in Figure 3. Differentially expressed lncRNAs were predicted to have the following main functions: keratinization, keratinocyte differentiation, and epidermal cell differentiation in the biological process category; intermediate filament and the GO molecular components in the cellular component category; and G-protein-coupled amine peptide receptor activity, hormonal activity, transcription factor activity, and RNA polymerase/distal enhancer sequence-specific binding in the molecular function category (Figure 3(a)). KEGG pathway analysis of all differentially expressed lncRNAs significant at p<0.05 revealed that the affected pathways included olfactory transduction, neuroactive ligandreceptor interaction, estrogen signaling, the cholinergic synapse, renin secretion, and cytokine-cytokine receptor interaction (Figure 3(b)). In order to explore individual lncRNA function in the development of lung fibrosis, we subjected the following 10 lncRNAs to individual GO and KEGG analysis: 18S F, RPL41P3, RN7SL541P, RP11.537H15.3, RN7SL293P, RP4-620F22.2, Metazoa-SRP, RP11-609D21.3, RP11-138I1.4, and RN7SL783P. GO analysis revealed most of these selected individual lncRNAs of involvement in chemokine binding, CoA-transferase activity, and Gprotein-coupled serotonin receptor activity, all of which are associated with SIPF development (Figure S1-S5); intriguingly, some lncRNAs were predicted to have immune response such as RP11.537H15.3 affected immunoglobulin production and adaptive immune response. The KEGG data showed that these selected individual lncRNAs participated in multiple pathways, including the herpes simplex infection pathway, a hepatitis C pathway involving inflammation, and the PPAR signaling pathway.

. . Verification of Expression Changes in lncRNAs.
To verify the differential expression of lncRNAs, we quantified the expression levels of 10 lncRNAs in 20 cases and 20 controls, "---" presents that no symbols are available in the microarray databases, but they were actually presenting dysregulation in the case group. . . LncRNA-mRNA Network Analysis. To analyze the lncRNA-mRNA coexpression interaction network, we sought genes that were significantly enriched in both cases and controls. A total of 2,053 such genes were identified, and lncRNA-mRNA coexpression networks were constructed; this identified five lncRNAs (RPL41P3, RP11-609D21.3, RP4-620F22.2, AP001610.5, and RN7SL783P) that determined the essential differences between cases and controls ( Figures  S6). These lncRNAs and their coexpressed mRNAs, may indicate how silica exerts biological effects and may also serve as useful therapeutic targets. For example, AP001610.5 was coexpressed with the following mRNAs/genes: nuraro, IFIT , TAGLN , bawma, rasuru, IFIT , CT A , nymawbo, and CT A . RN7SL783P was coexpressed with many more mRNAs/genes, including pulmonary fibrosis-associated genes such as POYNEY and PODUBU. RP4-620F22.2, RPL41P3, and RP11-605F22.1 were also associated with many mRNAs (Table S2).

Discussion
Microarray and bioinformatics techniques now allow us to detect changes in lncRNA expression; most work to date has focused on miRNAs [27][28][29][30][31][32][33]. lncRNAs are longer sequence and have more complex functions than miRNAs [34]. For the first time, we herein identified mRNAs and lncRNAs that were abnormally expressed in human lung fibrosis induced by silica exposure, which exhibits complex pathological changes; no effective therapy is available. Biomarkers discovered via mRNAs and lncRNAs analysis might facilitate diagnosis and therapy [35,36]. Differentially expressed mRNAs and lncRNAs have been reported to be associated with development of lung fibrosis. was regulated at the mRNA level [38]. Certain mRNAs associated with lung inflammation (including that encoding IL-8) were appeared to be a regulator of other miRNAs biogenesis in lung epithelial cells, showing they were potential candidates for anti-inflammatory therapeutics for cystic fibrosis [39,40]. We found that mRNAs encoding LOYBOY, DEEZOY, and SNEYSPORBY were upregulated, and those encoding CPA , FCER A, and CCR were downregulated; these mRNAs may be involved in silica-induced lung fibrosis. Although lncRNAs may be potential therapeutic targets in fibrosis [41], the only relevant report is that of Sun et al.; in the paraquat-induced mouse lung fibrosis model, 513 lncRNAs were upregulated and 204 were downregulated, which affected cell differentiation, epithelial morphogenesis, and immune response, all of which are closely associated with the epidermal-mesenchymal transition (EMT) [24]. Our study identified 366 differentially expressed lncRNAs in human sample in which the numbers were lower than Sun et al. study. This might be due to the species variation.
Go and KEGG analyses are commonly used to predict biological functions and/or signaling pathways affected by mRNAs and lncRNAs. On GO analysis of significant changes in mRNAs (fold change >2; p<0.05), in terms of biological proses, cellular components, and molecular function, the most-affected processes were interleukin-12 secretion, regulation of intracellular mRNA localization, effects on the dimeric IgA immunoglobulin complex, C-C chemokine binding, inflammation response, immune response, ubiquitin-like protein-specific protease activity, and thiol-dependent ubiquitinyl hydrolase activity. Here, the results which the inflammation response affected by mRNAs are consistent with the previous report; however, mRNAs affecting the ubiquitination and immune response are new insights into the mRNAs function in the lung fibrosis. This suggests mRNAs might participate in lung fibrosis process by regulating immune response and ubiquitination. Ubiquitin is a highly conserved low-molecular-weight protein and ubiquitination represents a very common posttranslational modification. Extracellular ubiquitin regulates the immune response and exhibits anti-inflammatory and neuroprotective activities. Zhou et al. reported that silica activated macrophages and increased circRNA HECTD1 levels via ubiquitination [42]. Tsubouchi et al. suggested that, in a bleomycin-induced lung fibrosis mouse model, azithromycin (AZM; a second-generation antibacterial macrolide) inhibited autophagy via ubiquitination of NOX4 [43]. Nan et al. found that ubiquitination of the transcription factors Smad2 and Smad3 was changed via the action of ubiquitin carboxylterminal hydrolase-L5 during the pathogenesis of idiopathic pulmonary fibrosis [44]. However, the role of ubiquitination in the development of SIPF remains poorly understood and further studies are thus required. GO analysis indicated that both mRNAs and lncRNAs participate in the immune response in lung fibrosis. Balloy V. et al. assessed transcriptomes of lncRNA between cystic fibrosis and found a specific cystic fibrosis signature of lncRNA expression was associated with immune response of patients [45]. A study by Sookoian S. et al. found that lncRNA MALAT1 participated in the development of liver fibrosis [46]. However, studies regarding the underlying mechanism of lncRNAs participating in the development of silica-induced lung fibrosis remain further investigated. Liu et al. created an lncRNA/mRNA coexpression network for human pulmonary epithelial cells; 24 mRNAs were coexpressed and interacted with 33 lncRNAs; each mRNA interacted with a few lncRNAs and each lncRNA with a few mRNAs [47]. We included RPL41P3, RP11-609D21.39, RP4-620F22.1, AP001610.5, and RN7SL783P in a coexpression network with related mRNAs. For RPL41P3, GO analysis revealed involvement in chemokine binding, CoA-transferase activity, and G-protein-coupled serotonin receptor activity, all of which are associated with SIPF development ( Figure S1); lncRNA-RP11-609D21.3 was principally associated with carbohydrate transporter, endoribonuclease, and ribonuclease A activities, and with spliceosome activity, ferroptosis, and Kaposi's sarcoma-associated herpesvirus infection ( Figure S2). The lncRNARP4-620F22.2 was associated with sphinganine-1-phosphate metabolism, penetration of the zona pellucida, diol and sphingoid metabolism, sperm fusion to the egg plasma membrane, and viral RNA genome replication. Cell component analysis indicated that this lncRNA affected the exterior and marginal plasma membrane; the molecular functions involved were dihydrosphingosine-1-phosphate phosphatase, sphingosine-1-phosphate phosphatase, and hyaluronoglucosaminidase, as well as sequence-specific single-stranded DNA binding and hexosaminidase ( Figure S3). The biological processes in which Metazoa-SRP lncRNA was involved included type I interferon signaling pathway, the cellular response to type I interferon, negative regulation of helicase activity, the viral defense response and its negative regulation, intracellular viral protein transport, intracellular symbiont protein transport, and cellular responses to interferon-and exogenous dsRNA. The KEGG data showed that the lncRNA participated in multiple pathways, including the herpes simplex infection pathway, a hepatitis C pathway involving inflammation, and the PPAR signaling pathway ( Figure  S4). The lncRNA RN7SL783P molecular functions included melanocortin receptor binding and adenylate cyclase and phosphorus-oxygen lyase. The KEGG data showed that the lncRNA was active in melanogenesis, estrogen signaling, gastric acid secretion, GnRH signaling, circadian entrainment, aldosterone synthesis, and secretion, as well as inflammatory regulation of TRP channels, glucagon signaling, vascular smooth muscle contraction, oocyte meiosis apelin signaling, adrenergic signaling in cardiomyocytes, oxytocin signaling, cGMP-PKG signaling, calcium signaling, cAMP signaling, rap1 signaling, and phototransduction ( Figure S5).
Our work provides an initial understanding of the roles played by lncRNAs in SIPF. The lncRNA levels determined via quantitative reverse transcription (qRT-PCR) were consistent with the microarray data. Our results suggest that (1) silica exposure modulates the expression of both lncRNAs and genes and (2) some silica-induced genes may be involved in ubiquitination, and some lncRNAs may affect immune cell interactions. The mechanisms by which lncRNA and mRNA coexpression mediates SIPF require further study.
However, our study had certain limitations. First, the lncRNAs measured via qRT-PCR constituted only a small proportion of those identified via microarray; some lncRNAs may have been missed such that further studies are required. Second, we recruited only five cases; studies with larger samples are required. However, we are the first to use human peripheral blood samples to identify lncRNAs and mRNAs in patients with SIPF.

Conclusions
Gene and lncRNA expression profiles were dysregulated in the peripheral blood of SIPF patients, suggesting future diagnostic and therapeutic biomarkers. LncRNA dysregulation may play an important role in disease development, given its plethora of effects on posttranscriptional mRNA regulation. However, further biology function work is required.

Data Availability
The data used to support the findings of this study are included within the supplementary information file.

Additional Points
Key Messages. What is already known about this subject? Silica exposure in the workplace is the potential risk factor for workers, which has been identified to be the cause of silicosis, leading to damage workers health, with many complications including tuberculosis. Preliminary studies suggest that early discovery and diagnosis of biomarkers would benefit from providing interventions on improving patients' life quality. What are the new findings? To the best of our knowledge, this study focusing on the uncovering of the transcriptomes signature of patients with silicosis is firstly reported. Previous similar studies were conducted using animal model or cell model rather than using peripheral blood samples. In total, 1069 differentially expressed mRNAs and 366 lncRNAs were identified. Four genes named FEM1B, TRIM39, TRIM32, and KLHL15 were enriched significantly with ubiquitination and immune response. How might it impact clinical practice in the foreseeable future? This study indicated that these lncRNAs and mRNAs may play a critical biological role in the progression and development of silica-induced silicosis and are thus candidate therapeutic targets and potential diagnosis biomarkers. Meanwhile, this transcriptomes study might provide more baseline understanding the molecular mechanisms of silicosis which could pave the way for a further functional study of underlying regulatory mechanisms.

Ethical Approval
The study was performed in accordance with all relevant tenets of the Declaration of Helsinki and was approved by the Ethics Committee of Xiangya School of Public Health (Approval no. XYGW-2018-11).