Cross-Species Transcriptomics Analysis Highlights Conserved Molecular Responses to Per- and Polyfluoroalkyl Substances

Abstract

In recent decades, per- and polyfluoroalkyl substances (PFASs) have garnered widespread public attention due to their persistence in the environment and detrimental effects on the health of living organisms, spurring the generation of several transcriptome-centered investigations to understand the biological basis of their mechanism. In this study, we collected 2144 publicly available samples from seven distinct animal species to examine the molecular responses to PFAS exposure and to determine if there are conserved responses. Our comparative transcriptional analysis revealed that exposure to PFAS is conserved across different tissues, molecules and species. We identified and reported several genes exhibiting consistent and evolutionarily conserved transcriptional response to PFASs, such as ESR1, HADHA and ID1, as well as several pathways including lipid metabolism, immune response and hormone pathways. This study provides the first evidence that distinct PFAS molecules induce comparable transcriptional changes and affect the same metabolic processes across inter-species borders. Our findings have significant implications for understanding the impact of PFAS exposure on living organisms and the environment. We believe that this study offers a novel perspective on the molecular responses to PFAS exposure and provides a foundation for future research into developing strategies for mitigating the detrimental effects of these substances in the ecosystem.

Keywords: PFAS; cross-species correlation; metabolism; transcriptome; transcriptomics.

Figure 1

Heatmap displaying the correlation among 110 different PFASs vs. control differential expression contrasts. The color gradient ranges from blue (denoting negative correlation) to red (denoting positive correlation), with darker colors indicating higher correlation values. Each colored dot indicates the correlation value between any two contrasts of the final signature matrix (Supplementary Table S3). The upper bar denotes the tissue of origin of each contrast.

Figure 2

(A) Correlation plot of M. musculus and H. sapiens exposed to the PFOA molecule [38,39,44]. As an example, the greatest correlation achieved between human and mouse contrasts (0.36) is highlighted by a black box. (B) Scatterplot showing the correlation among the contrasts of mouse and human highlighted by a black box in the previous plot. The highlighted genes are the most significant genes driving the correlation between the two species, defined by significant transcriptional change (p ≤ 0.001) in response to PFAS exposure in both species.

Figure 3

Comparison between PFAS exposure signatures in four different fish species [41,43,48,49]. For the Atlantic cod (Gadus morhua) dataset, one concentration for each PFAS molecule was selected. The full analysis including all concentrations and contrasts is displayed in Supplementary Figure S4.

Figure 4

Plot showing the integrated response to PFASs across 110 contrasts. Each point represents a gene. The x-axis indicates the integrated signature value (obtained by integrating signatures across the dataset using the Stouffer method). The y-axis indicates the standard deviation of the signature across the dataset. In red and orange, genes with the highest positive integrated signature (i.e., conserved PFAS-induced up-regulation across species), in blue and cyan, genes with the highest negative integrated signature (i.e., conserved PFAS-induced down-regulation across species). Genes in orange or cyan are also characterized by signature standard deviation above 10, indicating heavier fluctuations across the dataset (see also Figure 5 and Figure S5).

Figure 5

Line graph indicating the levels of expression of selected genes in response to PFAS molecules in different species. Each line is one gene: the genes shown here are the most consistently up- or down-regulated with low standard deviation, as extracted from the red and blue points of Figure 4. The x-axis reports all the 110 contrasts analyzed in the integrated dataset, grouped by species. The y-axis reports the signature for each gene, representing the significance (and sign) of the gene’s transcriptional response to PFASs. The horizontal lines delimit the p-value thresholds of 0.05.

Figure 6

Heatmap showing the 20 most significant pathways that are up-regulated and down-regulated across species. The blue-red color scale is proportional to the strength of the calculated pathway NES. White cells indicate contrasts with insufficient (<5) pathway genes to reliably calculate GSEA. The bottom bar indicates the species of each contrast in color code. The p-adjusted on the left side indicates the integrated p-value of pathway enrichment calculated across species.

Figure 7

Bar plot indicating the integrated normalized enrichment score of predicted metabolic changes upon exposure to PFASs in all human contrasts. TAG: triacylglycerol. LPC: lysophosphatidylcholine. PC: phosphatidylcholine. NAD: nicotinamide adenine dinucleotide.