Scientific Papers

Malignant peritoneal mesotheliomas of rats induced by multiwalled carbon nanotubes and amosite asbestos: transcriptome and epigenetic profiles | Particle and Fibre Toxicology

Description of Image

Genome-wide transcriptome analysis shows similarities and differences in tumors induced by MWCNTs and amosite asbestos

We carried out a genome-wide transcriptome analysis using Affymetrix microarrays on 11 tumors induced by MWCNTs or amosite asbestos, and 3 control peritoneal tissues. The principal component analysis (PCA) showed that the transcriptomes of the tumors grouped together, distinct from the control tissues (Fig. 1A). Hierarchical cluster analysis (unsupervised or on a few selected genes), also showed distinct clustering of tumor transcriptomes from those of control tissues (Fig. 1B, C). Visualization of the transcriptome data by Venn diagrams revealed both common and unique differentially expressed genes (DEGs) in the tumors induced by MWCNTs and amosite asbestos (Fig. 1D).

Fig. 1
figure 1

Transcriptome profiling of tumors induced by MWCNTs and amosite asbestos. A Principal component analysis (PCA) shows clear separation of tumors from the control peritoneal tissues. B Unsupervised hierarchical clustering analysis and heat map of 3062 filtered genes at fold change < − 2 or > 2, ANOVA P < 0.001, and FDR P < 0.001. C Hierarchical clustering and heat map on a few selected genes encoding for proteins such as mesothelin, osteopontin, caveolin, integrin, syndecan, chemokine, and interleukin. D Venn diagrams display the quantity of genes that are common and unique among the datasets of MWCNT B, MWCNT C, and amosite asbestos. The genes in the Venn diagrams were filtered using fold change < − 2 or > 2, ANOVA P < 0.05, and FDR P < 0.05. Quality control criteria and bioinformatics tools were according to Transcriptome Analysis Console (TAC 4.0.2, Thermo Fisher Scientific)

For further analysis of differential comparison of expressed genes between the tumors and control peritoneal tissues, we classified the tumors according to inducer or tumor type. The total number of DEGs in the tumors did not vary much between inducers (Table 1). On average, we found 9806 DEGs, of these 68.6% showed increased expression (upregulated), while 31.4% showed decreased expression (downregulated). When comparing results by tumor type, we found similar results. On average, there were 10,612 DEGs, of which 69.4% were upregulated and 30.6% were downregulated. The transcriptome data, when visualized through Venn diagrams, also showed both common and unique DEGs in the different tumor types (see Additional file 5: Figure S1A). From the total of 13,383 DEGs, 7866 (57.9%) were common to all three types of mesotheliomas. Sarcomatoid tumors had the highest number of unique DEGs with 1301 (9.6%), followed by epithelioid tumors with 1045 DEGs (7.7%), and biphasic tumors with 848 DEGs (6.2%). Furthermore, shown in Additional file 5: Figure S1B–D are the sample signals obtained for the genes encoding for forehead box M1 (Foxm1), mesothelin (Msln), and secreted phosphoprotein 1/osteopontin (Spp1) in the different tumor types and control peritoneal tissues. Sample signals, which represent the measured intensities of gene expression for each sample, were higher in all the tumor types than in the control peritoneal tissues.

Table 1 Summary of genome-wide transcriptome analysis in tumors by inducer or tumor type relative to control peritoneal tissues

To better understand the biological significance of the microarray data, we applied several bioinformatic approaches using Ingenuity Pathway Analysis (IPA). First, we searched for DEGs in the transcriptome datasets, which are implicated in mesothelioma or mesothelioma formation. In our analysis, we identified 38 DEGs that showed consistent changes in expression across all transcriptome datasets. This was true regardless of whether the datasets were grouped by inducer or by tumor type (Fig. 2A, Additional file 1: Table S1). For instance, mesothelin (Msln), osteopontin (Spp1), forkhead box MI (Foxm1), thymidylate synthetase (Tyms), and Wilms’ tumor suppressor gene (Wt1), were among those upregulated genes. Notably, ADAM metallopeptidase domain 10 (Adam10) and inhibin subunit beta A (Inhba) were also upregulated, and their upregulation was predicted to lead to the activation of mesothelioma formation (Fig. 2A).

Fig. 2
figure 2

