Small RNA-Seq Transcriptome Profiling of Mesothelial and Mesothelioma Cell Lines Revealed microRNA Dysregulation after Exposure to Asbestos-like Fibers

Environmental exposure to fibers of respirable size has been identified as a risk for public health. Experimental evidence has revealed that a variety of fibers, including fluoro-edenite, can develop chronic respiratory diseases and elicit carcinogenic effects in humans. Fluoro-edenite (FE) is a silicate mineral first found in Biancavilla (Sicily, Italy) in 1997. Environmental exposure to its fibers has been correlated with a cluster of malignant pleural mesotheliomas. This neoplasm represents a public health problem due to its long latency and to its aggression not alerted by specific symptoms. Having several biomarkers providing us with data on the health state of those exposed to FE fibers or allowing an early diagnosis on malignant pleural mesothelioma, still asymptomatic patients, would be a remarkable goal. To these purposes, we reported the miRNA transcriptome in human normal mesothelial cell line (MeT-5A) and in the human malignant mesothelioma cell line (JU77) exposed and not exposed to FE fibers. The results showed a difference in the number of deregulated miRNAs between tumor and nontumor samples both exposed and not exposed to FE fibers. As a matter of fact, the effect of exposure to FE fibers is more evident in the expression of miRNA in the tumor samples than in the nontumor samples. In the present paper, several pathways involved in the pathogenesis of malignant pleural mesothelioma have been analyzed. We especially noticed the involvement of pathways that have important functions in inflammatory processes, angiogenesis, apoptosis, and necrosis. Besides this amount of data, further studies will be designed for the selection of the most significant miRNAs to test and validate their diagnostic potential, alone or in combination with other protein biomarkers, in high-risk individuals’ liquid biopsy to have a noninvasive tool of diagnosis for this neoplasm.

Cells in confluent condition were separated from the culture flask (SPL Life Sciences; Korea) using 0.25% trypsin in 2.21 mM EDTA solution (Corning; Manassas, VA, USA) and counted using Bürker chamber by Trypan Blue Stain 0.4% (Gibco by Life Technologies; New York, NY, USA). The cells used for the experiments were between II/III passages.

In Vitro Treatments
FE fibers were obtained from the Biancavilla quarry (Sicily, Italy). These were sampled using magnifiers, needles and tweezers. Subsequently the fibers were weighed, sterilized under UV light for 10 min, suspended in a known volume of RPMI 1640 medium, and sonicated through Omni-Ruptor 4000 Ultrasonic Homogenizer (OMNI International Inc.; Kennesaw, GA, USA) for 10 min. The stock solution was then diluted appropriately to obtain the different concentrations for in vitro treatments.
The in vitro functional experiments were preceded by the determination of the doseresponse curves for both cell lines. MeT-5A were plated onto 96-well plates (Thermo Fisher Scientific; Roskilde, Denmark) for the dose-response curve at the density of 6 × 10 3 cells/50 μL, while JU77 were plated at the density of 4 × 10 3
Cells in confluent condition were separated from the culture flask (SPL Life Sciences; Korea) using 0.25% trypsin in 2.21 mM EDTA solution (Corning; Manassas, VA, USA) and counted using Bürker chamber by Trypan Blue Stain 0.4% (Gibco by Life Technologies; New York, NY, USA). The cells used for the experiments were between II/III passages.

