Abstract

In recent years, there has been a growing emphasis on research investigating the transmission across generations of the effects of exposure to environmental factors, which may predispose to chronic diseases. It has been hypothesized that the propagation of such effects is mediated by alterations in gene regulatory elements, such as DNA methylation, histone modifications, and non-coding RNAs. Studies in Drosophila and mice have demonstrated that the compartmentalization of eukaryotic genomes into heterochromatin and euchromatin can mediate multigenerational metabolism-disrupting effects elicited by exposure to various metabolism disruptors. These findings suggest that eukaryotic nuclear genomes may possess the capacity to integrate the impact of environmental cues in a metastable manner that is phenotypically relevant. Here, we present the results of a murine model to assess whether preconception exposure to three metabolism disruptors of distinct nature, including dietary factors and environmental toxicants, intended to emulate the complexity of human exposures, results in metabolic alterations in the offspring of exposed individuals. Our findings align with our central hypothesis but also open an unexpected avenue to explore whether preconception exposure to metabolism disruptors predisposes the offspring of exposed individuals to not only typical metabolic diseases such as obesity but also to complex metabolic-psychiatric conditions such as anorexia.

Introduction

Over the past two decades, substantial evidence has emerged demonstrating the significant role of environmental agents, including chemical, biological, physical, or psychosocial factors, in determining the prevalence of chronic noncommunicable diseases within human populations (1). Obesity is an example of a chronic metabolic condition whose prevalence has dramatically increased in the last 50 years, and only ~2.7% of cases can be attributed to genetic factors (2). Exposure to metabolism disruptors, including pesticides, stress, and microorganisms, increases the susceptibility to developing chronic metabolic conditions, such as obesity, which can be propagated across multiple generations (3). These exposures can alter lipid metabolism, energy balance, appetite regulation, and glucose homeostasis by interfering with hormone signaling pathways, including those involving insulin, leptin, and nuclear receptors (4).

Animal model research has been instrumental in assessing the significance of the impact of environmental exposures on health outcomes for exposed individuals and their subsequent generations, as well as determining their underlying mechanisms (5). Beginning in 2013, we demonstrated that exposure to the biocide tributyltin (TBT) in female mice one week prior to conception and during the in utero development and lactation of their offspring resulted in sexually dimorphic metabolic alterations in at least four offspring generations (F1, F2, F3, and F4) (6–9). To elucidate the underlying mechanism mediating the transgenerational metabolic disruption induced by TBT, we devised a series of integrative analyses of several traits related to the epigenomic function (7, 8). We compared epigenetic marks such as DNA methylation at a genomic scale and other structural and genomic traits known to be regulated by epigenetic modifications, including chromatin accessibility and gene expression, in the somatic tissues and mature male gametes of the F3 and F4 descendants of F0 females that had been exposed to TBT and controls (7, 8). In most of our analyses, we detected a dichotomy for TBT-dependent alterations with regard to the heterochromatic and euchromatic compartments of the mouse genome; changes in a particular trait were significantly biased in one direction in heterochromatin and the opposite direction in euchromatin (7, 8). In conjunction with our investigation of the distribution within the mouse genome of genes regulating metabolism and chromatin organization, we proposed that the most plausible mechanism underlying the TBT-dependent multigenerational effects is that TBT triggered a perturbation of chromatin organization with the ability to self-reconstruct through development and across generations and predispose the progeny of exposed females to metabolic disorders (7, 8). Although our analyses did not directly address which layers of chromatin organization were altered upon TBT exposure, our findings were consistent with the possibility that such perturbation encompassed an alteration of the heterochromatin/euchromatin compartmentalization (HEC) (7, 8).

To the best of our knowledge, the ability of HEC to function as an epigenetic mechanism of genome regulation capable of mediating multigenerational metabolic disruptions has only been demonstrated in a single study on Drosophila melanogaster. In 2014, Öst et al. observed a similar dichotomy in gene expression alterations for heterochromatic and euchromatic genes in the offspring of male Drosophila melanogaster subjected to diets with varying sugar content (10). Exposure to these diets also led to alterations in metabolic traits and the expression of heterochromatin markers (10). The association between the dichotomy of perturbations in epigenomic function traits and metabolic disorders in the descendants of individuals exposed to environmental modulators of metabolism in two independent studies could suggest that eukaryotic HEC is generally susceptible to modulation by environmental exposure to metabolism disruptors in a metastable manner, resulting in metabolic effects in the descendants of exposed individuals (7, 8, 10). Notably, these two studies differ in the nature, duration, and sex of exposure, the number of generations separating exposed and descendants, and even the epigenetic characteristics of the species under study. Drosophila melanogaster has very low levels of DNA methylation when compared to other species such as Mus musculus (11).