Genes implicated in mesothelioma or mesothelioma formation. A Thirty-eight differentially expressed genes (DEGs) that showed consistent expression changes across the rat transcriptome datasets, regardless of inducer or tumor type. Overlay gene expressions (fold changes) represent those of MWCNT C. B The same set of genes is overlaid with gene expressions (fold changes) from the dataset GSE51024, which pertains to human malignant pleural mesotheliomas. Genes were filtered by fold change < − 1.5 or > 1.5, ANOVA P < 0.05, and FDR P < 0.05. Red (upregulated), green (downregulated), gray (did not meet at least one filter), orange (predicted activation)

To determine if the findings in rat mesotheliomas could be translated to human malignant mesotheliomas, we searched the Gene Expression Omnibus (GEO) database for relevant and publicly available gene expression datasets. However, we did not find any transcriptome datasets on human peritoneal mesotheliomas, specifically that have been obtained using Affymetrix microarrays and for which raw data are available. One possible explanation is that human pleural mesothelioma is more prevalent than peritoneal mesothelioma, resulting in more research being conducted on the former. Consequently, findings from pleural mesothelioma studies are often extrapolated to peritoneal mesothelioma [26]. After applying the same bioinformatics tools, we compared our rat transcriptome datasets with that of GSE51024 [27], which comprised 55 human malignant pleural mesotheliomas (MPM), and 41 lung parenchyma samples. Interestingly, we found that 17 of the 38 DEGs in the rat datasets were similarly expressed in the human MPM dataset (Fig. 2B, Additional file 2: Table S2). Furthermore, to determine whether similar results could be obtained using a different type of human tumor, we compared the rat transcriptome datasets with that of GSE149507 [28]. This dataset consisted of 18 pairs of small cell lung cancer (SCLC) tumors and adjacent lung tissues. In contrast to the human MPM dataset, where 17 of the 38 DEGs showed similar expression, only 9 of the 38 DEGs showed similar expression in the human SCLC tumor dataset (Additional file 6: Figure S2). These 9 DEGs, which were also similarly expressed in the human MPM dataset, included genes such as Spp1 and Foxm1 but not Msln, Adam10, and Inhba.

Second, we carried out Core and Comparative Analysis to identify biological changes across different transcriptome datasets of tumors. A graphical summary presented as a network of the key findings from the Core Analysis of the MWCNT B dataset is shown in Fig. 3A. This summary shows a relationship of a subset of the most significant entities predicted by the analysis, including canonical pathways, upstream regulators, diseases, and biological functions. In the transcriptome of tumors induced by MWCNT B, cytokines such as tumor necrosis factor (Tnf), interleukin 6 (Il6), interleukin 1 beta (Il1b), interleukin 33 (Il33), and interferon gamma (Ifng) were represented to be activated. Significant canonical pathways included Molecular Mechanisms of Cancer, ILK Signaling, Integrin Signaling, RAC Signaling, and TREM1 Signaling. Additionally, the molecular function colony formation was shown to be activated through the activation of forkhead box MI (Foxm1), SRC proto-oncogene (Src), and TEA domain transcription factor 4 (Tead 4).

Fig. 3
figure 3

Bioinformatic analyses of tumor transcriptomes induced by MWCNTs and amosite asbestos using Ingenuity Pathway Analysis (IPA). A Graphical summary obtained by Core Analysis of the MWCNT B transcriptome dataset. The summary shows a relationship of a subset of the most significant entities predicted in the analysis, which include canonical pathways, upstream regulators, diseases, and biological functions. IPA graph legends: orange (activation); blue (inhibition); solid line (direct interaction); dashed line (indirect interaction); dotted line (inferred edge); for signaling pathways: an arrow pointing from A to B signifies that A causes B to be activated. B After conducting Comparative Analysis, shown are highly-scoring predicted canonical pathways in the datasets as classified by inducers. Heat maps depict z-scores, in which orange (positive z-score) represents activation. The canonical pathways are sorted by hierarchical clustering and z-scores. The column clusters represent the inducers

After performing Comparative Analysis of transcriptome datasets grouped by inducers, we found that among the highly scoring common significant canonical pathways were Phagosome Formation, IL8 Signaling, Integrin Signaling, RAC Signaling, and TREM1 Signaling (Fig. 3B). Furthermore, when we carried out hierarchical clustering of these highly scoring canonical pathways and their z-scores, we observed that MWCNTs clustered differently from amosite asbestos. Characteristics of cancer cells and metastasis were among the highly scored molecular functions predicted to be regulated, regardless of classification by inducer (Fig. 4A) or tumor type (Fig. 4B). For instance, functions such as cell movement, migration, homing, invasion, and spreading were activated, while apoptosis was inhibited. Certain functions, however, were more activated in the tumors induced by MWCNTs than amosite asbestos (Fig. 4A). Hierarchical clustering showed clustering that reflected the tumor induction potential of the MWCNTs, and amosite asbestos observed during the in vivo study in rats [17]. Hierarchical clustering also reflected that sarcomatoid tumors are more aggressive than biphasic and epithelioid tumors (Fig. 4B).