In Vitro Treatments
FE fibers were obtained from the Biancavilla quarry (Sicily, Italy). These were sampled using magnifiers, needles and tweezers. Subsequently the fibers were weighed, sterilized under UV light for 10 min, suspended in a known volume of RPMI 1640 medium, and sonicated through Omni-Ruptor 4000 Ultrasonic Homogenizer (OMNI International Inc.; Kennesaw, GA, USA) for 10 min. The stock solution was then diluted appropriately to obtain the different concentrations for in vitro treatments.
The in vitro functional experiments were preceded by the determination of the doseresponse curves for both cell lines. MeT-5A were plated onto 96-well plates (Thermo Fisher Scientific; Roskilde, Denmark) for the dose-response curve at the density of 6 × 10 3 cells/ 50 µL, while JU77 were plated at the density of 4 × 10 3 cells/50 µL.  [10]. Results have been analyzed using PRISM GraphPad 7.00 and data were represented as the mean ± SD. An unpaired Student's t-test was used to compare data between the two groups. A value of p < 0.05 was considered statistically significant.
MeT-5A was plated at the density of 1 × 10 6 cells while JU77 was plated at the density of 8.5 × 10 5 cells onto 100 × 20 mm Petri Dishes (Eppendorf; Hamburg, Germany) for these functional experiments. After 24 h-incubation, the medium of both cell lines has been replaced with FE fibers solutions to final concentrations of 50 and 10 µg/mL [23]. MeT-5A and JU77 cells grown in normal medium were used as controls. After 48 h from FE fibers exposure, pellets were collected in duplicate.
After removing the supernatant, cells were harvested on ice by scraping in cold DPBS. Cells are then centrifuged at 0.2× g for 5 min at 4 • C and suspended in 1 mL cold DPBS. The cell solution was transferred to Eppendorf tubes and cells were centrifuged at 0.8× g for 5 min at 4 • C [23]. After removing the supernatant, samples were stored to −80 • C until RNA isolation.

RNA Isolation
Total RNA containing small non-coding RNA was extracted from the cell lines using miRNeasy Mini Kit (QIAGEN; Venlo, The Netherlands) according to the manufacturer's recommended protocols (miRNeasy Mini Handbook 11/2020). The RNA was quantified by the absorbance ratio at λ = 260/280 nm through NanoDrop (ND 1000) UV-Vis spectrophotometer [12]. All samples were diluted at the final concentration of 50 ng/µL for the Small RNA-Seq.

Small RNA-Seq
The QIAseq miRNA library kit (QIAGEN, Hilden, Germany) was used for small RNA-Seq library preparation following the manufacturer's instructions. RNA samples were quantified and quality tested by Agilent 2100 Bioanalyzer RNA assay (Agilent Technologies, Santa Clara, CA, USA). Libraries were then checked with both Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) and Agilent Bioanalyzer DNA assay or Caliper (PerkinElmer, Waltham, MA, USA). Finally, libraries were prepared for sequencing and sequenced on single-end 150 bp mode on NovaSeq 6000 (Illumina, San Diego, CA, USA) [23].

Small RNA-Seq Analysis
Low-quality reads and adapters were trimmed using Trim Galore, which is a wrapper of FASTQC [24] and Cutadapt [25]. After that, trimmed reads were aligned onto the reference human genome (HG38 version) using HISAT2 [26]. The generated SAM files were first converted into BAM files using Samtools; secondly, mapped reads were quantified by feature Counts (parameters: -d 14 -primary) using a custom GTF annotation file containing the genomic coordinates of miRNAs (reported in miRBase [27]). The obtained raw count values were normalized to scale the raw library sizes in trimmed mean of M values (TMM) by using edgeR [28] and all miRNAs, whose geometric mean of TMM values across all samples was less than one, were removed from the analysis, because they were nonexpressed or expressed at a very low level. Finally, the filtered count matrix was used for the differential expression analysis using LIMMA [29]. miRNA with a |Log2FC| > 0.58 and an adjusted p-value (Benjamini-Hochberg correction) <0.05 were considered differentially expressed. Finally, the impact of differentially expressed miRNAs on biological pathways was evaluated by using MITHrIL [30]. In this case, we selected all miRNAs with significant adjusted p-values without a Log2FC cutoff in order to evaluate also the impact of slightly differentially expressed miRNAs on biological pathways. All the above-mentioned steps of the analysis were performed using RNAdetector [31] software.

Choice of Exposure Time
The in vitro functional experiments were preceded by the determination of the doseresponse curves for both cell lines. The results showed a significant difference between 6 h treatments and all other time points (24, 48, and 72 h) at all tested concentrations exposure to FE fibers. Finally, an exposure time to the fibers equal to 48 h was chosen because longer exposure times did not cause statistically significant differences in cell viability (Figures 2 and 3). biological pathways was evaluated by using MITHrIL [30]. In this case, we selected all miRNAs with significant adjusted p-values without a Log2FC cutoff in order to evaluate also the impact of slightly differentially expressed miRNAs on biological pathways. All the above-mentioned steps of the analysis were performed using RNAdetector [31] software.