In this study, we aimed to determine whether preconception exposure of female mice to three distinct metabolism disruptors alters the metabolism of their F1 offspring through chromatin organization modifications, particularly in HEC. Historically, research on the effects of multigenerational effects resulting from female environmental exposures has predominantly focused on exposure paradigms that encompassed the pre- and perinatal development of the offspring of exposed females (12). Consequently, there has been limited attention given to the effects of female exposures that occurred solely before conception (12). However, existing knowledge regarding the developmental establishment of HEC and the relevance of human preconception exposures that elicit multigenerational chronic metabolic diseases suggests that the preconception exposure window can be a relevant window of susceptibility for environmental exposures disrupting the early establishment of HEC in the offspring of exposed individuals.

On the one hand, research conducted across multiple eukaryotic species employing diverse methodologies has demonstrated that the earliest indications of genome compartmentalization and heterochromatin formation emerge very early in metazoan development, predating the full activation of the zygotic genome that characterizes the maternal-to-zygotic transition (MZT) (13, 14). Furthermore, the study of reporters of heterochromatic function in Drosophila melanogaster has suggested that adult heterochromatin structures are established very early in development and are faithfully propagated throughout development (15). Recent theoretical modeling for the propagation of heterochromatin-related epigenetic marks indicates that the spatial compartmentalization of heterochromatin and the limited availability of enzymes essential for its maintenance are necessary for the faithful propagation through cell differentiation and development of heterochromatin-related structures (16). Given that during the MZT transition, the zygote is largely transcriptionally silent (17), all processes occurring along MZT, including the establishment of heterochromatin, occur at the expense of the limited amounts of material deposited in the oocyte and to a lesser extent in the sperm. Considering these factors, it appears plausible that any environmental exposure resulting in alterations of the gamete material required for the establishment of the zygotic heterochromatin subsequently led to perturbations of heterochromatin that are propagated throughout development, thereby mediating phenotypes that reflect the biased localization of genes with specific functionalities in relation to heterochromatin (8, 10).

On the other hand, although human multigenerational studies are limited, the few that exist support the hypothesis that exposure to environmental stressors before conception contributes to chronic conditions in their offspring of at least two generations. The Dutch Famine stands as a pivotal human study that elucidates the significance of identifying critical windows of susceptibility to environmental stressors and chronic disease prevalence in subsequent generations (18). Numerous analyses of this cohort have demonstrated that low calorie intake of women during early, mid, and late gestation can result in chronic metabolic, psychological, and immunological disorders in their progeny (18). Notably, the descendants of the early exposure group exhibited a substantial increase in obesity, coronary heart disease, and elevated blood cholesterol levels compared to those whose mothers were exposed to famine during mid and late gestation (18). Given that the early exposure to famine spanned 10 weeks before and 10 weeks after conception, the higher disease incidence in the descendants of this group could support the hypothesis that the preconception window is particularly susceptible to the action of environmental exposures with effects in the next generation. Other human studies have shown that alterations in food availability or tobacco use during puberty have an effect on their grand offspring (19, 20). Altogether, these findings underscore the profound impact of preconception and early gestation events on transgenerational health and disease outcomes.

The selection of the metabolism disruptors for our study was based on the recognition that humans are subjected to a multitude of environmental stressors of varying nature from the moment of conception until death. This concept is encapsulated within the exposome framework (21), which encompasses the intricate interplay between life experiences and exposures that contribute to modulating susceptibility to disease. We selected three metabolism disruptors, TBT, inorganic arsenic (iAs), and diet, to model real-life human exposures of very different nature. TBT is an example of a human-made metabolism-disrupting chemical commonly found in soil, dust, and water and in human samples from liver, placenta, and blood, that has been associated with multigenerational obesity (6–9, 22). iAs is a naturally occurring element ubiquitously found in soils, sediment, and groundwater, and it poses major threats to global public health as around 200 million people worldwide are exposed to high concentrations (23). Prenatal exposure to arsenic has been associated with increased adiposity and increased risk of type-2 diabetes in human populations (24). Maternal preconception exposure to environmentally relevant doses of iAs in mice led to transgenerational metabolic perturbations, but little is known about potential epigenetic alterations associated with such transgenerational phenotypes (25). Lastly, total Western diet (TWD), whose micro- and macro-nutrient content represents the diet 50% of the U.S. population takes, is an example of a lifestyle metabolism disruptor (26). Although the separate exposure to these three metabolism disruptors will not conveniently model the complexity of the exposome, it can serve as a stepping stone in our understanding of whether environmental factors humans are exposed to and that have been associated with metabolic disorders do so through similar mechanisms of action.