Fig. 4
figure 4

Comparative Analysis of diseases and molecular functions of transcriptomes of tumors induced by MWCNTs and amosite asbestos. A Datasets are classified by inducers. B Datasets are classified by type of mesothelioma. Heat maps depict z-scores, in which orange (positive z-score) represents activation, while blue (negative z-score) represents inhibition. No color means the function being analyzed is not involved in the dataset. The column clusters represent the inducers or type of mesothelioma. Data analysis and heat map generation were performed using Ingenuity Pathways Analysis (IPA)

Furthermore, we carried out Upstream Analysis to identify master regulators across the transcriptome datasets. The activation or inhibition of a particular upstream molecule would lead to the gene expression changes observed in the datasets. Shown in Fig. 5A are examples of highly significant upstream regulators. These include tumor suppressor protein p53 (Tp53), colony stimulating factor 2 (Csf2), and E2F transcription factor 1 (E2f2). All these regulators were also predicted to be activated. We also performed Causal Network Analysis to further determine higher-level master regulators, which could have also influenced the gene expression changes in the transcriptome datasets of the tumors. In the analysis, Gadd45b was identified as a top master regulator. It was predicted to be activated with a z-score of 15. 299, 15.407, 14.769 and 14.275 for MWCNT B, MWCNT C, MWCNT D, and amosite asbestos datasets, respectively. The causal network with Gadd45b as a master regulator involved a depth of 2 and 11 other intervening regulators, that could helped explain the up- and downregulation of 630 genes in MWCNT B dataset (Fig. 5B, C, Additional file 3: Table S3). The canonical pathways that shared the highest overlap with the regulators in this causal network were Senescence Pathway (8 of 11) and Hepatic Fibrosis Signaling Pathway (7 of 11).

Fig. 5
figure 5

Upstream Analysis to identify novel master upstream regulators and causal networks in the transcriptome datasets of tumors induced by MWCNTs and amosite asbestos. A Highly scoring upstream regulators, ranked by P-values from right-tailed Fisher’s exact test. B Predicted upstream regulators of causal networks. Heat map depicts z-scores, in which orange (positive z-score) represents activation, while blue (negative z-score) represents inhibition. No color means the pathway or function being analyzed is not involved in the dataset. C A causal network with Gadd45b as an upstream regulator with 11 intervening regulators to explain differential expression of 630 genes in the transcriptome dataset MWCNT B. Shown are 13 target genes directly regulated by Jund. Graph legends: orange (activation); blue (inhibition); solid line (direct interaction); dashed line (indirect interaction); upregulated (red), downregulated (green), gray (did not meet at least one filter). Data analysis and graph generation were done using Ingenuity Pathways Analysis (IPA)

RT-qPCR also reveals gene expression changes in cancer-related genes in tumors induced by MWCNTs and amosite asbestos

To confirm the accuracy and robustness of our microarray data from the genome-wide transcriptome analysis, we performed RT-qPCR on a few selected genes that have been implicated in carcinogenesis. These genes were found be either downregulated: Hras-like suppressor (Hrasls), Nuclear receptor subfamily 4, group A, member 1 (Nr4a1), Fibroblast growth factor receptor 4 (Fgfr4), and Ret-proto-oncogene (Ret), or upregulated: Rho family GTPase3 (Rnd3/RhoE) and Growth arrest and DNA damage inducible beta (Gadd45b), by the microarray analysis. We analyzed 4 control peritoneal tissues and 36 tumors induced by the three MWCNTs or amosite asbestos. The RT-qPCR results corroborated the microarray data and provided further insight into the role of these genes in carcinogenesis. These results will be described next.

Hrasls is also known as Phospholipase A and acyltransferase 1 (Plaat1) among other aliases and belongs to the HRASLS (PLAAT) family consisting of five members (1–5) in humans, and three (1,3, and 5) in mice and rats [29]. The enzymes encoded by these genes are important in diverse biological functions, including tumor suppression, and organelle degradation [30]. Hrasls was downregulated in 36 tumors relative to the 4 control tissues (Fig. 6A).

Fig. 6
figure 6