Choice of Exposure Time
The in vitro functional experiments were preceded by the determination of the doseresponse curves for both cell lines. The results showed a significant difference between 6 h treatments and all other time points (24, 48, and 72 h) at all tested concentrations exposure to FE fibers. Finally, an exposure time to the fibers equal to 48 h was chosen because longer exposure times did not cause statistically significant differences in cell viability (Figures 2 and 3).   biological pathways was evaluated by using MITHrIL [30]. In this case, we selected all miRNAs with significant adjusted p-values without a Log2FC cutoff in order to evaluate also the impact of slightly differentially expressed miRNAs on biological pathways. All the above-mentioned steps of the analysis were performed using RNAdetector [31] software.

Choice of Exposure Time
The in vitro functional experiments were preceded by the determination of the doseresponse curves for both cell lines. The results showed a significant difference between 6 h treatments and all other time points (24, 48, and 72 h) at all tested concentrations exposure to FE fibers. Finally, an exposure time to the fibers equal to 48 h was chosen because longer exposure times did not cause statistically significant differences in cell viability (Figures 2 and 3).

Differential Expression Analysis
In order to identify dysregulation in miRNAs induced by FE fibers, we performed an RNA-Seq transcriptome profiling of unexposed and exposed normal mesothelial (MeT5A) and malignant mesothelioma (JU77) cell lines. The differentially expressed miRNAs were: 333 untreated JU77 vs. MeT5A (145 up-regulated and 188 down-regulated), 323 in JU77 vs. MeT5A exposed to 10 µg/mL FE fibers (101 up-regulated and 222 down-regulated), and 325 in JU77 vs. MeT5A exposed to 50 µg/mL FE fibers (124 up-regulated and 201 down-regulated) ( Table 1). Among these, 14 were in common between the JU77 and MeT5A untreated and exposed to 10 µg/mL FE fibers (6 up-regulated and 8 down-regulated), 21 were in common between the JU77 and MeT5A untreated and exposed to 50 µg/mL FE fibers (18 up-regulated and 3 down-regulated), and 19 were in common between the JU77 and MeT5A exposed to 10 and 50 µg/mL FE fibers (2 up-regulated and 17 down-regulated). Among all samples, the differentially expressed miRNAs in common were 54 (Table 2). Table 1. Interspecies differentially expressed miRNAs.

Down-Regulated miRNAs in Common
JU77NTvsMeT5ANT and JU10vsMeT10 JU77NTvsMeT5ANT and JU50vsMeT50