MethodsMouse work

We conducted all mouse procedures at the University of California, Santa Cruz (UCSC) Vivarium. These procedures were reviewed and approved by the UCSC Institutional Animal Care and Use Committee (UCSC IACUC) as part of the animal protocol Chamr1908.

We purchased 160 female and 80 male C57BL/6J mice from Jackson Laboratory (strain #000664). To acclimate mice to the UCSC Vivarium, we scheduled their arrival a week before the start of the experiment. Upon arrival, we ear-punched and weighed the female mice and housed them in four-mouse cages. To minimize initial body weight disparities between the experimental groups, we randomly assigned the cages to five groups. We ranked the cages based on the cumulative body weights of the mice within each cage, from highest to lowest. Then, we divided the ranked list into eight subgroups of five cages each and randomly assigned each cage within each subgroup to one of the experimental groups.

We exposed 5-week-old female mice in each group to the following combinations of treated drinking water and diet for 3.5 weeks:

- DMSO (negative control) group: 0.1% dimethylsulfoxide (DMSO) in water and control diet (CD; Envigo Teklad Diets TD.140148).

- 5TBT group: 5 nM tributyltin (TBT) and 0.1% DMSO in water and CD.

- 50TBT group: 50 nM tributyltin and 0.1% DMSO in water and CD.

- IAS group: 10 µg/L sodium meta(arsenite) (inorganic arsenic, iAs) and 0.1% DMSO in water and CD.

- TWD group: 0.1% DMSO in water and Total Western Diet (TWD; Envigo Teklad Diets TD.110919).

The Agency for Toxic Substances and Disease Registry (ATSDR) determined that the no-observed adverse effect level (NOAEL) for TBT is 0.025 mg/Kg/day (27). 50 nM TBT in drinking water is equivalent to 0.005 mg/Kg/day, assuming an average mouse weight of 30 g and an average of 10 mL of water taken per day. Thus, 5 and 50 nM TBT represent 50 and 5 times lower than the NOAEL, respectively. We used 10 µg/L of inorganic arsenic in the drinking water, as that is the allowable level established by the U.S. Environmental Protection Agency (28). Low-fat CD and high-fat TWD have energy densities of 3.8 and 4.4 Kcal/g, respectively, and 17.2 and 34.5% of their energy content is provided by fat, respectively. TWD represents the diet 50% of U.S. people eat on a regular basis (26).

Immediately before starting the experiment and weekly thereafter, we weighed the mice before and after a four-hour fast. To protect female mice from any adverse effects of fasting on reproduction, we did not fast them at the end of the exposure period and before mating. We replenished the mouse feed weekly and water bottles semiweekly to maintain freshness and control consumption.

After 3.5 weeks of exposure at 8 weeks of age, we mated exposed female mice with unexposed same-age male mice for a week to produce F1 offspring. We checked females daily for successful matings by observing copulatory plugs. Mice with plugs were returned to their same cage mates from the exposure period for the rest of the pregnancy. Two days before birth, we relocated pregnant females to individual cages with ample bedding for nest-building. We monitored births and litter welfare using minimally invasive methods until the pups were old enough for further manipulations.

At 8–11 days of age, we toeclipped F1 mice. At 3 weeks of age, we weaned and selected at least 10 F1 mice per sex and exposure group. We housed mice of the same sex and exposure group together. We weighed F1 mice at weaning and weekly thereafter until the end of the experiment. We fed F1 mice CD between 3 and 7 weeks of age and TWD between 7 weeks of age and the end of the experiment. Before changing diets at 7 weeks of age, we weighed F1 mice before and after fasting for 12 hours. We measured fasting glucose using a Contour® blood glucose meter and strips (Ascensia Diabetes Care) with blood drawn from the mouse tails. At 12 weeks of age, we fasted F1 mice for 12 hours and euthanized them using an overdose of isoflurane followed by cardiac exsanguination. We determined fasting glucose levels from peripheral blood as previously indicated and harvested gonadal and inguinal white adipose tissue (gWAT and iWAT, respectively) depots and liver. We drew blood from the heart using EDTA-flushed syringes to minimize coagulation and collected cardiac blood samples in microcentrifuge tubes with 6 µL of a 100X protease inhibitor cocktail (Sigma Aldrich NC2042678). We centrifuged blood samples at 5,000 x g for 10 minutes at 4°C to separate plasma. We weighed gWAT, iWAT, and liver samples and snap-frozen them in dry ice. We maintained plasma, gWAT, iWAT, and liver samples frozen for downstream analyses. Dissections were performed by trained Chamorro-Garcia group members. To minimize daily workload, we divided F1 mice into five groups with equal representation of each exposure group and sex. Each group underwent dissection over five consecutive days. The lead experimenter (Diaz-Castillo) randomly assigned mice to dissectors who were unaware of the group or sex of the mice.

Metabolic analyses

We determined body weight by placing each mouse in an empty beaker atop a standard laboratory balance (Supplementary Data 6, 7, 12, 13, 15). To minimize age-related differences, especially in younger F1 mice, we corrected body weight data by the days of age of each mouse. We assessed the significance of body weight differences between exposure groups (5TBT, 50TBT, IAS, and TWD) and the control group (DMSO) using unmatched-measures Monte Carlo-Wilcoxon (uMCW) tests (see Statistics section) (Supplementary Data 21). We assessed the significance of body weight changes associated with fasting challenges between exposure and control groups using matched-bivariate Monte Carlo-Wilcoxon (mbMCW) tests (see Statistics section) (Supplementary Data 22).

To determine water consumption, we filled water bottles in each cage with 200 mL of treated water (entry water) and measured the volume of the remaining water after three or four days (exit water) (Supplementary Data 8). To determine food consumption, we weighed the food provided in each cage (entry food) and the remaining food after one week (exit food) (Supplementary Data 9, 10). We used mbMCW tests to assess the significance of differences in entry and exit food and water between exposure and control groups (Supplementary Data 22).

We corrected plasma glucose levels of fasted mice by body weight and the days of age of each mouse (Supplementary Data 14, 15). We used uMCW tests to assess the significance of fasting glucose differences between exposure and control groups (Supplementary Data 21).

We weighed gWAT, iWAT, and liver samples harvested from F1 mice using a precision balance (Supplementary Data 15). We corrected tissue weights by body weight and the days of age of each mouse. We used uMCW tests to assess the significance of tissue weight differences between exposure and control groups (Supplementary Data 22).

We submitted 50 µL of F1 plasma samples from euthanized mice to Eve Technologies (Calgary, Canada) for analysis of 11 metabolites in the Mouse, Rat Metabolic Array (MRDMET): amylin, gastric inhibitory peptide (GIP), ghrelin, glucagon-like peptide 1 (GLP-1), insulin, leptin, peptide YY (PYY), glucagon, pancreatic peptide (PP), resistin, and connecting peptide (C-Peptide) (Supplementary Data 16). We set to zero any data outside the standard curve range or that could not be mathematically extrapolated. We corrected plasma levels of each metabolite by body weight and the days of age of each mouse. We used uMCW tests to assess the significance of differences in plasma levels between exposure and control groups (Supplementary Data 21). We used the prcomp function from the R package Stats (version 4.5) (29) to perform Principal Component Analysis (PCA) on metabolites with valid data for all female or male mice to determine the similarities between experimental groups. We used the fviz_eig and fviz_pca functions from the R package Factoextra (version 1.0) (30) to visualize PCA results.

Transcriptomic analyses

We extracted RNA from F1 gWAT and liver using a VWR® Micro Homogenizer (catalog number #10032-328) and Qiagen RNeasy Plus kits (catalog number #74034), following the manufacturer’s protocols. We submitted the RNA samples to the University of California, Irvine (UCI) Genomics High Throughput Facility for sequencing using Illumina TruSeq Stranded mRNA kits and an Illumina NovaSeq 6000 sequencer. Libraries obtained from 100 RNA samples (2 sexes x 2 tissues x 5 experimental groups x 5 replicates per group) were sequenced up to three times to yield approximately 30 million 150-bp paired reads per sample. We used the following accessory files and three sets of informatics tools to analyze the sequencing results.

Accessory files for RNA-seq analyses

We downloaded FASTA files containing sequences of all major autosomes and the X chromosome, and the file “mm39.ncbiRefSeq.gtf.gz” with primary gene transcript annotation for the June 2020 GRCm39/mm39 mouse genome assembly from the UCSC Genome Browser (https://hgdownload.soe.ucsc.edu/downloads.html#mouse) (31). We retrieved the most recent core Gene Ontology from the Gene Ontology website (https://geneontology.org/) (32, 33).

Linux/Python tools

We used isoSegmenter (https://github.com/bunop/isoSegmenter) (34) to define isochores and isochore classes for the mouse genome. We executed isoSegmenter with the FASTA files for individual chromosomes downloaded from the UCSC Genome Browser. isoSegmenter operates in three phases: segmenting the provided sequence into non-overlapping 100 kb windows, assigning windows to five classes based on GC content, and defining isochores by concatenating juxtaposed windows of the same class. The five isochore classes and their respective GC contents are: L1 (below 37%), L2 (37%-41%), H1 (41%-46%), H2 (46%-53%), and H3 (exceeding 53%).

Galaxy Platform tools

We uploaded the FASTQ files obtained from the UCI Genomics High Throughput Facility to the Galaxy Platform website (usegalaxy.org) (35). We used FastQC (version 0.12) to assess read quality, multiQC (version 1.27) (36) to generate multi-sample reports, Cutadapt (version 5.0) (37) to remove adapter sequences, RNA Star (version 2.7) (38) to map reads to mouse genes, featureCounts (version 2.0) (39) to count reads, and Intersect (version 1.0) to determine the overlap of genes and isochores in the mouse genome. We used the annotation file “mm39.ncbiRefSeq.gtf.gz” to map reads to splice junctions. After completing all Galaxy Platform operations, we downloaded a table with RNA-seq gene-wise read counts for each gWAT and liver sample and isochore class associations (Supplementary Data 17).

R tools

We conducted all transcriptomic analyses using R (version 4.5) (29) packages in RStudio (Version 2024.12.1 + 563) (40). We used base R functions and the data.table package (version 1.17) (41) for data operations unless stated otherwise.

We used the R package ComplexHeatmap (version 2.24) (42) to visualize expression patterns across samples of the 100,000 most expressed genes. We aggregated raw counts for each gene across all samples and identified genes with the top 100,000 cumulative counts. We normalized gene read counts to counts per million (cpm) by dividing raw read counts for each gene by the sum of raw read counts for all genes in each sample and multiplying the result by 1 million. To rescale normalized counts to the same 0–1 range, we used the formula: (x - min(x))/(max(x) - min(x)), where x represents cpm values for each gene and sample, and min(x) and max(x) denote the minimum and maximum cpm values for each gene across samples, respectively.

We used uMCW tests to assess the significance of gene expression differences between exposure and control groups (see Statistics section) (Supplementary Data 23). We extracted read count tables for each of the 16 exposure versus control contrast combinations (2 sexes x 2 tissues x 4 exposure groups). We considered a gene to be expressed if at least 2 sequencing reads mapped it in at least 2 samples in each contrast and discarded non-expressed genes. We normalized gene read counts to cpm as previously indicated.

To determine the association of genes with extreme expression differences between exposure and control groups with specific functionalities, we performed pre-ranked Gene Set Enrichment Analyses (GSEA) (43) using the R package fgsea (version 1.34) (44) and the most updated mouse Gene Ontology (GO) annotation obtained from the R package msigdbr (version 10.0) (45). We restricted our analyses to GO sets of the Biological Process (BP) category with at least 15 and less than 500 genes (Supplementary Data 24). For each contrast, we ranked genes using the uMCW bias indexes (BIs) from highest to lowest, with the top genes being more expressed in the exposure group, and the bottom genes being more expressed in controls. GSEA calculates an enrichment score (ES) measuring the overrepresentation of a GO-BP set at either extreme of the pre-ranked list of genes, a normalized enrichment score (NES) to account for set size, and a P value calculated by comparing observed ES with a null distribution of ES obtained by permuting gene ranks (43).

To assess the significance of concerted changes in gene expression across the entire transcriptome and specific genomic regions such as chromosomes, individual isochores, and isochore classes, we performed biased-measures Monte Carlo-Wilcoxon (bMCW) tests (see Statistics section) (Supplementary Data 25, 27). We restricted these analyses to genes on the main chromosomes of the nucleus for both sexes by filtering out genes in the mitochondrial genome, Y chromosome, and unassembled chromosome segments. We also filtered out non-expressed genes and normalized read counts for each sample as previously indicated.

We determined the relative fraction of genes overlapping isochores of each class for genes in each autosome and the X chromosome, and for genes associated with specific GO-BP terms, for each exposure versus control group contrast (Supplementary Data 26). For each chromosome, GO-BP term, and isochore class, we calculated relative gene fractions as log10((x/n)/(X/N)), where x and X represent the number of genes in each chromosome or associated with each GO-BP term overlapping isochores of each class or any class, respectively, and n and N represent the number of genes for the whole transcriptome overlapping isochores of each class or any class, respectively.

We used the R package smplot2 (version 0.2) (46) to visualize the outcomes of Pearson correlation analyses between leptin plasma levels and Lep gene expression in gWAT, and between Lep and Lnc-Lep gene expressions in gWAT.

Other informatic tools

We used R packages data.table (version 1.17) (41), ggplot2 (version 3.5) (47), ggplotify (version 0.1) (48), ggrepel (version 0.9) (49), ggtext (version 0.1) (50), patchwork (version 1.3) (51), RColorBrewer (version 1.1) (52), scales (version 1.3) (53), and svglite (version 2.1) (54). The Supplementary Code file includes R code to reproduce our analyses and figures.

Statistics

To facilitate the integrative assessment of statistical significance in the differences between exposure and control groups for litter, metabolic, and transcriptomic traits, we conducted most of our analyses using Monte Carlo-Wilcoxon (MCW) tests. Initially, we devised MCW tests to ascertain whether a set of matched-paired quantitative measures exhibited a statistically significant bias in the same direction when compared with what would be expected by chance (7, 8, 55–57). Recently, we have reformulated MCW tests to interrogate four distinct data structures and developed the R package MCWtests (58) to publicly distribute the functions that facilitate conducting such tests. A comprehensive description of MCW testing rationale and specific MCW tests can be found at https://diazcastillo.github.io/MCWtests/index.html. Briefly, all MCW tests follow two basic steps. Firstly, to quantitatively determine the magnitude and direction of the bias between two conditions for the measure being analyzed, MCW tests compute a bias index (BI) that spans the range of 1 to -1, indicating that the measure is completely biased in each conceivable direction. Subsequently, to ascertain the statistical significance of the BIs computed for the user-provided dataset (observed BIs), a series of expected-by-chance BIs are generated by repeatedly rearranging the original dataset and computing BIs for each iteration. Pupper and Plower values are subsequently calculated as the proportions of expected-by-chance BIs that exhibit values equal to or higher than and equal to or less than the observed BIs, respectively.

MCW tests are particularly suitable for highly integrative studies like this one. Like other non-parametric approaches, MCW tests do not require original data to follow a specific distribution or undergo mathematical transformations. Also, since BI values always range from -1 to 1, comparing results from MCW tests conducted with data from different ranges, scales, or even traits becomes straightforward.

In our study, we used three of the four MCW tests offered by the R package MCWtests: unmatched-measures MCW (uMCW), matched-measures bivariate MCW (mbMCW), and biased-measures MCW (bMCW) tests.

uMCW tests assess whether two sets (e.g., a and b) of unmatched measures exhibit significant bias in the same direction. The uMCW testing process involves the following steps. First, all possible disjoint data pairs using measures from both sets are drawn. Second, the second measure in each pair is subtracted from the first measure. Third, measure pairs are ranked based on the absolute value of their differences, with the lowest absolute difference receiving the lowest rank. Measure pairs whose subtraction equals 0 are assigned a 0 rank. Measure pairs with the same difference values are assigned the lowest rank. Fourth, ranks are assigned a sign based on the sign of the measure differences calculated in the third step. Fifth, signed ranks for each subtraction type (e.g., a-b and b-a) are aggregated. Sixth, bias indexes (BIs) are calculated as the sum of signed ranks divided by the maximum value this sum would have if all measures in the first set (e.g., a) were greater than those in the second set (e.g., b). Finally, the significance of observed BIs is determined by comparing them with a collection of expected-by-chance BIs calculated after repeatedly rearranging the measure assortment between the two sets, and calculating Pupper and Plower values as previously indicated. uMCW tests can follow two paths based on the user-defined parameter max_rearrangements. If the number of potentially distinctive data rearrangements is less than max_rearrangements, uMCW tests will draw all potentially distinctive rearrangements, and Pupper and Plower values will represent an exact estimation of the significance of observed BIs. If the number of potentially distinctive data rearrangements is greater than max_rearrangements, uMCW tests will perform a number of rearrangements equal to max_rearrangements, and Pupper and Plower values will represent an approximated estimation of the significance of observed BIs.

mbMCW tests assess whether two sets (e.g., x and y) of inherently matched-paired measures are significantly differentially biased in the same direction. The mbMCW testing process involves the following steps. First, for each matched-pair of measures, the two potential subtractions are calculated (e.g., a-b and b-a). Second, for each subtraction type (e.g., a-b and b-a), measure pairs are ranked, with the lowest absolute difference receiving the lowest rank. Measure pairs whose subtraction equals 0 are assigned a 0 rank. Measure pairs with the same difference values are assigned the lowest rank. Third, ranks are assigned a sign using the sign of the measure differences calculated in the first step. Fourth, signed ranks for each subtraction type (e.g., a-b and b-a) and paired measure set (e.g., x and y) are aggregated. Fifth, bias indexes (BIs) are calculated as the sum of signed ranks divided by the maximum value that sum would have if all measure differences in the first set (e.g., x) were greater than those in the second set (e.g., y). Finally, the significance of observed BIs is determined by comparing them with a collection of expected-by-chance BIs calculated after repeatedly rearranging the measure matched-pair assortment between the two sets, and calculating Pupper and Plower values as previously indicated. mbMCW tests can also follow two paths to calculate exact and approximated estimations of the significance of observed BIs as indicated for uMCW tests.

bMCW tests are a combination of two tests that assess whether a set of measures of bias for a quantitative trait between two conditions and a specific subset of these bias measures are significantly biased in the same direction. The bMCW testing process involves the following steps. First, all bias measures are ranked using their absolute values, with the lowest measure receiving the lowest rank. Second, ranks are assigned a sign using the sign of the bias measure under study. Third, a whole-set bias index (wBI) is calculated by summing signed ranks for all elements in the dataset and dividing it by the maximum number this sum could have if all bias measures were positive. Third, for each subset of elements in the dataset under study, a subset bias index (sBI) is calculated by summing signed ranks for the elements in the subset and dividing it by the maximum value this sum would have if the elements in the subset had the most positive values in the whole dataset. Fourth, the significance of observed wBIs is determined by comparing them with a collection of expected-by-chance wBIs calculated after repeatedly rearranging the signs of all signed ranks, and calculating Pupper and Plower values as previously indicated. Finally, the significance of observed sBIs is determined by comparing them with a collection of expected-by-chance sBIs calculated after repeatedly rearranging the number of elements in the whole dataset that are associated with the subset under study, and calculating Pupper and Plower values as previously indicated. bMCW tests can also follow two paths to calculate exact and approximated estimations of the significance of observed wBIs and sBIs as indicated for uMCW tests.

We used uMCW tests to determine whether mouse body weights, fasting glucose, tissue weights, plasma metabolite levels, litter size, sex ratio, and gene expression for individual genes were significantly biased in the same direction for each exposure group compared to controls (Supplementary Data 21, 23). We used mbMCW tests to determine the significance of fasting body weights and water and food consumption for each exposure group compared to controls. For fasting body weights, we conducted mbMCW tests for mouse pre-fasting and post-fasting body weights data pairs for each experimental group (Supplementary Data 22). For water and food consumption, we conducted mbMCW tests for mouse cage entry and exit, and water and food data pairs for each experimental group. We used bMCW tests to determine whether gene expression differences for each exposure versus control contrast for genes in the whole transcriptome, overlapping specific isochores, overlapping isochores of the same class or located within autosomes and the X chromosome were significantly biased in the same direction (Supplementary Data 25, 27). The measure for gene expression difference we used for bMCW tests was the uMCW test BIs obtained for each gene and exposure versus control contrast. For all MCW tests, we set max_rearrangements parameter to 10,000.

We used principal component analysis (PCA), gene set enrichment analysis (GSEA), and Pearson correlations, as previously described in the Metabolic and Transcriptomics section. We used the same threshold of significance of P = 0.05 for all analyses.

ResultsDirect exposures to TBT, inorganic arsenic, and the Total Western Diet cause different metabolic disruptions

We randomly divided 160 female C57BL/6J mice into five groups. Each group received one of the following combinations of treated drinking water and diet for 3.5 weeks, starting at 5 weeks of age: DMSO and CD (negative control), 5 nM TBT and CD (5TBT group), 50 nM TBT and CD (50TBT group), 10 µg/L sodium (meta)arsenite and CD (IAS group), and DMSO and TWD (TWD group) (Figure 1a). We repeatedly measured body weight, fasting body weight, water, and food consumption to assess the efficacy of these exposures and determine if 5TBT, 50TBT, IAS, and TWD exposures caused metabolic disruptions compared to DMSO controls (Figure 1).

To determine whether body weights, fasting body weights, water consumption, and food consumption were significantly different between 5TBT, 50TBT, IAS, and TWD groups when compared with the control DMSO group, we used Monte Carlo-Wilcoxon (MCW) tests (see Methods for more details). MCW tests compute a bias index (BI) spanning 1 to -1 to quantify and determine the direction of bias between two conditions. To assess the statistical significance of the BIs computed for the user-provided dataset (observed BIs), a series of expected-by-chance BIs are generated by repeatedly rearranging the dataset and computing BIs for each iteration. Pupper and Plower values are calculated as the proportions of expected-by-chance BIs that exhibit values equal to or higher than and equal to or less than the observed BIs, respectively. MCW tests include four different tests to interrogate different data structures. We used unmatched-measures MCW (uMCW) tests to determine whether mouse body weights were significantly biased in the same direction for each exposure group compared to controls (Supplementary Data 21). We used matched-measures bivariate MCW (mbMCW) tests to determine the significance of fasting body weights and water and food consumption for each exposure group compared to controls (Supplementary Data 22).

Although 5TBT, 50TBT, IAS, and TWD exposures appeared to disrupt the metabolism of exposed females, these disruptions manifested differently depending on the exposure. 5TBT exposure resulted in a significant reduction in body weight immediately after exposure, which was subsequently mitigated (Figure 1c). 50TBT exposure also resulted in a reduction in body weight, but such a difference was not significant at any timepoint (Figure 1c). No consistent patterns were observed for fasting body weight or water and food consumption (Figure 1d) for either TBT exposure. While TBT mice exhibited a tendency to consume more food and drink less water compared to DMSO controls (Figures 1e–g), water consumption was only significantly decreased for the 5TBT group in a single timepoint (Figure 1e).

Direct effects of exposing F0 female mice to TBT, inorganic arsenic, and the Total Western Diet. (a) Experimental timeline detailing primary mouse operations for exposure to TBT, iAs, and TWD in female mice, and the timepoint structure used to assess variations in body weight, fasting body weight, water, and food consumption between exposure groups and controls. (b) Body weight dynamics for each mouse within each experimental group, from entry at UCSC Vivarium to timepoint T3.5. (c) Comparison of body weight normalized by age (g/days old) between exposure and control groups using uMCW tests. (d) Comparison of fasting body weight normalized by age (g/days old) between exposure and control groups using mbMCW tests. (e) Comparison of water consumption (mL) between exposure and control groups using mbMCW tests. (f) Comparison of food mass consumption between exposure and control groups using mbMCW tests. (g) Comparison of food caloric consumption (kcal) between exposure and control groups using mbMCW tests. Number of replicates (N) for each experimental group is indicated in each panel.

IAS exposure also led to a substantial reduction in body weight immediately after exposure, which persisted throughout the entire exposure period (Figure 1c). Additionally, during the initial week of exposure, but not thereafter, fasting body weight in IAS females was significantly higher compared to DMSO controls (Figure 1d). This observation suggests that the reduction in body weight observed in IAS females may not have been equally responsive to the fasting challenge as the controls during the time when the body weight differences were most pronounced. Furthermore, IAS females appeared to consume more food and drink less than controls throughout the exposure period, although these differences were rarely significant (Figures 1e–g).

The group exhibiting the most pronounced and substantial differences was the TWD group. These female mice demonstrated significantly higher body weights, lower fasting body weights, increased water consumption, and decreased food consumption compared to controls throughout the exposure period (Figures 1c-g). The body weight excess in TWD mice is more readily mobilized even during a mild fasting regimen, which aligns with the substantial increase in water consumption observed in TWD mice when compared to controls (Figures 1c–e).

All these patterns are consistent with TBT and IAS exposures causing similar decreases in body weight, with the IAS group weight loss being more pronounced, and TWD exposure causing an increase in body weight. The fact that these four exposures lead to distinct metabolic disruptions provides a great opportunity to comparatively assess if they lead to distinctive metabolic disruptions in the offspring of exposed females.

Preconception exposures of female mice to TBT, inorganic arsenic, and the Total Western Diet induce sexually dimorphic metabolic disruption in their offspring

At 8 weeks of age and after