A bulk RNA sequencing was performed to identify transcriptional alterations in the liver. Global transcriptome screening identified 11,928 and 16,826 differentially expressed transcripts in the liver of BPA lineage males and females compared to the control group. Based on the volcano plot, the livers of the BPA lineage males (Supplementary Fig. 1A) displayed a significantly decreased number of unique DEGs than the livers of the BPA lineage females (Supplementary Fig. 1B). Using strict selection criteria (5 < log2FC< -4 and for females and 2.5 < log2FC< -2 for males), mRNA biomarkers associated with the transgenerational NAFLD phenotype were identified in males (Fig. 1A) and females (Fig. 1B) of the BPA lineage. As compared to the control group (Fig. 1A), elovl5, igfbp1, tlr5, hck, and apcs were significantly upregulated, and elovl1, fbxo4, invs, prkaca, and pheta2 were significantly downregulated in the liver of the BPA lineage males. The top 10 up- and down-regulated genes associated with direct exposure to BPA in the male liver were provided in Supplementary Table 1A and B71. In the female liver, vtg3, fabp7, cacna2d4, esr1, aldh18a1, mttp, fas, and pycr1 were significantly upregulated while pparα, acox3, cdlk2, enpep, trak1 and ugt3a1 were significantly downregulated in the BPA lineage than the control lineage (Fig. 1B). A list of significantly up- and down-regulated DEGs in females due to direct BPA exposure was provided in Supplementary Table 1C and D71. As the metabolic pathways are common among medaka, mice, and humans43, the DEGs (biomarkers) induced by BPA exposure (direct, intragenerational) in the mouse liver were compared between the two vertebrate species- mice and medaka. Results suggest that transgenerational liver disease biomarkers are significantly different from those caused by direct exposure to BPA.
Transcription factors-mediated cancerous pathways were found in the livers from BPA lineage females
To gain insight into how transcription factors play a role in BPA-induced transgenerational NAFLD, we examined the expression of several transcription factors and their roles in triggering disease-specific pathways in the liver of BPA lineage medaka. In the livers of BPA lineage male, srf, sall1, rad21, elf1, mef2c were collectively formed transcription factors (TFs) network controlling expression of genes associated with NAFLD phenotype (Supplementary Fig. 2A and Supplementary Table 2). TF networks and their expression profiles in the liver of BPA lineage females were shown in Supplementary Fig. 2B and heatmaps (Fig. 1C), respectively. In the liver of the BPA lineage females, transcription factors, mainly xbp1, myc, atf4, stat 4, and stat 6 were upregulated, whereas hnf4a, ppargc1a, klf5, sp1, and vdr were downregulated (Fig. 1C). The number of disease-specific TFs found in the livers of BPA lineage females was significantly increased compared to the livers of BPA lineage males. Transcription factors found in the female livers of the BPA lineage were associated with transcriptional dysregulation in cancer, estrogen signaling pathways, Jak-STAT signaling pathways, Th1 and Th2 cell differentiation, and TGF beta signaling pathways in cancer (Fig. 1D). Ingenuity pathway analysis (IPA) determined 50 cancerous genes (Supplementary Table 3) expressed in the female livers of the BPA lineage compared to the control lineage. Estrogen receptor-mediated (esr1 and esr2) activation of myc, rara, and ccnd1 and downstream activation of the cross-linking network were provided in Fig. 1E and Supplementary Fig. 3.
Both male and female livers of the BPA lineage expressed upregulated genes related to the innate immune pathway, LDL/HDL-mediated lipid transport, and downregulated genes related to lipid digestion
To determine the up and downregulated DEGs and their associated pathways contributing to the NAFLD pathogenesis, significantly altered genes in the liver of BPA lineage were screened by considering significance criteria |log2FC| > 0.5, and FDR < 0.1. In total, there were 115 and 1012 significantly upregulated genes and 387 and 702 significantly downregulated genes in the liver of BPA lineage males (Fig. 2A) and females (Fig. 2B), respectively. Out of all upregulated genes, BPA lineage had 11 common in both sexes, 104 male-specific, and 1001 female-specific DEGs (Fig. 2A and Supplementary Table 4A). Out of all downregulated genes, 85 were common in both sexes, 302 were male-specific, and 617 were female-specific (Fig. 2B and Supplementary Table 4B). As fatty liver disease phenotypes developed in both males and females from BPA lineages, we explored pathways triggered by up- and down-regulated genes common in both sexes. According to KEGG analysis, commonly upregulated genes were enriched in chemokine and chemokine receptor activation, triggering innate immunity (Fig. 2C), suggesting immunogenic pathways being impacted. In addition, down-regulated genes common in males and females indicated impact in lipid digestion, mobilization and transport, lipoprotein metabolism, HDL/LDL-mediated lipid transport, and bile acid synthesis (Fig. 2D). Supplementary Tables 5 and 6 demonstrated the enrichment of total pathways associated with DEGs in the BPA lineage males and females. The present results indicated that up- and downregulated genes in both sexes of BPA lineage fish contributed to dysregulation of immunogenic and lipid clearance pathways in the liver.
Sexually dimorphic expression of genes associated with fat metabolism and immunogenic genes in the BPA lineage fish.
To ascertain the role of fat metabolizing genes in unequal and sexually dimorphic fat deposition in the liver of the BPA lineage fish, we screened for the expression pattern of DEGs associated with cholesterol metabolism, fatty acid transport, lipogenesis, and lipolysis. Heatmap was used to demonstrate the expression pattern of fat metabolizing genes in the liver of BPA lineage males (Fig. 3A) and females (Fig. 3B). The expression of genes encoding enzymes for the lipogenic cycle, mainly sqle, mttp, cers5, and aco1 was significantly downregulated in the liver of the BPA lineage males compared to the control females (Fig. 3A). However, lipogenic genes such as scd, mttp, and abcg1 were significantly upregulated in the liver of BPA lineage females compared to controls, except for pnpla3 (Fig. 3B). This suggested that lipogenesis was elevated in females' livers than males of the BPA lineage. In contrast, genes encoding enzymes controlling the lipolytic cycle were significantly downregulated in both males and females of the BPA lineage. Downregulated lipolytic genes included acadm, cpt1ab, and pparδ in the males and pparα, pparδ, cyp 450, and crot in the females, suggesting that the lipolytic cycle was disrupted in the liver of BPA lineage fish. The mRNAs for srebf1, srebf2, pparγ, and pparggc1a were significantly decreased in the liver of BPA lineage females compared to the control females. Genes associated with fatty acid transport mechanisms, such as osbp, cetp, lrpap1, and slc27a1 in the male liver and apoab and fabp3 in the BPA lineage female livers, were significantly downregulated, suggesting restricted fat release from the liver in the BPA lineage fish. In contrast to aberrant expression of fat transport genes, the receptor protein CD36, which controls the intake of extrahepatic fat granules, was significantly increased in the liver of BPA lineage females compared to control females. The receptor protein cd36, which controls the intake of extrahepatic fat granules, was significantly upregulated in the liver of BPA lineage females, suggesting an increased intake of fatty acids rather than an effective fat clearance process. In the liver of both BPA lineage males and females, the DEGs associated with cholesterol metabolism were downregulated, but soat2 was upregulated in the liver of the BPA lineage males compared to control males. BPA lineage males showed downregulation of acot7, elovl6, and elovl5, while elovl1 was upregulated. The expression of tlr5, tnfaip3, cxcl12, hla-dqa1, rsl24d1, and cxcr4 were considerably increased in the liver of BPA lineage males compared to control males (Fig. 3A). cxcr4, cxcl12 were upregulated and irf3, caspase3, caspase8, ap1m2, il10, tlr8 and irak1bp1 were upregulated in the liver of BPA lineage females compared to controls (Fig. 3B).
The liver of BPA lineage females showed abnormal expression of genes associated with estrogen signaling, epigenetic processes, aberrant cell cycle regulation, activation of kinases, and cytokines
In the BPA lineage females, estrogen signaling pathway genes such as adcy9, cacna1s, cacna1d, and cacng8 were downregulated except for cacna2d4 (Fig. 3B). In addition, genes encoding for estrogen-dependent receptors, including esrp1, esr1, cav2, and cavin1b, were significantly upregulated, whereas esr2 and esrp2 were downregulated (Fig. 4A). A list of genes involved in estrogen receptor-mediated activation of downstream genes is shown in Supplementary Fig. 3. To explore the involvement of epigenetic processes in disease phenotype, we examined genes encoding enzymes involved in epigenetic modification in the livers of the BPA lineage females, which showed increased NAFLD phenotype. The genes related to DNA demethylation, mainly tet1, tet2, and tet3, were significantly downregulated, and all subtypes of dnmts, pcna, and uhrf1 were significantly upregulated in BPA lineage females than control females (Fig. 4B). The genes encoding histone deacetylation, mainly hdac1, hdac5, hdac7b, hdac11, hdac4 were significantly downregulated, but hdac9b, hdac8, hdac3, hdac12, hdac6, and sirt were upregulated (Fig. 4B). This suggested that epigenetic dysregulation in the liver of BPA lineage females due to ancestral BPA exposure. The genes encoding kinases, including map3k21, ripk2, mapk1, and mtor were differentially expressed (Fig. 4C), and cytokine genes were significantly upregulated in the liver of BPA lineage females than control livers (Fig. 4D). The cell cycle regulation genes wee1, aruka, rrm2, and rrm2 were significantly upregulated in the liver of BPA-lineage females compared to control females. The cell cycle regulation genes wee1, aruka, rrm2, yme1|1 and rrm2 were significantly upregulated in BPA-lineage female livers compared to control female livers (Fig. 4E) and these genes were associated with G1/S, G2/M transition, APC/C mediated degradation of cell cycle protein, DNA replication (Supplementary Fig. 4A).
Impact on ER and mitochondrial genes in developing oxidative stress response in the liver of BPA lineage females
When compared between females in the BPA and control lineages, endoplasmic reticulum-mediated stress response genes, mainly Xbp1, hspa5, hyou1, atf4, atf6, and ddit3 were significantly upregulated in the liver of the BPA lineage females (Fig. 4F). Among all ER stress genes, atf4, hspa5, ddit3, and hyou1 were upstream hub genes that control many downstream genes, such as mapk14, jun, ankrd1, and fgf19 (Supplementary Fig. 4B). These genes are associated with the production of ROS, unfolded protein, and ER-mediated stress (Supplementary Fig. 4B). To investigate mitochondrial gene expression associated with energy metabolism, we found a total of 25 DEGs linked to the electron transport chain (ETC). Among the genes encoding protein complex in ETC, cox6b2, cox7a2, nufa-3,-4,-5,-8, nduf-1,-2,-3,-6,-9 were significantly upregulated, while cox5a was only downregulated in BPA lineage females, suggesting dysregulation of mitochondrial energy metabolism pathway ( Fig. 4G and Supplementary Figs. 5 and 6). The genes encoding nuclear receptors were significantly downregulated except esr1 (Fig. 4G). Furthermore, oxidative stress-induced genes, mainly gpx7, gpx4, and gpx3 were highly upregulated in BPA lineage females (Fig. 3B). The present results indicated activation of oxidative stress in the liver of BPA lineage females, possibly mediated by mitochondrial and endoplasmic reticulum-mediated stress.
DEG-mediated pathways, Gene Ontology (GO), and IPA analysis revealed disease-related pathways triggered in females from the BPA lineage
In the BPA lineage males, upregulated genes in the DEG list were enriched in Th1 and Th2 cell differentiation, cysteine and methionine metabolism, the intestinal immune network for IgA production, and thiamine metabolism pathways (Fig. 5A), suggesting dysregulation of immunogenic response. As revealed by the KEGG pathway analysis, male DEGs were associated with metabolic pathways, protein processing in the endoplasmic reticulum, cholesterol metabolism, and PPAR signaling (Supplementary Fig. 7A). In the BPA lineage females, upregulated DEGs were typically associated with cell cycle and checkpoint, p53-independent DNA damage, unfolded protein response by activation of XBP1, respiratory electron transport, and apoptosis regulation pathways, suggesting cellular stress resulting in the development of severe phenotypic traits in females from the BPA lineage (Fig. 5B and Supplementary Fig. 7B).
Downregulated DEGs in the males from the BPA lineages were mainly linked to fatty acid metabolism, steroid biosynthesis, fatty acid elongation, and metabolic pathways (Fig. 5C). In contrast, beta-oxidation, PPARα regulation, lipid and lipoprotein metabolism, gluconeogenesis, and transport of fatty acid pathways were associated with downregulated DEGs, indicating perturbed lipolysis and fatty acid clearance mechanism in the liver of BPA lineage females (Fig. 5D). In total, 15 pathways common to both males and females, 2 pathways specific to males, and 198 pathways specific to females were found to be activated in the liver of BPA lineage fish, suggesting sex-specific associations of disease pathways in liver disease development (Supplementary Fig. 7C). Metabolic pathways, cholesterol metabolism, fatty acid elongation, and PPAR signaling pathways were found to be mutually triggered due to ancestral BPA exposure regardless of sex (Supplementary Fig. 7C). According to KEGG pathway analysis, NAFLD, cytokine-cytokine receptor interaction, leukocyte migration, antigen processing and presentation, Jak-STAT signaling and MAPK, AMPK, Insulin, thyroid hormone signaling pathways were found to be associated with up- and downregulated DEGs, respectively in the liver of BPA lineage females compared to control females (Supplementary Fig. 8).
Next, we performed gene ontology analysis with top-ranked significant DEGs in males and females from the BPA lineage. In the liver, Gene Ontology Enrichment (GEO) analysis revealed a total of 12 (male) and 94 (female) molecular functions (MF), 85 (male) and 678 (female) biological processes (BP), and 23 (male) and 183 (female) cellular components (CC) in BPA lineage fish (Supplementary Figs. 9 and 10). In the liver of BPA lineage males, cholesterol metabolic and biosynthetic process (GO:0016126 and GO:0008203) fatty acid synthase and elongase activity (GO:0004312 and GO:0009922) were significantly enriched. Cellular protein modification and metabolism (GO:0006464 and GO:0044267), ribosome biogenesis (GO:0042254), kinase activity (GO:0016301), mitochondrial membrane (GO:0005743 and GO:0031966), focal adhesion (GO:0005925) were detected in the liver of BPA lineage females. The results of IPA pointed out disease-specific canonical pathways, such as liver hyperplasia, hepatocellular carcinoma, liver steatosis, and liver necrosis in the liver of BPA lineage females (Supplementary Fig. 11). Adipogenesis, triacylglycerol degradation, and retinol biosynthesis were found in the liver of BPA lineage males (Supplementary Fig. 12), whereas liver necrosis, and hepatocellular carcinoma were found in the liver of BPA lineage females as per IPA. Upstream regulators hnf4a, tp53, and xbp1 expressed in the liver of the BPA lineage females (Supplementary Fig. 12).
Comparative analyses of biological processes altered in the liver due to direct vs indirect BPA exposure across humans, mice, and medaka
To investigate liver-specific pathways associated with direct and indirect BPA exposure, a comparison of KEGG pathway analyses was performed, as shown in the Venn diagram. The results showed that directly exposed mice and indirectly exposed medaka shared eight common pathways induced by BPA (Fig. 6A). In males, metabolic pathways, cholesterol metabolism, PPAR signaling pathway, cholesterol metabolism, lysosome, and carbon metabolism pathways were common in databases obtained from direct and indirect BPA exposure. In females, KEGG pathway analysis (Fig. 6B) found directly exposed mice and indirectly exposed medaka shared ninety-two common pathways induced by BPA. Common pathways in females represented oxidative phosphorylation, metabolic pathways, ribosomes, chemical carcinogenesis, and proteasomes in the liver. However, MTOR, AMPK signaling, insulin resistance, nonalcoholic fatty liver disease (NAFLD), and pathways in cancer were uniquely triggered in the female liver due to ancestral BPA exposure. Results indicated that in females, similar pathways were affected by the direct and indirect exposure to BPA. Compared to BPA lineage males and females, human and mouse transcription factors associated with steatosis were found most commonly in the liver of BPA lineage females (Supplementary Table 7), suggesting disease-specific transcription factors in humans and mice were significantly expressed in the female liver of BPA lineage.
Comparison of human NAFLD-NASH signature genes with those in female medaka from BPA lineage.
Ancestral exposure led to NAFLD-NASH phenotypes in the liver of BPA lineage females44. Here we screened for NAFLD-NASH signature genes in the livers of BPA lineage females and compared them with gene signatures associated with NAFLD-NASH in human patients. The BPA lineage females and human NAFLD patients shared 88 DEGs (Fig. 7A). Additionally, 20% of DEGs were mutually upregulated between medaka and human NAFLD patients (Fig. 7B). Among the total shared DEGs, 27.69% of DEGs were mutually downregulated. The total mutual DEGs triggered in the female liver of BPA lineage and human NAFLD patients were illustrated in Fig. 7C. A total of 52.31% of DEGs expressed in BPA lineage females were expressed in human patients (Fig. 7B). The commonly expressed DEGs in medaka and human NAFLD patients were found to be associated with immune response, cholesterol and lipid metabolic process (Fig. 7D). The common downregulated genes were abcb11, aldh6a1, igf1, cyp1a1, and apof. Among the shared upregulated genes, atf3, cxcr4, plin1, gins1, and alpk2 were common in BPA-lineage females and human patients (Supplementary Table 8). In BPA lineage females, gene signaling networks identified hsp90ab1, myc, ctnb1, tp53, and mitf as potential upstream genes involved in human NAFLD (Supplementary Fig. 13). In total, 14 pathways associated with NAFLD pathogenesis were found to be common between BPA lineage females and human NAFLD patients (Fig. 7E). Among the common pathways, PPAR signaling and metabolic pathways were highly significant in human NAFLD patients (Fig. 7F) and BPA lineage females (Fig. 7G), respectively. Genes associated with gluconeogenesis, lipid metabolism, mitochondria, ceramide metabolism, insulin pathway, inflammatory response, cell adhesion, coagulation, cytokines controlling liver fat, and advanced liver disease in humans were comparable to liver genes of BPA lineage females suggesting an association of similar DEGs in NAFLD pathogenesis (Supplementary Table 9). In addition, 39 common DEGs were detected in patients with NASH and BPA lineage females, and only 15 DEGs showed ancestral BPA-specific expression (Fig. 7H). Among common NASH-specific DEGs, igf1 was found to be an upstream regulator connected to downstream genes, including hnf1a, cyp7a1, mmp13, fbn1, and linked to adipogenesis in the liver (Supplementary Fig. 14). The synthesis of extremely long-chain fatty acyl CoAs, fatty acid triacylglycerol and ketone body metabolism, bile acid, and salt metabolism were found to be linked to the common NASH genes network (Supplementary Fig. 15).