RT-qPCR results of selected cancer-associated genes in tumors induced by MWCNTs or amosite asbestos. As compared to the control peritoneal tissues, mRNA expression was downregulated in Hrasls, Nr4a1, Fgfr4 and Ret (AD) or upregulated in Rnd3 and Gadd45b (E, F) in all or most of 36 tumors. Scatter dot plots and statistical results were obtained using GraphPad Prism 9. Statistical significance was determined at P < 0.05, t–test for unpaired values, two-tailed. Bar depicts mean and standard error of the mean (SEM). Asterisks depict statistical significance at **** (P < 0.0001), *** (P < 0. 001), **(P < 0. 01), * (P < 0. 05)

Nr4a1 encodes an orphan nuclear receptor that functions as a ligand-independent transcription factor. In human cancer, NR4A1 has a conflicting role, acting both as a tumor suppressor and an oncogenic driver [31]. For instance, in certain samples of human breast cancer, NR4A1 gene and protein expressions were decreased to suggest a tumor suppressor role [32]. Nr4a1 was downregulated in all 36 tumors relative to the 4 control tissues (Fig. 6B).

Fgfr4 encodes one of the four tyrosine kinase receptors for fibroblast growth factors (FGFs). FGF and their receptors (known as FGF/FGFR signaling) regulate various cellular functions, including tumor progression [33]. In a wide range of tumors in patients, FGFR4 is frequently overexpressed and exhibits sequence variations [34]. Fgfr4 was downregulated in 94% (34 of 36) of the tumors relative to the 4 control tissues (Fig. 6C). This discrepancy may be accounted to other genetic and or epigenetic variations present in the rat tumors, as will be shown later.

Ret, which was generated by recombination between two unlinked DNA fragments, encodes a transmembrane receptor with tyrosine kinase (RTK) activity [35]. RET rearrangements have been detected in a variety of human cancers. Ret was downregulated in 92% (33 of 36) of the tumors relative to the 4 control tissues (Fig. 6D).

Rnd3/RhoE belongs to a family of genes, which encode small G proteins with important role in tumor initiation and progression [36]. Rho GTPases regulate proliferation and apoptosis, metabolism, senescence, and cancer cell stemness. Rnd3/RhoE has been implicated in cancer cell motility and its mRNA and protein expression can be upregulated by DNA damage-inducing stimuli [37]. Rnd3 was upregulated in 36 tumors relative to the 4 control tissues (Fig. 6E).

Gadd45b is a member of the growth arrest DNA damage-inducible gene family. Members of this family have been implicated in many cellular functions such as DNA repair, cell cycle control, genotoxic stress response, and tumorigenesis [38]. In human cancer, GADD45B has been reported either as tumor suppressor gene or as oncogene [39]. Gadd45b was upregulated in 75% (27 of 36) of the tumors relative to the 4 control tissues (Fig. 6F).

Tumors exhibit CpG and non-CpG DNA methylation in gene promoter regions

Our RT-qPCR analyses showed that the mRNA expression of Hrasls, Ret, Nr4a1, and Fgfr4 was downregulated in all or most of the tumors we analyzed. To determine if DNA methylation, an epigenetic change, could be responsible for the repression of their transcription, we conducted bisulfite sequencing along the promoter regions of these genes in 11 tumors and 3 control peritoneal tissues. For each gene, we sequenced on average 8–10 clones from either tumor or control tissue, to identify methylated sites in the amplified fragments. For Hrasls, we analyzed a 431-bp fragment (28 CpGs, c.-376 to c.-46, relative to translation start site, ATG). Although overall percent CpG methylation between tumor and control tissues was not statistically significant, we obtained higher percentage of CpG methylation in certain tumors (Fig. 7A, B). Percent CpG methylation ranged from 3.13 to 12.3% in tumors, as compared to 3.03–4.64% in control.

Fig. 7
figure 7

DNA methylation in the promoter regions of genes that are downregulated in tumors induced by MWCNTs or amosite asbestos. A Examples of CpG methylation maps for the Hrasls and Ret genes obtained after Sanger bisulfite sequencing. Shown are methylation patterns of 8–10 sequenced clones from both control and tumor tissues. Methylated CpGs (filled lollipops); Unmethylated CpGs (unfilled lollipops). B, C Percent CpG methylation and the number of non-CpG methylation sites in 3 control and 11 tumor tissues in Hrasls, Ret, Nr4a1, and Fgfr4. Scatter dot plots and statistical results were obtained using GraphPad Prism 9. Statistical significance was determined at P < 0.05, t–test for unpaired values, two-tailed. Bar depicts mean and standard error of the mean (SEM)