Pathways Analysis
Once the differentially expressed miRNAs for each comparison were identified, we investigated the impact of their dysregulation in metabolic and signaling pathways by using MITHrIL [30]. MITHrIL fully exploits the topological information encoded by pathways when computing perturbation scores. Pathways are then modeled as complex graphs where each node is a biological element (protein-coding gene, miRNA, or metabolite) and each edge is an interaction between them [30]. Importantly, MITHrIL takes into account experimentally validated miRNA-mRNA interactions so as to predict their effects on biological pathways [30].
The results demonstrated clear patterns of positive and negative perturbation scores involving 96 different pathways between the JU77 and MeT5A untreated and exposed to 10 and 50 µg/mL FE fibers. Of these, 28 showed mildly positive perturbation scores, while 68 showed mildly strong negative perturbation scores. Among these latter pathways, strong, negative patterns were observed in 25 pathways, including apoptosis; regulation of lipolysis in adipocytes; endocytosis; ovarian steroidogenesis; focal adhesion; progesteronemediated oocyte maturation; axon guidance; adherence junction; longevity regulating pathway; regulation of actin cytoskeleton; sphingolipid signaling pathway; chemokine signaling pathway; leukocyte trans-endothelial migration; signaling pathways regulating pluripotency of stem cells; and specific signaling pathways of Rap 1, cGMP-PKG, AMPK, Hippo, PI3K-Akt, TGF-beta, ErbB, cAMP, Ras, and FoxO ( Figure 4).
The pathway analysis performed between untreated vs. FE fibers treated JU77 and MeT-5A and the main impacted pathways showed clear patterns of positive and negative correlations involving 121 different pathways. Of these, 66 showed mildly strong, positive correlations and, among these 25, showed strong, positive correlations including signaling pathways of FoxO, PI3K-Akt, Rap1, VEGF, p53, and Ras; the regulation of actin cytoskeleton; signaling pathways of chemokine, prolactin, sphingolipid, estrogen, insulin, thyroid hormone, and oxytocin; long-term depression; progesterone-mediated oocyte maturation; natural killer cell-mediated cytotoxicity; Gap junction; serotonergic and cholinergic synapse; axon guidance; signaling pathways regulating the pluripotency of stem cells; longevity regulating pathways; and focal adhesion. Mildly strong, negative correlations were observed for 55 patterns, and among these, 14 showed strong, negative correlations, including metabolic pathways of arginine and proline, alanine, aspartate, and glutamate; glycolysis/gluconeogenesis; signaling pathways of adipocytokine; NF-kappa B; Toll-like receptor; glucagon; HIF-1; Jak-STAT; osteoclast differentiation; glycosaminoglycan degradation; cytokine-cytokine receptor interaction; and the adherence junction. The correlations between these pathways were quite obvious for the JU77 cell line untreated vs. treated with FE fibers ( Figure 5).  The pathway analysis performed between untreated vs. FE fibers treated JU77 and 5A and the main impacted pathways showed clear patterns of positive and neg correlations involving 121 different pathways. Of these, 66 showed mildly strong, po correlations and, among these 25, showed strong, positive correlations including sign pathways of FoxO, PI3K-Akt, Rap1, VEGF, p53, and Ras; the regulation of actin cytoske signaling pathways of chemokine, prolactin, sphingolipid, estrogen, insulin, th hormone, and oxytocin; long-term depression; progesterone-mediated oocyte matur natural killer cell-mediated cytotoxicity; Gap junction; serotonergic and cholinergic syn axon guidance; signaling pathways regulating the pluripotency of stem cells; long receptor; glucagon; HIF-1; Jak-STAT; osteoclast differentiation; glycosaminoglycan degradation; cytokine-cytokine receptor interaction; and the adherence junction. The correlations between these pathways were quite obvious for the JU77 cell line untreated vs. treated with FE fibers ( Figure 5).