For Ret, we analyzed a 468-bp fragment (37 CpGs, − 508 to − 56 relative to transcription start site, TSS). Similarly, there was no statistically significant differences in overall percent CpG methylation between tumor and control tissues. However, some tumors had a higher percentage of CpG methylation (Fig. 7A, B). Percent CpG methylation ranged from 2.95 to 13.78% in tumors, as compared to 2.7–4.05% in control. For Nr4a1, we analyzed a 500-bp fragment (47 CpGs, c.-530 to c.-89, relative to translation start site, ATG). In both tumor and control tissues, we observed very low methylation, averaging 2% (Fig. 7B). For Fgfr4, we analyzed a 391-bp fragment (19 CpGs, − 306 to 77, relative to transcription start site, TSS). Percent CpG methylation in both control and tumor tissues was 10% and below (Fig. 7B).

During bisulfite sequencing, we could also analyze for methylation at non-CpG sites (CpA, CpT, and CpC). Overall, there was no statistically significant differences in non-CpG methylation sites between tumor and control tissues (Fig. 7C). Nonetheless, in the promoter regions of Nr4a1, Fgfr4 and Ret, we observed more tumors with higher number of non-CpG methylation sites than in control tissues (Fig. 7C). To summarize, bisulfite sequencing of the promoter regions of Hrasls, Ret, Nr4a1 and Fgfr4 showed a higher percentage of CpG methylation and a higher number of non-CpG methylation sites in specific tumors when compared to control tissues. There was also observed inter-tumor heterogeneity in the DNA methylation of promoter regions, irrespective of inducer or tumor type (Additional file 7: Figure S3A–C).

Interestingly, along the promoter region of Fgfr4 in control tissues, we identified four sites involving a CpG (− 266 CpG), 2 sequence variations (− 226 T > C; − 124 T > C), and a non-CpG (-117CpT) relative to TSS. We observed concurrent methylation of the cytosine of all four sites in 70% (21 of 30) of clones in control tissues, as compared to 0.9% (1 of 108) of clones in tumors (Fig. 8A–C). This result suggests hypomethylation of these sites in tumor tissues. Furthermore, we identified a sequence variation (-38A > G) that was present in 62% (67 of 108) of clones in tumor tissues. In contrast, this variation was present in only 3% (1 of 30) of clones in control tissues.

Fig. 8
figure 8

Hypomethylation of four sites in Fgfr4 in tumors induced by MWCNTs and amosite asbestos. A Partial alignment of bisulfite sequences in control peritoneal tissues along the promoter region of Fgfr4, showing frequent methylation of all four sites: a CpG (− 266CpG), 2 sequence variations (− 226T > C; − 124T > C), and a non-CpG site (-117CpT). Methylated cytosines (blue boxes) are not converted into uracil by bisulfite treatment and remain unchanged. In contrast, similar alignment of bisulfite sequences in clones from tumors induced by MWCNT D did not show these frequently methylated sites. Instead, there is a high frequency of methylated cytosines in non-CpG sites. B Sequence electropherograms showing the four methylated sites in clones from a control peritoneal tissue (BF3). C Methylation frequency for the four sites in tumors and control tissues

Global DNA and RNA methylation levels are higher in tumors than in control tissues

Local and global changes in methylation can lead to cancer development. These changes are not limited to DNA but can also occur in RNA and histone proteins [40]. To determine further the role of methylation in the development of the tumors after exposure to MWCNTs, we analyzed global DNA and RNA methylation levels in 36 tumors and 4 control peritoneal tissues (Fig. 9). We found that percent 5mC DNA, m6A RNA, and 5mC RNA methylation levels were overall higher in tumors, regardless of the inducer or tumor type, than in control peritoneal tissues. The presence of 5mC methylation in RNA was confirmed by performing bisulfite sequencing on 28S ribosomal RNA. This was undertaken in 35 clones from 5 MWCNT-induced rat tumors, and in 13 clones from human lung carcinoma A549 cells (see Section “Materials and Methods”).

Fig. 9
figure 9

Global DNA and RNA methylation levels in tumors induced by MWCNTs or amosite asbestos. A, B Percent 5mC DNA methylation. C, D Percent m6A RNA methylation. E, F Percent 5mC RNA methylation. Results show global DNA and RNA methylation levels are higher in tumors than in control tissues. Scatter dot plots and statistical results were obtained using GraphPad Prism 9. Statistical significance was determined at P < 0.05, t–test for unpaired values, two-tailed. Bar depicts mean and standard error of the mean (SEM). Asterisks depict statistical significance at *** (P < 0. 001), **(P < 0. 01), * (P < 0. 05)

Description of Image

Source link