Discussion
Small RNA-Seq transcriptome profiling of healthy mesothelium and MPM in vitro has been evaluated to highlight the deregulated miRNAs and the various pathways involved in an aggressive cancer such as MPM.
Certainly, the results showed a big difference in the number of deregulated miRNAs between tumor and nontumor samples both exposed and not exposed to FE fibers. As a matter of fact, the effect of exposure to FE fibers is more evident in the expression of miRNAs in tumor samples than in nontumor samples.
In the present paper, we analyzed several pathways that are involved in the pathogenesis of MPM. It is very interesting to point out the strong perturbation scores involving the above-mentioned pathways in the MPM vs. healthy mesothelial cell line. Among these, certain pathways were involved that play important roles in inflammatory processes and angiogenesis. Inflammation has a central role, since mesothelioma is a multicentric neoplasm that originates from inflammatory foci. Inflammation has been correlated with cancer, enhancing the development of malignancies [32]. In particular, the chemokine and TGF-beta signaling pathways lead to acute and chronic inflammation, the latter resulting in several fiber-associated pulmonary and pleural diseases [16,33]. The inflammasome is responsible for the activation of inflammatory processes via multiple mechanisms [33] that trigger a cell death process called proptosis, characterized both by apoptosis and necrosis. Cell death mechanisms and the release of chemokines and cytokines may help cancer regress and resist toxicity by fibers and cell growth [34] in inflammasome-dependent and -independent pathways [35]. Indeed, apoptosis after exposure to FE fibers is a mechanism meant to remove irreparably damaged cells, which causes genetic changes that predispose cells to neoplastic transformation [36]. Specific signaling pathways such as sphingolipid, FoxO, and Hippo pathways are involved in many important signal transduction processes, such as cell proliferation and apoptosis [37][38][39]. Dysregulation of the Hippo signaling pathway is highly conserved by phosphorylating and inhibiting the transcription coactivators YAP and TAZ, key regulators of proliferation and apoptosis. On the contrary, dephosphorylated YAP/TAZ moves into the nucleus and activates gene transcription through binding to the TEAD family and other transcription factors. Such changes in gene expression promote cell proliferation and stem cell/progenitor cell self-renewal but inhibit apoptosis, thereby promoting tissue regeneration, and tumorigenesis [38]. An experimental model demonstrated the activation of YAP caused by ATG7 deletion [40], which is an important transcription activator in malignant mesothelioma [41]. In our recent research, ATG7's high expression represents a promising prognostic tool for patients with MPM [42]; thus, it would be interesting to explore whether there is an inverse correlation between ATG7 and YAP in malignant mesothelioma. The leukocyte migration to the site of injury is orchestrated by chemokines. Neutrophils are the first to be recruited onto the injury site, followed by monocytes, which differentiate into macrophages. Once activated, macrophages are the main source of growth factors and cytokines that affect the local microenvironment. Mast cells also contribute to inflammatory mediators, such as histamine, cytokines, and proteases, as well as lipid mediators [43]. In previous studies of our research group [16,44], we demonstrated the involvement of cytokines IL-18 and IL-1beta in the inflammasome activation process, suggesting that these immune modulators are involved in the pathogenic mechanisms triggered by FE fibers. Acute and chronic inflammations often generate common molecular mediators [37]. The strong inflammation, with an even increased angiogenesis, causes the fatal outcome of the neoplasm. The release of angiogenic cytokines, including TGF-beta and VEGF, occurs during the angiogenesis process in the malignant mesothelioma progression [45]. In particular, VEGF represents the main angiogenic cytokine involved in this cancer [46,47], modulating also the development of pleural effusion and ascites through an enhanced vascular permeability [47]. VEGF activation plays an essential role in increasing the survival of normal cells exposed to carcinogenic agents. FE fibers are able to induce functional modifications of parameters with crucial roles in cancer development and progression [11]. The synthesis of VEGF and beta-catenin, two critical steps of epithelial cell activation pathways, is affected by FE fibers exposure shown by an abnormal cellular status with up-regulated cell activities and a risk of neoplastic transformation [48]. Furthermore, the influence of FE fibers on cell motility has been demonstrated through a dysregulated and altered distribution of actin network [48]. Focal adhesion, which forms mechanical links between the cytoskeleton and extracellular matrix (ECM) and adherence junction, are clearly involved in malignant mesothelioma pathogenesis. About that, our recent study on FE exposure in lung fibroblasts suggested an ECM remodeling that can give rise to profibrotic cellular phenotypes and the tumor microenvironment [49]. The signaling pathway of natural killer (NK) cells mediated cytotoxicity has been shown to be involved in the development of malignant mesothelioma. However, the absence of NK cells does not alter tumor growth rates, suggesting that they cannot function as effector cells in this microenvironment. However, after local IL-2 and/or anti-CD40 antibody treatment of mesothelioma tumors, NK cells help acquire and/or maintain systemic immunity and long-term effector/memory responses [50].
Specific signaling pathways have been found to be involved in malignant mesothelioma. Among these, Ras and p53, the commonly mutated genes associated with cancer, are rarely targeted in malignant mesothelioma [51]. Ras has not been shown to be altered in mesothelioma cell lines [52,53]. However, several receptor tyrosine kinase pathways have been shown to be triggered in mesothelioma, including the epidermal growth factor receptor (EGFR), insulin-like growth factor receptor (IGFR), and c-Met [54][55][56], all of which activate Ras signaling. According to these results, several studies have already suggested that the PI3K-Akt pathway is hyperactivated in mesothelioma cell lines [51,57], resulting in the gain or loss of function of its downstream proteins, 4E-BP1 and pS6, both crucial to the regulation of protein synthesis [58]. However, the prognostic role of the PI3K pathway in MPM is not as yet quite defined [59].

Conclusions
Our goal was the validation of the most promising results in a subset of patients chronically exposed to FE fibers, using the liquid biopsy, to provide a minimally invasive screening tool for the secondary prevention of MPM. Early detection of circulating tumor biomarkers represents one of the most promising strategies to enhance the survival of cancer patients by increasing treatment efficiency [60]. Besides this large amount of data, further studies will be necessary in order to select the most significant miRNAs to test and validate their diagnostic potential, alone or in combination with other protein biomarkers, in high-risk individuals.