Paternal microbiome perturbations impact offspring fitness

Abstract

The gut microbiota operates at the interface of host–environment interactions to influence human homoeostasis and metabolic networks^1,2,3,4. Environmental factors that unbalance gut microbial ecosystems can therefore shape physiological and disease-associated responses across somatic tissues^5,6,7,8,9. However, the systemic impact of the gut microbiome on the germline—and consequently on the F₁ offspring it gives rise to—is unexplored¹⁰. Here we show that the gut microbiota act as a key interface between paternal preconception environment and intergenerational health in mice. Perturbations to the gut microbiota of prospective fathers increase the probability of their offspring presenting with low birth weight, severe growth restriction and premature mortality. Transmission of disease risk occurs via the germline and is provoked by pervasive gut microbiome perturbations, including non-absorbable antibiotics or osmotic laxatives, but is rescued by restoring the paternal microbiota before conception. This effect is linked with a dynamic response to induced dysbiosis in the male reproductive system, including impaired leptin signalling, altered testicular metabolite profiles and remapped small RNA payloads in sperm. As a result, dysbiotic fathers trigger an elevated risk of in utero placental insufficiency, revealing a placental origin of mammalian intergenerational effects. Our study defines a regulatory ‘gut–germline axis’ in males, which is sensitive to environmental exposures and programmes offspring fitness through impacting placenta function.

Main

Sperm transmit heritable information to the next generation in the form of genetic (DNA sequence) and epigenetic (non-DNA sequence-based) material^11,12. Evidence across phyla indicates that the epigenetic component, including chromatin states, small RNAs and macromolecules, has the potential to be modified by the preconception environment and influence offspring phenotype^{13,14,15,16,17,18,19,20,21}. Nonetheless, the extent and underlying mechanisms of paternally inherited epigenetic effects in mammals remain opaque²², while how environmental factors converge and signal to germ cells is also unclear. The gut microbiota is increasingly appreciated to play a principal role in integrating environmental signals into host responses^1,9. Indeed, the mammalian lifecycle has evolved in the presence of gut microbial communities, which have assumed an important role in metabolic, hormonal and immune function, implying that hosts exploit a healthy microbiome to maximize fitness. The composition of gut microbial ecosystems is nonetheless profoundly shaped by environmental factors, such as diet and medication^23,24,25, and thus loss of gut biodiversity associated with modernization poses a risk to human health. Yet, despite accumulating evidence that an imbalance (dysbiosis) of the microbiome triggers physiological responses across somatic tissues^{5,6,7,8,26,27,28,29}, little is known of the effect of microbiome perturbations on the germline.

To investigate this, we established an inducible model of gut microbiota dysbiosis in isogenic male mice using ad lib non-absorbable antibiotics (nABX) (Fig. 1a). These nABX cannot cross the gastrointestinal epithelium and thus any reproductive responses reflect acute perturbation of gut microbial communities rather than systemic drug exposure³⁰. Applying 16S ribosomal RNA sequencing showed that 6 weeks of low-dose nABX treatment (denoted 6 wk) led to a marked reduction in diversity, abundance and richness of gut microbiota (P = 0.000003, Wilcoxon rank sum), which is reversible and progressively recovers after 8 weeks of nABX withdrawal (Fig. 1b and Extended Data Fig. 1a,b). The dysbiosis after 6 weeks of nABX had no significant effects on male weight, fertility or survival (Extended Data Fig. 1c–h). Moreover, nABX were undetectable in circulating serum and testis using mass spectrometry, confirming their specific action in the gut (Extended Data Fig. 1i–l).

Fig. 1: Paternal gut dysbiosis probabilistically triggers major F₁ phenotypes.

a, Schematic showing the strategy for induced paternal dysbiosis and recovery using nABX. b, Quantification of microbial taxa richness in males after 6 weeks of nABX treatment and during the recovery (rec) time course, by 16S rRNA sequencing (CON t0 = 18, 6 wk = 18, 6 wk + 4 rec = 14, 6 wk + 8 rec = 11; nABX t0 = 19, 6 wk = 12, 6 wk + 4 rec = 13, 6 wk + 8 rec = 12 individuals per timepoint). Bar represents median, whiskers 1.5× interquartile range. c, Body weight of F₁ offspring at postnatal days P3 and P15 according to paternal nABX treatment. P value by two-tailed nested (hierarchical) t-test (CON n = 172, nested into N = 26 fathers; nABX n = 181 nested into N = 28 fathers). Bar indicates mean. d, Representative images of SGR phenotype in F₁ offspring from dysbiotic fathers (nABX-treated). e, Forest plot showing the log OR of risk for abnormal body-weight classes in offspring that survive to P15. Null effect is represented by a vertical line for which OR value is 1. Whiskers indicate 95% CI, P value by two-sided Chi-square test (mortality P = 0.0001; SGR P = 0.044). f, Kaplan–Meier plot showing postnatal survival of F₁ progeny depending on paternal nABX treatment regime (CON n = 179; nABX n = 199). P value by Mantel–Cox (log-rank) test. g, PCA of transcriptomes from F₁ brains derived from control or nABX sires. The nABX offspring are stratified by normal or SGR phenotype. SGR not observed in control offspring. h, Heatmap showing expression of top upregulated and downregulated genes in F₁ SGR brains, from independent litters sired by nABX-treated fathers. i, Left, OR of F₁ susceptibility for abnormal body weight and mortality when sired by males with dysbiosis induced by 6 weeks of treatment with general antibiotics (avaABX). Whiskers indicate 95% CI, P value by two-sided Chi-square test (mortality P = 0.014; SGR P = 0.038). Right, Kaplan–Meier plot showing postnatal survival of F₁ progeny from avaABX-treated males. j, Left, OR of F₁ risk for abnormal body weight when sired by males with dysbiosis induced by 6 weeks of bowel cleansing with PEG laxative. Whiskers indicate 95% CI, P value by two-sided Chi-square test (mortality P = 0.013; SGR P = 0.014). Right, Kaplan–Meier plot showing survival of F₁ progeny from PEG-treated males. n, offspring; N, litters; ND, not detected.

Source Data

Offspring of dysbiotic fathers

To assess the impact of induced microbiota dysbiosis on offspring, we mated nABX-treated males with naive (untreated) females and scored F₁ phenotypes, applying a rigorous nested statistical analysis for intergenerational significance. We found that offspring sired by nABX fathers had significantly lower neonatal birth weight (P3) relative to offspring from control fathers (P = 0.023, nested unpaired t-test; control (CON) n = 172 (26 litters), nABX n = 181 (28 litters)) (Fig. 1c). Both female (P = 0.017) and male (P = 0.029) offspring were affected, while litter size was constant (Extended Data Fig. 2a,b). Moreover, mean body weight of F₁ offspring sired by nABX males remained significantly lower throughout postnatal development (P = 0.013 at P15 (preweaning); P = 0.015 at P21 (weaning), nested unpaired t-test) (Fig. 1c and Extended Data Fig. 2c,d).

Among offspring fathered by dysbiotic males, we also observed a major but partially penetrant postnatal phenotype that manifested as severe growth restriction (SGR; body-weight Z-score < −3), which was not observed among control offspring (Fig. 1d and Extended Data Fig. 2e,f). Quantification of this showed a significantly increased odds-ratio (OR = 3.52; P = 0.044, Chi-square) of SGR amongst nABX-derived progeny (Fig. 1e). This was further reflected by a negative skew and excess kurtosis of nABX offspring body weights (skew −1.98; Rku 5.0). Most strikingly, however, F₁ offspring sired by nABX-treated males were associated with a highly significant increased rate of postnatal mortality (P = 0.0002, Mantel–Cox test) relative to offspring from control sires (Fig. 1f). This occurred preferentially amongst SGR offspring, suggesting that elevated mortality is linked with increased F₁ susceptibility to growth restriction.

Transcriptome profiling of SGR offspring fathered by independent nABX males showed that they clustered together by principal component analysis (PCA) and separately from offspring of control sires, indicating a reproducible F₁ molecular response (Fig. 1g). Indeed, 2,973 and 1,563 differentially expressed genes (DEGs) are detected in F₁ brain and brown adipose tissue (BAT), respectively. DEGs are preferentially enriched for reactome pathways relating to metabolic processes (top terms: metabolism P = 2.51⁻⁹; metabolism of lipids P = 0.000018) and are robust between SGR offspring from independent litters (Fig. 1h and Extended Data Fig. 2g,j). These data support an intergenerational impact of nABX-mediated paternal dysbiosis on offspring growth, metabolic networks and survival. Importantly, these phenotypes arise as probabilistic rather than deterministic responses to paternal status, and thus manifest as altered risk.

We next asked whether orthologous strategies of paternal microbiota perturbation also elicit an F₁ response. First, we used an alternative combination of antibiotics (avaABX) and observed that paternal ab lib administration of avaABX replicated the increased susceptibility of F₁ offspring to reduced body weight (SGR OR = 7.0; P = 0.038, Chi-square) and increased mortality (P = 0.014, Mantel–Cox test) (Fig. 1i). Second, we perturbed the paternal gut microbiota without any exposure to antimicrobial drugs, by performing gastrointestinal cleansing with osmotic laxatives (polyethylene glycol (PEG)), which induces widespread dysbiosis^31,32. We observed progeny fathered by PEG-treated males had significantly lower F₁ body weight (P = 0.021, nested unpaired t-test; CON n = 76 (13 litters), nABX n = 76 (13 litters)), with increased SGR susceptibility (OR = 5.8; P = 0.0142 Chi-square) and premature mortality (P = 0.013, Mantel–Cox test) (Fig. 1j). Thus, multiple distinct perturbations of the gut microbiota in prospective fathers increase the risk of offspring presenting with developmental impairment and premature mortality, supporting a direct link between induced paternal dysbiosis and offspring fitness.

Reversibility of paternal effects

Next, we investigated whether paternal recovery from gut dysbiosis can rescue the F₁ phenotypic effects. Following 4 weeks of nABX withdrawal (6 wk + 4 rec (recovery)), when a significant perturbation of the paternal microbiome persists, F₁ offspring again showed significantly lower neonatal weight (P = 0.012, nested unpaired t-test; CON n = 160 (24 litters), nABX n = 146 (22 litters)), increased susceptibility to SGR (OR = 8.1; P = 0.020 Chi-square) and impaired growth trajectory (Fig. 2a and Extended Data Fig. 3a,b). However, upon recovery of the paternal microbiome after 8 weeks of nABX withdrawal (6 wk + 8 rec) (Fig. 1b and Extended Data Fig. 1a,b), we observed a concurrent recovery of the F₁ neonatal weight phenotype (P = 0.55, nested unpaired t-test; CON n = 87 (13 litters), nABX n = 89 (13 litters)) and normal developmental growth (SGR OR = 0.97; P = 0.98, Chi-square) (Fig. 2b and Extended Data Fig. 3a,b). Moreover, progeny of dysbiotic 6 wk + 4 rec nABX fathers reproduced the significantly elevated F₁ mortality rate (P = 0.0009, Mantel–Cox test) (Fig. 2c) but progeny from the same fathers following microbiota recovery at 6 wk + 8 rec had no excess mortality (P = 0.73) (Fig. 2d,e). These data imply that the increased probability of F₁ effects persists during the period of paternal dysbiosis (6 wk, 6 wk + 4 rec) but is reversible and reverts coincident with recovery of the father’s gut microbiome (6 wk + 8 rec) and a spermatogenic cycle (about 5 weeks) (Fig. 2e).

Fig. 2: Recovery of paternal gut microbiota rescues susceptibility to F₁ phenotypes.

a,b, Growth curves of F₁ offspring sired by males during a time course of nABX withdrawal which still retain gut microbiota dysbiosis (6 wk + 4 rec) (a) and offspring of the same males after microbiota recovery (6 wk + 8 rec) (b), as well as aged-matched control sires. Line indicates mean, P value by nested (hierarchical) two-tailed t-test. c,d, Kaplan–Meier plot showing postnatal survival of F₁ progeny sired by dysbiotic (c) or recovered nABX males (d). P value by Mantel–Cox (log-rank) test. e, Forest plots showing OR of F₁ susceptibility to abnormal body weight and premature mortality when sired by dysbiotic (6 wk, 6 wk + 4 rec) or recovered males (6wk + 8 rec). Null effect is represented by a vertical line for which the OR value is 1. Whiskers indicate 95% CI, P value by two-sided Chi-square test (mortality: 6 wk P = 0.0001, 6 wk + 4 rec P = 0.0014, 6 wk + 8 rec P = 0.76; SGR: 6 wk P = 0.044, 6 wk + 4 rec P = 0.020, 6 wk + 8 rec P = 0.99). f, PCA of transcriptomes from F₁ SGR adipose tissue sired by 6 wk or 6 wk + 4 rec nABX or control males, each from several independent litters. g, Gut microbiota richness of offspring (left) is not affected by paternal gut microbial status and does not correlate with F₁ phenotype. Paternal microbiota richness (right) does correlate with F₁ offspring phenotype. Shaded areas indicate 95% CI. h,i, Left, neonatal F₁ body weight following IVF of isogenic oocytes using control or nABX-treated sperm donors. Right, F₁ body-weight distribution at P15 following IVF. Data are independently derived from CD1 strain surrogate mothers (CON n = 66 offspring, nested into N = 12 fathers; nABX n = 75, nested into N = 12) (h) or BL6 strain surrogate mothers (CON n = 33 offspring, nested into N = 8 fathers; nABX n = 41, N = 7 fathers) (i). P value by nested (hierarchical) two-tailed t-test.

Source Data

Of note, transcriptomics showed SGR offspring sired by 6 week nABX fathers clustered with independent offspring from 6 wk + 4 rec fathers and exhibited highly similar gene ontology enrichments (Fig. 2f and Extended Data Fig. 3c,d). This suggests that affected offspring conceived during the window of paternal dysbiosis funnel into a consistent molecular phenotype, implicating a common underlying aetiology. Moreover, de novo genome sequencing of SGR offspring showed no elevated mutational load, nor explanatory differences in structural variants, single-nucleotide polymorphisms (SNPs) or small insertions and deletions (INDELs) relative to controls (Extended Data Fig. 4a,b). We were further unable to detect transmitted F₂ effects (Extended Data Fig. 4c–f). These data establish that F₁ phenotypes arising from dysbiotic fathers are not due to inheritance of genetic differences and do not propagate beyond the first generation.

Mode of intergenerational transmission

To examine the modality of intergenerational inheritance, we first asked whether there is paternal transmission of a dysbiotic gut microbiome³³. We found that postpartum mothers exhibited no significant compositional changes in their microbiome (P = 0.77, alpha diversity Wilcoxon; paternal CON n = 8, nABX n = 11) and did not cluster by paternal exposure, judged by faecal 16S profiling (Extended Data Fig. 5a). We further detected no significant effects on the maternal microbiome in the days after mating with nABX males, nor differences in offspring or seminal microbiota, suggesting against transmission of altered paternal microbiomes to mother and/or offspring (Extended Data Fig. 5b–h). Indeed, offspring phenotype correlated with the microbiome of its father rather than its own microbiota (Fig. 2g), with some specific taxa showing preferential associations (Extended Data Fig. 5e). We also examined potential coprophagic effects of residual nABX during parental cohousing, noting a small but non-significant trend which resolved in days (Extended Data Fig. 5f–h). To empirically test whether this or other indirect (non-germline) paternal factors impact offspring phenotype, we performed a series of cohousing experiments whereby females were maintained with a control or nABX male in its environment but then mated with independent treatment-naive males. We found no difference in birth weight (P = 0.85; nested t-test), growth (P = 0.98) or survival (P = 0.35, Mantel–Cox) of offspring, irrespective of previous maternal exposure to dysbiotic nABX males (Extended Data Fig. 6a–d). We conclude that neither paternal transfer of altered microbiota nor indirect maternal responses underlie F₁ effects.

We therefore asked whether F₁ phenotypes are transmitted specifically through paternal gametes, by performing in vitro fertilization (IVF). Here, isogenic oocytes are fertilized with sperm from either nABX-treated or control males and then implanted into CD1 high-quality surrogate dams, precluding any parental contact (Extended Data Fig. 6e). We found IVF progeny derived from dysbiotic sperm donors had significantly reduced neonatal birth weight (P = 0.034, nested unpaired t-test; CON n = 65, nABX n = 80), impaired postnatal growth (P = 0.047) and elevated SGR incidence relative to controls (Fig. 2h and Extended Data Fig. 6f). Independent IVF using BL6 recipient dams, which are relatively poor surrogates, reproduced F₁ birth weight effects from nABX sperm donors with greater effect size (Fig. 2i and Extended Data Fig. 6g) (P = 0.050, nested unpaired t-test; CON n = 33, nABX n = 41). These data suggest that paternally induced F₁ phenotypes arise in independent in utero genetic backgrounds and are transmitted primarily through the gametes and copurifying molecules.

The gut–germline axis

Transmission through the germline prompted us to investigate the physiological changes in the father’s reproductive system induced by acute gut microbiome dysbiosis. We first observed that dysbiotic males from 6 weeks of nABX exposure had significantly smaller testes by mass than did controls (P = 0.001, unpaired t-test; CON n = 31, nABX n = 32), which correlated with lower sperm count (Fig. 3a and Extended Data Fig. 7a,b). Histological analysis showed architectural changes in a subset of seminiferous tubules, including vacuoles formed by partial loss of germ cells, which was not observed in control testes (Fig. 3b and Extended Data Fig. 7a). Indeed, nABX males exhibited a significantly increased number of abnormal testis tubules (P = 0.032, nested Mann–Whitney; CON mean 0.64%, nABX mean 3.84%) and reduced epithelial thickness (P = 0.016; nested unpaired t-test) (Fig. 3c,d and Extended Data Fig. 7b–d). These data indicate that testicular physiology is impacted by gut microbiota perturbation.

Fig. 3: Testicular responses to gut microbiota dysbiosis indicate a regulatory gut–germline axis.

a, Boxplot showing testis mass to body weight ratio after 6 weeks of nABX treatment in males. Bar represents median and whiskers indicate 5th–95th percentile. P value by unpaired two-tailed t-test (CON n = 31; nABX n = 32). b, Representative haematoxylin and eosin stained histological sections of testes from control and nABX males. Seminiferous tubules from dysbiotic males show incidence of vacuoles formed by loss of epithelium and absence of mitotic (spermatogonial) compartments. Asterisks indicate abnormal tubules. c, Quantification of abnormal testis tubules in control and dysbiotic males. P value by nested Mann–Whitney test (CON n = 54 sections, nested into N = 4 males; nABX n = 72, N = 5; P = 0.032). Bar represents median and whiskers indicate 5th–95th percentile. d, Quantification of epithelial thickness of seminiferous tubules. P value by nested t-test (CON n = 854 tubules, N = 4 males; nABX n = 1,061, N = 4; P = 0.016). Bar represents median and whiskers indicate 5th–95th percentile. e, PCA of untargeted metabolomics profiles from independent testis of control or dysbiotic males (6 wk, 6 wk + 4 rec) and after gut microbiome recovery (6 wk + 8 rec). f, Volcano plot highlighting differentially abundant metabolites in testes after 6 weeks of nABX and during recovery to restore the gut microbiome. P value by two-tailed t-test adjusted for multiple testing. g, MA ((M (log₂ ratio) and A (mean average)) plot showing gene expression changes in testes from independent (n = 5) nABX males. h, Quantitation of leptin hormone level in testes (left) and circulating plasma (right) of dysbiotic males after 6 weeks of nABX, by ELISA. Bar indicates mean with 95% CI whiskers. P value by unpaired two-tailed t-test (CON n = 9; nABX n = 9). i, Testes features in males that are leptin deficient (ob/ob) for 6 weeks and wild-type controls (WT n = 3; ob/ob n = 3). Left, haematoxylin and eosin stained tubule sections; stars indicate abnormal histo-architecture. Right, testis weight, P value by unpaired two-tailed t-test and bar indicates mean with 95% CI. j, Transcriptome PCA of blastocysts derived from leptin-deficient or control fathers. Data points represent single-embryos from duplicate independent IVF experiments, each with triplicate independent fathers and littermate fathers as controls. Scale bars, 100 μm (left panel) or 50 μm (right panel) (b); 50 μm (i). FC, fold change.

Source Data

To characterize the reproductive system response to dysbiosis at the molecular level, we first performed untargeted metabolomic profiling of father’s testes, annotating 3,803 features³⁴. Global PCA of annotated metabolites showed testes cluster according to gut microbiota status, with clear separation between nABX and control males at both 6 wk and 6 wk + 4 rec, which corresponds to the period of transmitted F₁ effects (Fig. 3e). By contrast, testicular metabolomic profiles at 6 wk + 8 rec are indistinguishable, indicating that metabolites dynamically recover coincident with restoration of the gut microbiota and with reversal of transmitted F₁ effects. We identified 68 significant differentially abundant metabolites in testes of dysbiotic males, which included the fatty acid, anandamide, which acts in the endocannabinoid pathway and the signalling lipid, sphingosine-1-phosphate (S1P) (Fig. 3f and Extended Data Fig. 8a). Differentially abundant metabolites are specifically enriched for sphingolipids and glycerophospholipids, as well as endocannabinoids, which have all been implicated in germ cell function³⁵ (Extended Data Fig. 8b). We further investigated transcriptomic profiles in testes of dysbiotic males, observing limited expression changes at both bulk and single-cell levels (Fig. 3g and Extended Data Figs. 8c–f and 9a–g). Gene set enrichment indicated that glycerophospholipids- and steroidogenesis-related genes were preferentially dysregulated, consistent with altered metabolomic profiles (Extended Data Fig. 8e–g). The most sensitive gene was Leptin, however (Fig. 3g and Extended Data Fig. 9h), which encodes a hormone produced primarily by adipocytes but also by germ cells, with key roles in energy homoeostasis and reproduction^36,37.

The cumulative evidence indicates a significant change in the testicular environment in response to gut microbiome perturbation, including altered metabolite profiles, physiology and hormones. This suggests the existence of a gut–germline axis in mammals which carries an important homoeostatic function, analogous to recent reports in Drosophila^38,39.

Among the paternal responses to nABX was strong dysregulation of Leptin (Fig. 3g). Validation with ELISA showed that nABX-mediated dysbiosis leads to significantly reduced levels of leptin hormone, both systemically in circulating blood (P = 0.0001, unpaired t-test; CON n = 3; nABX n = 3) and specifically in testes (P = 0.0003; CON n = 9; nABX n = 9) (Fig. 3h). To investigate the implications of this, we used 6-week-old ob/ob mice (Leptin-null), which models the timespan of leptin deficiency induced by 6 weeks of nABX exposure. Such young ob/ob males had reduced testes weight (P = 0.014) and tubule abnormalities interspersed with normal histo-architecture, similarly to dysbiotic 6-week nABX males (Fig. 3i and Extended Data Fig. 10a–c). We therefore assessed the intergenerational impact of short-term leptin deficiency. Although natural matings failed, IVF using 6-week ob/ob sperm fertilized wild-type oocytes as efficiently as control sperm, implying underlying germ cell viability (Extended Data Fig. 10d). Single-embryo transcriptomics of the resulting blastocysts showed that offspring sired by leptin-deficient fathers strongly cluster away from embryos sired by control fathers (Fig. 3j and Extended Data Fig. 10e,f). Indeed, we detected more than 500 F₁ high-confidence DEGs, with chromatin pathways preferentially deregulated (Extended Data Fig. 10g–i). Importantly, Leptin itself is not expressed during early development ruling out a haploinsufficient effect of F₁ heterozygosity (Extended Data Fig. 10g). Taken together, these data show that leptin is systemically dysregulated by induced microbiome dysbiosis and that direct perturbation of paternal Leptin before conception has an intergenerational legacy on offspring gene expression programmes, implicating leptin as an important signalling component in the gut–germline axis.

We next sought to understand the impact of the gut–germline axis on mature gametes by searching for molecular changes in sperm. We initially charted DNA methylation at base resolution. Independent methylomes (n = 5) of purified sperm from nABX males were highly comparable to controls, with no change in DNAme globally or at genomic features (Fig. 4a and Extended Data Fig. 11a,b). We identified only 21 differentially methylated regions (DMRs) genome-wide (logistic regression P(adjusted (adj)) < 0.05 and >20% absolute change, 50 CpG tiles) (Fig. 4a and Extended Data Fig. 11c,d), which typically overlapped epivariable CpG shore regions, implying that they may reflect natural variation. Genomic imprints were all correctly established in both control and nABX sperm (Extended Data Fig. 11e). We next assayed sperm-borne small RNA with high-quality libraries from independent sperm collections (n = 9 males, pooled into N = 3). We detected significant abundance changes in several microRNA, including miR-141 (P(adj) = 1.15 × 10⁻⁹) and miR-200a (P(adj) = 0.008), which act together to regulate epithelial–mesenchymal transition and placental development (Fig. 4b and Extended Data Fig. 11f,g)^40,41. We also observed changes in the abundance of 5′ transfer RNA fragments (tRF) and in particular upregulation of tRF-Gly-GCC in dysbiotic males (Fig. 4b and Extended Data Fig. 11f–h), which has been implicated in intergenerational effects^42,43. Quantitative differences in small RNA abundance were confirmed in independent samples by TaqMan quantitative PCR (Extended Data Fig. 11i). Overall, although DNA methylation is relatively stable, the composition of small RNAs in sperm is modified in response to nABX-mediated dysbiosis. Considered with the altered metabolite and hormonal profiles, this suggests that compound changes in macromolecule composition are transmitted to offspring.

Fig. 4: Paternal dysbiosis induces F₁ placental insufficiency.

a, Left, heatmap showing DNA methylation levels across genomic features in sperm from control or nABX-treated males (n = 5 methylomes per condition) by whole-genome bisulfite-seq. Right, scatter plot of genome-wide DNA methylation (50 CpG tiles), with differential tiles highlighted. b, Heatmaps showing differential abundance of selected miRNA (left) and tRNA fragments (right) in pooled purified sperm from independent (n = 9) control or nABX males. c, PCA of transcriptomes from embryonic (E13.5) brain and placenta according to paternal exposure. d, PCA of placenta transcriptomes at E18.5 from independent litters. e, Volcano plot showing DEGs in E18.5 placenta sired by 6-week nABX-treated males. P value adjusted for multiple testing. f, Expression of key genes for placental development in E18.5 placenta. Bar indicates mean and whiskers are s.d. Each data point is an independent placenta (CON n = 5 placenta (3 litters); nABX n = 5 placenta (3 litters)). g, Ratio of fetal mass to placental mass at E18.5, depending on preconception paternal condition. Bar indicates mean, P value by unpaired two-tailed t-test (CON n = 82, nABX n = 53). h, Expression of biomarkers of pre-eclampsia (PE) in E18.5 placenta derived from nABX fathers. Bar indicates mean and whiskers are s.d. P value by multiple-testing corrected DESeq2. CON n = 5 placenta (3 litters); nABX n = 5 placenta (3 litters). i, Representative placenta sired by control or dysbiotic males (nABX), stained for DAPI (blue) and VE-cadherin (red) to demarcate the labyrinth zone (LZ). Quantification below shows LZ as a ratio of total area. Bars indicate mean ± 95% CI. Each data point is an independent placenta tissue section (CON n = 4 placenta (4 litters); nABX n = 6 placenta (6 litters)). JZ, junctional zone. P value by unpaired two-tailed t-test. j, Levels of placental growth factor (PLGF) protein (left), a marker of PE and the sFLT1/PLGF ratio (right), in F₁ placenta depending on paternal regime. P value by two-tailed t-test. Bars indicate median. Each data point is an independent placenta (CON n = 12 placenta (7 litters); nABX n = 12 placenta (6 litters); avaABX n = 8 placenta (3 litters)). Asterisks indicate P values of 0.05 or less* and 0.0001 or less***. Scale bar, 100 μm. Prom, promoter.

Source Data

Placental responses to dysbiotic fathers

To understand the mechanism through which sperm influence offspring phenotype, we sought to identify the initial source of embryonic defects. At E13.5 (about mid-gestation) we found no DEGs in embryos sired by nABX-treated males relative to control and embryonic transcriptomes could not be distinguished by PCA (Fig. 4c). By contrast, transcriptomes of placenta at E13.5 from independent matings clustered strongly depending on the paternal nABX regime (Fig. 4c), with 538 DEGs, which included downregulation of several Prolactin genes (Extended Data Fig. 12a,b). To examine this further, we assayed mature placenta at E18.5 and observed their transcriptomes clustered according to paternal microbiome status, with 348 high-confidence DEGs (Fig. 4d and Extended Data Fig. 12c). Upregulated DEGs are enriched for steroid metabolism (P(adj) = 2.31 × 10⁻⁷) whereas downregulated DEGs are associated with glycolysis (P(adj) = 0.00039) (Fig. 4e). Notably, they also include downregulation of several factors important for placenta development, such as Hand1 and Syna, suggestive of impaired placental ontogeny (Fig. 4f). To investigate this possibility further, we scored the ratio of placental mass relative to embryo at E18.5. We found a significant change in the F₁ fetal-to-placental ratio when sired by nABX males (P = 0.0004, unpaired t-test; CON n = 82 nABX n = 53), which is driven specifically by reduced placenta mass (Fig. 4g). This is consistent with a placental defect derived from dysbiotic fathers.

To decipher the molecular aetiology underlying this, we noted that the top ten DEGs in nABX-derived placenta included several clinical markers for human placental insufficiency disorders, such as pre-eclampsia (Fig. 4d and Extended Data Fig. 12c). For example, Plgf was significantly downregulated in mature placentae sired by dysbiotic fathers, whereas upregulated pre-eclampsia markers include the Flt:Plgf ratio, Clu and Afp (Fig. 4h). To further investigate the potential for dysbiotic fathers to induce placental insufficiency, we examined placenta structure. Here, placenta derived from nABX fathers were associated with a significantly reduced labyrinth zone (P = 0.0098) (Fig. 4i), which is a frequent underlying cause of placental disorders⁴⁴. Indeed, we further found significantly impaired vascularization (P = 0.0076) and increased placental infarction (P = 0.0296) (Extended Data Fig. 12d–g). Finally, we assayed the level of placental growth factor (PLGF) hormone, which when low, is a prime diagnostic marker for pre-eclampsia in humans⁴⁵. We found that PLGF was significantly lower in placenta when offspring were sired by dysbiotic males induced by either nABX (P = 0.034) or avaABX (P = 0.0006), while the sFLT/PLGF ratio was significantly elevated (Fig. 4j).

Discussion

Taken together, cumulative evidence suggests that environmentally induced perturbations to the gut microbiota elicit a significant reproductive response in prospective fathers. This supports a regulatory gut–germline axis, which when perturbed can propagate to affect offspring disease risk, at least in part by impacting forthcoming placental function. The gut microbiota may therefore act as a principal interface wherein distinct environmental inputs, such as antibiotic regimes or diet, can converge and signal to male germ cells directly or indirectly, ultimately exerting influence on progeny. Such paternal F₁ effects manifest probabilistically however, with multifactorial molecular mechanisms acting in/on sperm probably interacting to underpin transmission. Although we identify altered lipid metabolites, hormones such as leptin and small RNA payloads here, future work is necessary to deconvolve the phenotypically relevant modalities of inheritance and their applicability beyond murine models. However, our observation that paternal conditioning can impact placental ontogeny provides a mechanistic grounding for mammalian intergenerational effects⁴⁶. Moreover, because restoration of paternal gut microbiota before conception rescues emergent F₁ phenotypes, the effect is reversible and thus remedially tractable. Given the prevalence of lifestyle and antibiotic practices that (re)shape human microbial communities^23,24,25,47, this could prove an area of interest towards mitigating against adverse pregnancy outcomes. More generally, our data highlight the importance of understanding how environmental factors can modify complex biological systems across scales; from direct molecular responses to intergenerational disease susceptibility.

Methods

Animal husbandry

All experiments involving mice were carried out in accordance with the approved protocol and guidelines by the laboratory animal management and ethics committee of the European Molecular Biology Laboratory (EMBL) under license no. 20190708_JH and the Italian Ministry of Health under authorization code no. 308/2021-PR. The inbred C57BL/6J strain was used as the primary mouse model, whereas CD1 IGS or C57BL/6J dams were used as surrogate mothers for IVF. C57BL/6J mice deficient in leptin (ob/ob mice) were used to study the intergenerational effects of paternal testicular leptin deficiency. Mice were housed under a 12 h light/dark cycle (from 07:00 to 19:00), with ad libitum access to regular diet and water. Standard chow contained 18.5% protein, 5.3% fat, 4% fibre and other nutritional additives in pellet form and is suitable for long-term maintenance, breeding, lactation and gestation periods (NFM18, Mucedola).

Experimental strategy for induced dysbiosis in animals

At 5 weeks of age, specific pathogen free (SPF) male C57BL/6J inbred mice were transferred to a conventional colony facility. The male mice were divided into separate cages (one mouse per cage) and randomly assigned to two groups: (1) control group, which received regular diet and sterile filtered (0.2 μm) water for six successive weeks (abbreviated as CON mice); (2) treated group, which received regular diet and sterile filtered (0.2 μm) water supplemented with a cocktail of FDA-approved nABX (neomycin 2.5 mg ml⁻¹, bacitracin 2.5 mg ml⁻¹ and pimaricin 1.25 μg ml⁻¹; abbreviated as nABX) for 6 successive weeks. The water bottle and nABX cocktail were replaced freshly every 5 days. At the end of 6 weeks of treatment, male mice from both groups were mated with treatment-naive 6-week-old C57BL/6J inbred female mice. Subsequent to mating, male mice were maintained for an extra 8 weeks without nABX (recovery period) and they had ad libitum access to regular diet and sterile filtered (0.2 μm) water. During the 8 week recovery period, each male mated at two time points: at 4 weeks after nABX withdrawal and at 8 weeks postwithdrawal of nABX (denoted as 6wk + 4rec and 6wk + 8rec, respectively). All the female mice received regular diet and water throughout gestation and lactation period. Of note, the dosage of each antibiotic used is lower than the standard human-equivalent dosage and half of the standard dose administered to mice in previous studies and therefore constitutes a relatively low dosage regime designed to perturb rather than ablate gut microbiota communities^30,48,49.

An equivalent experimental setup and 6 week time course was used for alternative strategies to induce acute gut dysbiosis. First, ad libitum administration of a cocktail of absorbable antibiotics (ampicillin 0.5 mg ml⁻¹, vancomycin 0.25 mg ml⁻¹ and amphotericin B 12.5 μg ml⁻¹; abbreviated as avaABX) in filtered (0.2 μm) water. Second, through bowel cleansing using 5% PEG (PEG 4000, EMD Millipore). Herein, the sire males were given avaABX or 5% PEG in drinking water for 6 weeks to provoke acute gut microbiota dysbiosis, whereas the randomly assigned control group received only sterile filtered (0.2 μm) water. At the end of 6 weeks, male mice from each group mated with treatment-naive 6-week-old C57BL/6J inbred female mice. For each antibiotic, the administered dose was adjusted to half the standard antibiotic dose that was used in previous studies^50,51,52. Likewise, the PEG concentration was lower than the standard dose used for human bowel cleansing or for medicating constipation and in previous mouse studies³¹.

Male fertility parameters

To evaluate the effect of nABX-induced dysbiosis on breeding males, the testis-to-body weight ratio and standard male fertility parameters were characterized on the basis of a quantitative descriptive analysis which includes, sire body weight, testis weight, sperm count and fecundity potential:

• 1. Testes-to-body weight ratio: to assess the effect of nABX treatment, the body weights of sires treated for 6 weeks were measured and compared with those of age-matched CON sires. Furthermore, to assess the impact on reproductive organs, testes were carefully collected from CON and nABX-treated mice, weighed and the testes-to-body weight ratio difference between the two groups was calculated.

• 2. Sperm count: the method described by ref. ⁵³ was used for sperm counts. In brief, sperm suspensions were prepared by mincing cauda in 1 ml of 1× PBS (pH 7.2). The suspension was pipetted and filtered through 80 μm nylon mesh to remove tissue fragments. An aliquot (0.05 ml) from the sperm suspension (1 ml) was diluted with 1:40 PBS (pH 7.2) and mixed thoroughly. The sperm suspension was counted using Bürker chamber after discharging a few drops of the suspension. To determine sperm count per millilitre, sperm counts in each of the four squares (1 mm²) were averaged and multiplied by the dilution factor ×10⁴ to calculate the total number of sperm per mouse.

• 3. Fecundity potential: 11-week-old sire males (both CON and nABX-treated) were mated with treatment-naive 6-week-old C57BL/6J females. The females were checked for the presence of a copulatory plug over the first 4 days and moved to a holding cage once they were plug positive, then checked again at 3 weeks for pregnancy and litter size.

Offspring generation

To generate F₁ offspring both natural mating and artificial (IVF) were used. A coding system was used to ensure that researchers and animal husbandry staff were double-blind to the paternal treatment regime.

Natural mating

To generate F₁ offspring, each CON or nABX male mouse was placed with a single 6-week-old treatment-naive virgin female. During the mating period (1–4 days), mice were housed under standard optimized environmental conditions with ad libitum access to regular diet and water. Each morning females were checked for vaginal status, all females with positive signs of mating (vaginal plugs) were removed from the male cage and housed in a new cage, whereas all females that did not show a vaginal plug were returned to the male cages in the afternoon and examined on each successive morning. To minimize microbial transfers between the pairs housed together, the mating hours of the pairs was restricted from late afternoon to early morning (15:00–09:00) throughout the mating period.

In vitro fertilization

As described previously⁵⁴, this method of F₁ offspring generation follows five critical steps: (1) superovulation of oocyte-donor females; (2) preparation of fresh sperm from sire males; (3) collecting isogenic oocyte from donor females; (4) in vitro sperm–oocyte fertilization; (5) embryo transfer to pseudopregnant CD1 or C57BL/6J surrogate mothers. This process was performed as follows:

• 1. Superovulation: 4-week-old C57BL/6J oocyte-donor female mice were injected intraperitoneally with 0.2 ml of PMSG plus inhibin antiserum (IAS). Precisely 48 h later, hCG (0.15 ml) was injected intraperitoneally into oocyte donors. Herein, hCG injection was done according to the planned sperm–oocyte fertilization time (approximately 13–14 h after hCG injection).

• 2. Preparation of fresh sperm: cauda epididymis were dissected from 11-week-old mice, with any excess adipose tissue and blood cleaned away. Subsequently, cauda epididymis was transferred into a centre-well organ culture dish (Falcon) containing 120 μl of cryoprotective agent (CPA). Five to six incisions were made across each cauda epididymis with optimized scissors and the surface of the cauda was gently pressed to release the sperm. Fresh sperm were pre-incubated in CPA for 3 min at 37 °C with gentle stirring of the Petri dish every minute. A 10 μl sperm suspension was transferred into the centre drop of 90 μl of TYH-MBCD medium and the fresh sperm were capacitated by incubation for at least 10–30 min at 37 °C in the 5% CO₂ incubator. Dishes remained undisturbed until the sperm were moving rapidly and equilibrated in the medium.

• 3. Oocyte collection: a day before oocyte collection, 90 μl of HTF medium containing 1 mM GSH was covered with paraffin oil and placed overnight at 37 °C in the 5% CO₂ incubator to allow equilibration. The next day, the oviduct isolated from each superovulated female (ten donors in total) was divided into two, half used for fertilizing with sperm from CON mice and half with sperm from nABX-treated mice, which assured isogenic oocytes source and eliminated oocyte-donor quality as a confounding effect. Once each half of the oviduct ampullae were transferred into fertilization dish containing HTF drops, swollen ampullae were clipped with a 18G needle to release the cumulus masses into the oil and then transferred into the fertilization drop.

• 4. Sperm–oocyte fertilization: from capacitated sperm in TYH-MBCD drops, a 10 μl of sperm suspension was taken from the peripheral part of the pre-incubation drops which contain the most motile sperm. Aliquoted sperm were then transferred either from CON mice or nABX-treated mice into the centre of the drops in the fertilization dish containing oocytes. After insemination of the oocytes with sperm, sperm–oocyte fertilization dishes were cultured at 37 °C under 5% CO₂ for 3–4 h. After 3–4 h of incubation, the presumptive zygotes were washed three to four times with 150 μl drops of M2 medium and cultured overnight in four-well NUNC plates with 500 μl of KSOM medium supplemented with amino acids (KSOMaa) at 37 °C under 5% CO₂ and 5% O₂ incubator. After 18–24 h, the fertilized two-cell stage embryos were transferred to the oviducts of pseudopregnant CD1 or C57BL/6J female surrogate mice.

• 5. Embryo transfer: embryo transfer was conducted as described previously⁵⁵. Vasectomized males were mated with 6–8-week-old CD1 females or C57BL/6J (27–30 g) to produce 0.5 days postcoitum (dpc) pseudopregnant recipients. Plug-positive females were selected and anaesthetized by inhalation of isofluorane. To transfer two-cell stage embryos, the following subsequent steps were applied: the dorsal midline opened, ovary body wall incised and held outside the body, after breaking the bursa the tip of the capillary inserted into the infundibulum to transfer ±20 embryos per single pseudopregnant CD1 or C57BL/6J mouse and finally a gentle push applied through mouth pipette to dislodge embryos from the glass capillary tube into the infundibulum of one oviduct.

F₁ phenotype after co-housing to examine indirect effects

To separate the potential indirect impact on offspring of maintaining females with a dysbiotic male (such as microbial transfer to mothers or coprophagy of nABX faeces), we used a cohousing strategy coupled with control (independent) male mating. A male with either 6 weeks of previous nABX exposure or a control male was placed with one 6-week-old treatment-naive virgin female during the day (09:00–16:00), recapitulating full parental co-exposure. Males were removed from the cages at 16:00 to prevent mating during the night, with females remaining in the male cages to maximize the potential to detect microbial transfer or coprophagic effects. Each day at 16:00, females were checked for vaginal status, all females with positive signs of mating (vaginal plugs) were excluded from the study, whereas all females that were vaginal plug negative were held in the male cages. In the 4 day cohousing period, mice were maintained in standard optimized environmental conditions with ad libitum access to food and water. After 4 days of cohousing, all plug-negative females (exposed to nABX or CON males) were mated with conventionally raised 11-week-old C57BL/6J male mice (independent naïve males) to determine whether any indirect effects of previous nABX male exposure could be detected in F₁ phenotypes.

F₂ offspring generation

F₂ offspring were produced by intercrossing two siblings from the F₁ generation derived from either dysbiotic (nABX) or CON F₀ sires.

Phenotyping and tissue collection

Offspring body weight and survival metrics

A comparative growth phenotype analysis was applied to evaluate F₁ prenatal to postnatal development. Matings were setup using a blinded code system and F₁ phenotypes were recorded by individuals and/or husbandry staff without knowledge of the paternal condition. For postnatal analysis, the litter size at birth was recorded and still-born pups were monitored. The body weight of pups was measured every 3 days from postnatal day P3 until weaning age (P3–P21) and every week up until P56. F₁ offspring morbidity and mortality were also recorded daily in this age range (P0–P60). To account for outliers due to extreme litter size effects, only litters within 2 s.d. of the average litter size (6 ± 2 pups per litter) were included in the study (that is, litters of 4–8 pups), which corresponded to about 81% of total litters. Likewise, the postnatal body-weight trajectories of F₁ progeny generated by IVF were characterized and analysed using aforementioned methods. Note, F₀ pregnant females as well as F₁ offspring were fed a regular diet and sterilized water throughout gestation, lactation and the postweaning period. To phenotype the prenatal F₁ conceptus, derived from CON or nABX-treated males, pregnant dams were killed either at E13.5 or at E18.5. Placenta and fetus were carefully isolated, cleaned and weighed. Subsequently, embryonic tissues and placentae were isolated for transcriptome profiling. Placental tissue was further collected for immunofluorescent or ELISA analysis.

Tissue and sperm isolation for omics

Testes and epididymal sperm were isolated as described previously⁴², to ensure maximum purity and quality. Briefly, testes were carefully dissected from 11-week-old mice, with all fat pads removed, before being individually weighed, snap-frozen in liquid nitrogen and stored at −80 °C until processing for transcriptomics or ELISA assay. The contralateral testis was processed for histological staining. For sperm collections, cauda epididymis was dissected from mice and placed in M2 media at warmed 37 °C. Two small incisions were made at the proximal end of cauda using a 26G needle and holes were poked in the rest of the tissue to let the sperm swim out and epididymal fluid release. Sperm-containing media were incubated for 30 min at 37 °C, then transferred to fresh tubes and incubated for another 15 min at 37 °C. After incubation, sperm were collected by centrifugation at 2,000g for 5 min, followed by a 1× PBS wash and a second wash with somatic lysis buffer (contains 0.1% SDS and 0.5% Triton-100X diluted in distilled H₂O) for 20 min in ice, which is important to eliminate any potential somatic cell contamination. Somatic lysis buffer-treated sperm samples were then collected by centrifugation at 3,000g for 5 min and finally washed twice with 1× PBS and stored at −80 °C. Purification of sperm cells was confirmed by microscopic examination.

Testis single-cell isolation for 10× single-cell RNA sequencing

Single testicular cells were isolated using a two-step enzymatic digestion protocol described previously⁵⁶, with minor modifications. In brief, testes from nABX and CON males were excised, washed twice in PBS, the tunica albugineas removed, each testis transferred to 1 ml of lysis buffer (1 mg ml⁻¹ of collagenase A diluted in PBS) and triturated five times using a P1000 tip to disrupt the tissue and incubated for 5 min at 37 °C with gentle shaking at 300 rpm. The tubules were gravity sedimented, centrifuged for 5 min at 300g at room temperature. The supernatant was removed and remaining tissue washed with PBS and digested in TrypLE (1×) and DNase I. The suspension was triturated vigorously ten times, incubated at 37 °C for 5 min, followed by repeat trituration and incubation until the end of 15 min. The digested single-cells suspension was sequentially size-filtered and washed through 70 and 40 μm strainers and pelleted by centrifugation at 300g for 5 min. The pelleted cells were resuspended in 1 ml of cold PBS + 0.04% BSA, pooled into a Falcon tube (four testis per group) and counted using Countess II (average viability 85–90%). Target cell concentration was adjusted in the range of about 1,100 cells μl⁻¹ using cold PBS with 0.04% BSA and placed on ice until loading on the 10x Chromium chip.

Seminal fluid collection

Seminal vesicles were collected in a laminar flow-hood as described previously⁵⁷ with slight modifications. To prevent cross contamination, one experimental animal (CON or nABX) was placed in the hood per sample collection. To ensure maximum sterilization, the experimenter wore a sterile gown, gloves, mask and scrubbed their hands with 70% ethanol and followed the following procedures. In brief, the mice were placed in a laminar flow-hood that had been previously cleaned with 70% ethanol and 2% Virkon followed by 2 h of ultraviolet radiation. To collect seminal vesicles aseptically, mice were euthanized by cervical dislocation and their abdominal areas were cleaned with 70% ethanol and incised with sterilized tweezers and scissors. Seminal vesicles were then excised using new sterilized tweezers and scissors and placed in a sterilized 2 ml Eppendorf tube, snap-frozen in dry ice and stored at −80 °C until DNA extraction for seminal fluid microbiome profiling. Note, dissection tools used in this sample collection were used only once per mouse.

Blood collection

Blood samples were collected from saphenous vein microvette blood collection tubes (Microvette 100K3E SARSTEDT, catalogue no. 20.1278.100) and centrifuged at 5,000 rpm for 10 min at 4 °C; plasma was obtained and stored at −80 °C until assayed.

Placenta tissue collection

The placental tissues were collected according to the protocol described in ref. ⁵⁸. Briefly, samples were collected at E13.5 and E18.5, counting from E0.5 as the day the copulation plug was found. The conceptuses were removed from a pregnant mouse and the uterine horn containing a single conceptus was segmented. Using microscissors, the uterine muscle was gently peeled away from the antimesometrial side (the side that is least vascularized), the yolk sac was cut to expose the amniotic sac and the placenta was carefully dissected and separated from the umbilical cord and yolk sac to collect the entire placental disc without cut or tear remnants. Samples were snap-frozen in liquid nitrogen and stored at −80 °C until they were processed for RNA-seq or fixed in 4% paraformaldehyde for histology and immunohistochemical analysis.

Pre-implanatation embryos

To generate blastocyst-stage embryos, IVF was used. Briefly, 4-week-old C57BL/6J oocyte-donor female mice were superovulated with intraperitoneal injections of PMSG (0.1 ml) and hCG (0.1 ml) after 48 h. Capacitated sperm suspensions were prepared from either control or leptin-deficient (ob/ob) mice aged 6 weeks and transferred to the centre of prepared drops in a fertilization dish containing oocytes. Note control (wt/wt) and leptin-deficient (ob/ob) male sperm donors were littermates derived from heterozygous crosses (wt/ob) and thus near-isogenic. After insemination of the oocytes with isolated sperm, sperm–oocyte fertilization dishes were cultured at 37 °C under 5% CO₂ for 3–4 h. After 3–4 h of incubation, the presumptive zygotes were washed 3–4 times with 150 μl drops of M2 medium and cultured for 4.5 days in four-well NUNC plates with 500 μl of KSOM medium supplemented with amino acids (KSOMaa) at 37 °C in a 5% CO₂ and 5% O₂ incubator. After 4.5 days, each blastocyst-stage embryo was collected in 0.2 ml PCR tube strips (containing 2 µl of lysis buffer + 1 µl of dNTPs + 1 µl of oligo-dT primer) for single-embryo RNA-seq analysis (Smart-Seq). Note, lysis buffer composition: 0.5% Triton-X100 (vol/vol) in H₂O + 2 U μl⁻¹ of SUPERase In RNase inhibitor (20 U μl⁻¹; Thermo, AM2694) and Oligo-dT primer: 5′-AAGCAGTGGTATCAACGCAGAGTACT30VN-3′.

Histological staining

Testis tissue collected from 6-week-treated (nABX) or CON sire male mice were fixed in Bouin’s fixative, then treated through graded ethanol and clearing agent (xylene), followed by embedding in paraffin. Subsequently, the samples were sectioned serially at a thickness of 5–7 µm and stained with standard haematoxylin and eosin stain. Testicular morphometric analyses of the seminiferous tubules (that is, spermatogonium, spermatocyte, spermatid and Sertoli cells) and Leydig cells were performed using LMD 7000 laser microdissector microscope (Leica). The quantitative evaluation of the abnormal seminiferous tubules was done as follows: for each sample at least ten sections were analysed with more than 1,000 seminiferous tubules evaluated (normal or abnormal) and the percentage of abnormal tubules was calculated accordingly. For this analysis, the images were take using NanoZoomer S60 Digital slide scanner (Hamamatsu Photonics, C13210-01).

Immunofluorescent staining

Placenta samples collected at E18.5 were fixed in 4% paraformaldehyde at 4 °C, followed by embedding in paraffin and sectioned into 7-μm-thick slices. Sectioned samples were mounted on slides and dried at 40 °C on hot plate. The dried sections were subjected to heat-induced antigen retrieval using citrate buffer at pH 6.0, incubated overnight at 4 °C with anti-mouse VE-cadherin (Thermo Fisher Scientific; eBiosciences. catalogue no.14-1441-81; 2.5 μg ml⁻¹) and detected by anti-Rat-Alexa Fluor 568 (Thermo Fisher Scientific; catalogue no. A-11077; 4 μg ml⁻¹). Nuclei were counterstained with 4′,6- diamidino-2-phenylindole (DAPI) (Thermo Fisher Scientific; catalogue no. D1306; 5 μg ml⁻¹). Slides were washed before being mounted with ProLong Glass (Thermo Fisher Scientific; catalogue no. P36983) and imaged on a Leica THUNDER Imager Live Cell with a ×20 objective (numerical aperture = 0.8). QuPath v.0.2.1 image analysis software was used to measure areas of labyrinth zone and whole placenta.

Placental vessels at E13.5 and E18.5 were stained with isolectin B4. Briefly, sagittal 7 μm sections of placenta samples were deparaffinized with xylene and rehydrated through decreasing ethanol concentrations. Slide sections were subjected to heat-induced antigen retrieval for 10 min in 10 mM citrate pH 6.0 buffer. Thereafter, these sections were permeabilized with 0.3% Triton-X100, blocked in 5% donkey serum and incubated overnight at 4 °C with a rabbit Isolectin HRP (Sigma, catalogue no. L5391) at 2 μg ml⁻¹. Chromagen detection was performed using reagents from an ABC DAB kit (Vector Laboratories no. PK-6100) as per manufacturer’s instructions

Macromolecule analysis

ELISA assay

ELISA kits were used to determine the concentration of a specific target protein in biological samples. Frozen plasma and testis samples were used for the measurement of leptin (Lep) in F₀ sire males, whereas F₁ placenta tissues were used for the measurement of placental growth factor (Plgf) and soluble fms-like tyrosine kinase-1 (sFlt1). Mouse Leptin ELISA kit MOB00, mouse PlGF-2 ELISA Kit MP200 and mouse VEGFR1/Flt-1 ELISA Kit MVR100 (R&D Systems), were used for quantification of Lep, Plgf and sFlt1 following the manufacturer’s instruction. Reagents, standard curve dilutions and samples were prepared as described on the kits and all samples, standards and control were assayed in technical duplicate, with several biological replicates. To quantify the target protein concentrations, the data generated were interpolated by nonlinear regression model curve fit.

To quantify plasma leptin, an aliquot of plasma samples sourced from F₀ sire male blood and stored in −80 °C were thawed in ice, diluted 20-fold with Calibrator Diluent and assayed according to manufacturer’s instruction.

To quantify leptin, Plgf and sFlt1 in tissue, testis or placenta samples were thawed on ice and washed in PBS. This was followed by homogenization with pestle ‘A’ (about ten strokes) and pestle ‘B’ (about 20 strokes) in 2 ml of Dounce tissue grinder containing 1 ml of RIPA buffer (Sigma-Aldrich) mixed with protease inhibitor cocktail, kept in ice for 5 min with gentle agitation and centrifuged at 15,000 rpm for 20 min at 4 °C. After centrifugation an aliquot of clear supernatant was used for ELISA assay.

To screen for antibiotic residues and test whether nABX or its residues are detectable in distal tissues or the circulatory system we took independent and complementary approaches.

Mass spectrometry

Samples for the quantification of antibiotics in testis tissues were prepared as described in the section on Metabolomics below. Extracted and dried samples were resuspended in 20 μl of methanol:water (1:1) and 10 μl of samples were directly injected into an Agilent 6550 iFunnel qToF mass spectrometer. Calibration curves (twofold dilutions) were prepared using chemical standard resulting in the following ranges of injected amounts: bacteriocin (5.6–90.0 pmol), primaricin (7.0–112.5 pmol) and neomycin (5.3–85.0 pmol). Quantification of antibiotics was performed using the MassHunter Quantitative Analysis Software (Agilent Technologies, v.10.0) and limit of detection was determined by linear regression using R.

Functional residue screening

PremiTest kit (R-Biopharm) was used for the determination of nABX or avaABX residues in testis samples stored at −80 °C, which was performed in accordance with the manufacturer’s instructions for use in meat. The principle of the test is based on thermophilic bacterium spores (Bacillus stearothermophilus) which are highly sensitive to several antibiotics and sulfa compounds. Briefly, testis collected from 6-week-treated (nABX or avaABX) and untreated (CON) sire males were thawed in ice and washed in 1× PBS, followed by homogenization using 2 ml of Dounce tissue grinder containing 200 μl of 1× PBS with pestle ‘A’ (about ten strokes) and pestle ‘B’ (about strokes). Herein, sterilized water was used as negative control, whereas nABX cocktail diluted in water was used as a positive control at a concentration of the minimum residue detection limit (0.5 μg ml⁻¹). From the homogenized sample solution, 100 μl was pipetted onto the agar in the ampoule and allowed to stand at room temperature for 20 min for prediffusion. After the prediffusion, the sample solution was flushed out of the ampoule by washing twice with water and covered with foil to avoid evaporation. Subsequently, the test ampules were incubated in Eppendorf block heater at 64 ± 0.5 °C for approximately 3–3.5 h until the negative control turned from purple to yellow. At this stage, ampoules were removed from the block heater and results interpreted on the basis of the kit bicolour indicator.

Microbiota analysis

Faecal sample collection

Fresh faecal samples were collected from each individually housed F₀ sire males at several different time points: (1) before nABX administration at day 0; (2) at the end of the 6 week treatment before setting up mating; and (3) during the recovery period following nABX withdrawal. Faecal sample from dams and F₁ offspring were collected at weaning stage (P21). To collect fresh faecal samples, each parent or offspring for which faeces were to be collected was put into clean, bedding-free autoclaved cages with food and water for 2–3 h, faecal pellets collected with sterilized tweezer and stored immediately in −80 °C freezers until DNA extraction for microbiome profiling.

Microbial DNA extraction

All samples were stored at −80 °C until processing by 16S rRNA sequencing. Microbial DNA extraction from faecal specimens (about 200 mg from F₀ parent and about 100 mg from F₁ offspring) was performed using QIAamp PowerFecal Pro DNA Kit (QIAGEN) according to the manufacturer’s instructions. Briefly, faecal samples processed with the DNA extraction kit were added to a PowerBead beating tube and rapidly homogenized on a Vortex Adapter (QIAGEN catalogue no. 13000-V1-24) using Vortex-Genie 2 mixer (Scientific Industries). Once cells are lysed and potential inhibitors removed, total genomic DNA is captured on a silica membrane, washed and eluted for downstream gut microbiota profiling.

Determination of 16S rRNA gene copy number

The relative abundance of 16S rRNA gene copy number in mice faeces was determined from the standard curve generated on the basis of serial dilutions of the control sample by quantitative real-time PCR using F341 5′-CCTACGGGAGGCAGCAG-3′ and R534 5′-ATTACCGCGGCTGCTGG-3′ primers, as described previously⁵⁹.

Targeted amplification of the 16S rRNA V4 region (primer sequences F515 5′-GTGCCAGCMGCCGCGGTAA-3′ and R806 5′-GGACTACHVGGGTWTCTAAT-3′ (ref. ⁶⁰), was performed using the KAPA HiFi HotStart PCR mix (Roche) in a two-step barcoded PCR protocol (NEXTflex 16S V4 Amplicon-Seq Kit; Bioo Scientific) with minor modifications from the manufacturer’s instructions. PCR products were pooled, purified using size-selective SPRIselect magnetic beads (0.8 left-sized) and then sequenced at 2 × 250 base pairs (bp) on an Illumina MiSeq (Illumina) at the Genomics Core Facility, European Molecular Biology Laboratory, Heidelberg.

Raw 16S rRNA reads were trimmed, denoised and filtered to remove chimaeric PCR artefacts using DADA2 (ref. ⁶¹). The resulting amplicon sequence variants were then clustered into operational taxonomic units (OTUs) at 98% sequence similarity using an open-reference approach: reads were first mapped to a preclustered reference set of full-length 16S rRNA sequences at 98% similarity using MAPseq⁶². Reads that did not confidently map were aligned to bacterial and archaeal secondary structure-aware rRNA models using Infernal⁶³ and clustered into OTUs with 98% average linkage using hpc-clust⁶⁴, as described previously⁶⁵. The resulting OTU count tables were noise filtered by asserting that samples retained at least 1,000 reads and taxa were prevalent in at least two samples; these filters removed 58% of spurious OTUs but only 0.09% of total reads from the dataset.

Local sample diversities were calculated as OTU richness, exponential Shannon entropy and inverse Simpson index (corresponding to Hill diversities of order 0, 1 and 2 (ref. ⁶⁶) as average values of 100 rarefaction iterations to 5,000 reads per sample. Between-sample community diversity was calculated as Bray–Curtis dissimilarity. Trends in community composition were quantified using ordination methods (principal coordinate analysis, distance-based redundancy analysis) and tested using permutational multivariate analysis of variance (PERMANOVA⁶⁷, as implemented in the R package vegan⁶⁸.

Cohousing to measure microbiota cross-transfer during mating

To directly test whether microbiota cross-transfer occurred during the mating window period, we sampled faeces at two time points In this 4 day experiment, males from each group (CON and 6-week nABX) were cohoused at a 1:1 ratio in a single cage with a treatment-naive 6-week-old female C57BL/6J mouse continuously, regardless of plug positivity (single breeding pairs). During the mating period, both groups were maintained on regular diet and sterile water ad libitum. We sampled faeces (as above) at two time points: at day 0 precohousing to determine baseline microbiota and at day 4 postmating (denoted as pre- and post-mating, respectively) to examine the effect of female exposure/cohabitation with nABX males. Furthermore, to determine whether microbiota transfer occurs at mating, we sampled a copulatory plug from the female’s vaginal area directly following mating. Using swabs, control samples were collected from the working area and sterile sharp-toothed tweezers were used to collect plug samples.

RNA sequencing (transcriptomics)

mRNA-seq

For messenger RNA sequencing, total RNA was extracted from four different tissue types: sire testis, F₁ placenta, F₁ brain and F₁ (BAT), using RNeasy Protect Mini Kit (Qiagen) following the manufacturer’s instructions. The quantity and quality of total RNA was evaluated on a NanoDrop Spectrophotometer and Agilent TapeStation system (Agilent) to ensure RNA integrity number > 8.5. For library preparation, mRNA was first enriched using the NEBNext Poly(A) mRNA magnetic isolation module with 1 μg of input RNA. Purified mRNA was subsequently prepped into stranded libraries using the NEBnext Ultra II directional RNA library prep kit, following all manufacturer’s guidelines. Amplified libraries were multiplexed and sequenced on an Illumina NextSeq 500 (PE40).

Small RNA-seq

For sperm small RNA analysis, total RNA was isolated using miRVana microRNA isolation kit (Thermo Fisher Scientific). The procedure was performed according to the manufacturer’s instructions with a slight modification at the cell lysis step to ensure complete lysis of sperm. Briefly, purified sperm samples were thawed on ice and then the miRVana lysis buffer was added for 5 min at room temperature, with 40 mM dithiothreitol and 12 μl ml⁻¹ of proteinase K subsequently added into the sample. The sample was then incubated at 56 °C for 20 min on a thermoshaker at 450 rpm. After this step, miRNA homogenate was added to the lysate and the kit instructions recommended for total RNA isolation were followed. Small RNA libraries were prepared from 1 μg of this total RNA using the Illumina Truseq small RNA library preparation kit according to the manufacturer’s instructions. Libraries were sequenced on an Illumina Hiseq 2000.

Single-cell RNA-seq (10X)

Testicular single-cells resuspended in PBS + 0.04% BSA were adjusted at a concentration of about 1,100 cells μl⁻¹ for loading on the 10x Chromium chip. Cell capturing and library preparation was carried out as per kit instructions (Chromium Next GEM Single Cell 3′ v.3.1 (Dual Index) User Guide). In brief, about 6,000 cells were targeted for capture per sample; after complementary DNA synthesis, 12–14 cycles were used for library amplification. The resultant libraries were size selected, pooled and sequenced using Illumina NextSeq 2000 P2 flowcell (100CYC) sequencing protocol.

Hybridization in situ sequencing

The protocol was followed as described in ref. ⁶⁹. Briefly, mRNA transcripts were targeted from testis tissue cryosectioned at 10 μm thickness mounted on SuperFrost Plus adhesion slides. Sections were fixed in 3% formaldehyde, permeabilized with 0.1 M HCl for 5 min and washed with PBS. By applying SecureSeal Hybridization Chambers (Grace Bio-Labs) around tissue sections, mRNA was reverse transcribed overnight at 37 °C using random hexamers, RNase inhibitor and reverse transcriptase (BLIRT). Tissue sections were fixed for 40 min following reverse transcription, washed with PBS, phosphorylated padlock probes (PLPs) were hybridized at a final concentration of 10 nM/PLP and ligated with Tth Ligase (BLIRT) and RNaseH. After sections were washed with PBS, rolling circle amplification was performed overnight at 30 °C using phi29polymerase and exonuclease I (Thermo). After incubation, reagents were removed, samples washed twice in PBS and the SecureSeal chamber was removed with forceps, followed by hybridizing Bridge-probes (10 nM) for 1 h at room temperature in hybridization buffer (2× SSC, 20% formamide) in the dark, on a rocker. Following this, readout detection probes (0.1 μM) were hybridized in the dark at room temperature for 1 h in hybridization buffer, washed twice with PBS and stained with DAPI (0.5 μg ml⁻¹) for 5 min at room temperature. Sections were washed with PBS and mounted with about 20 μl of Fluoromount-G Mounting medium, overlayed with a cover slip and stored at 4 °C until imaged. Images were taken with a slideview vs200 slide scanner (Olympus) and analysed using TissUUmaps.

TaqMan assay for small RNA

Independent males from the small RNA-seq samples were used to extract purified sperm total RNA with the miRVana microRNA isolation kit, applying minor modifications to ensure complete lysis. Quantification of miRNA and tRNA was performed using predesigned (miR-141-3p) and custom-designed (tRNA-Gly-GCC-2) TaqMan assays, according to the manufacturer’s protocol (Applied Biosystems) using 7 ng of total RNA from CON (n = 4) and nABX (n = 4) mice. The quantitative PCR with reverse transcription was performed in 15 μl reactions using TaqMan Fast Advanced Master Mix, following the standard programme (10 min at 95 °C; 15 s at 95 °C; 1 min at 60 °C, for 40 cycles). snoRNA202 (Applied Biosystems) was used as the endogenous control for normalization of miRNA and tRNA levels. Relative quantities were normalized to endogenous control values and foldchange calculated by means of the 2^−∆∆Ct method. Amplification efficiency of TaqMan probes was confirmed by serial dilution of template runs.

Single-embryo RNA-seq (Smart-Seq)

Single-embryo RNA-seq was carried out using the Smart-Seq protocol, as described in ref. ⁷⁰. For quality control and data preprocessing, raw reads were aligned using STAR v.2.7.10a (ref. ⁷¹) at default parameter settings to mm10 (GRCm38) primary genome assembly. The data were quantified using the featureCounts module of Subread v.2.0.1 (ref. ⁷²). Differential expression and PCA were performed in R (v.4.1.2). Differential expression was performed by passing the raw counts into the DESeq2 (v.1.34.0) package⁷³. The vst function was used to normalize the counts at default settings and the top 500 variable genes were used for PCA. For gene ontology analysis, the Metascape web application (https://metascape.org/gp/index.html) was used to infer enrichment of gene ontology terms⁷⁴. The list of genes to be tested was uploaded to Metascape as a text file. The ‘Express Analysis’ option was selected after setting both the ‘Input as species’ and ’Analysis as species’ parameters to M. musculus.

DNA sequencing

Whole-genome bisulfite-seq

Extraction of total sperm DNA was performed using the DNeasy Blood & Tissue Kit (QIAGEN Hilden) optimized for purification of sperm genomic DNA. The DNA extraction was performed as follows. Frozen pure sperm samples stored at −80 °C were thawed in ice; 100 μl of DNA extraction buffer (20 mM Tris Cl pH 8, 20 mM EDTA, 200 mM NaCl, 4% SDS) containing 80 mM dithiothreitol and 12 μl ml⁻¹ of proteinase K was added to the sample; and then incubated at 56 °C on Eppendorf thermomixer until complete sperm lysis was assured (about 1 h), shaking occasionally during incubation to disperse the sample. After incubation, the user-developed protocol DY03 (QIAGEN) for purification of total DNA from animal sperm was applied. BS-Seq libraries were constructed according to the manufacturer’s instructions using TruSeq DNA Methylation Library Prep Kit (Illumina). Amplified libraries were multiplexed and sequenced on an Illumina NextSeq 500 (PE75).

De novo gDNA-seq

The gDNA was extracted from control or nABX-derived F₁ offspring (normal or SGR) liver using DNeasy Blood & Tissue Kit (Qiagen) following the manufacturer’s instructions. The quantity and quality of DNA was evaluated on a NanoDrop Spectrophotometer and Agilent TapeStation system (Agilent) to ensure DNA integrity number > 7. For library preparation NEBNext Ultra II DNA Library Prep kit was used, following all manufacturer’s guidelines. Amplified libraries were multiplexed and sequenced on an Illumina Hiseq 4000 (PE150) to appropriate depth.

Metabolomics

Untargeted metabolomics measurements

All chemicals for liquid chromatography–mass spectrometry (LC–MS) analysis including water and acetonitrile (LC–MS grade) were purchased from Fisher Scientific. Standards for online mass calibration were purchased from Agilent Technologies.

Sample preparation

Testis samples were collected from both CON and nABX-treated sire males across the three breeding time points: immediately after 6 weeks of treatment (11-week-old mice), 4 weeks after treatment withdrawal (15-week-old mice) and 8 weeks after treatment withdrawal (19-week-old mice). Immediately after collection, samples were snap-frozen in liquid nitrogen and stored at −80 °C. On retrieval, 300 μl of ice-cold solvent mixture (acetonitrile:methanol:water, 2:2:1) and two Tungsten Carbide Beads, 3 mm from Qiagen were added to each testicle sample. Samples were homogenized using Qiagen tissueLyser II at 30% strength for 5 min. The lysed samples were centrifuged at 12700 RCF at 4 °C for 10 min. An equal volume (125 μl) of extraction supernatant was transferred to two Nunc 96-well, V-shape plates and evaporated in the Genevac centrifugal concentrator (Genevac) at 25 °C for 2 h. Concentrated samples were stored at −80 °C and resuspended in 20 μl of same solvent mixture for LC–MS analysis.

LC–MS measurements

Chromatographic separation was performed using an Agilent InfinityLab Poroshell 120 EC-C18, 3.0 × 150 mm², 2.7 μm column and an Agilent 1290 Infinity II LC system coupled to a 6550 iFunnel qToF mass spectrometer. Column temperature was maintained at 45 °C with a flow rate of 0.4 ml min⁻¹. The following mobile phases were used: mobile phase A—water with 0.1% formic acid; and mobile phase B—acetonitrile with 0.1% formic acid. The 5 μl of sample were injected at 5% mobile phase B, maintained for 0.10 min, followed by a linear gradient to 30% B in 0.5 min, followed by a linear gradient to 95% B in 10 min and maintained at 95% B for 1 min. The column was allowed to re-equilibrate with starting conditions for 3 min before each sample injection. The mass spectrometer was operated in both positive and negative scanning mode (50–1,700 m/z) with the following source parameters: VCap, 3,500 V; nozzle voltage, 2,000 V; gas temperature, 275 °C; drying gas, 13 l min⁻¹; nebulizer, 45 psi; sheath gas temperature, 275 °C; sheath gas flow, 12 l min⁻¹; fragmentor, 130 V; and skimmer, 0 V. Online mass calibration was performed using a second ionization source and a constant flow (10 µl min⁻¹) of reference mass solvent (119.0363 and 1033.9881 m/z for negative operation mode and 121.0509 and 922.0098 for positive operation mode, respectively). Each sample was measured separately in both positive and negative ionization modes.

Metabolic feature extraction

The MassHunter Qualitative Analysis Software (Agilent Technologies, v.10.0) was used to extract metabolic feature from the acquired LC–MS data. The following settings were applied: peak filter of absolute height: 5,000 counts, limit assigned charge states to 1, only ±H⁺ charged molecules were included with compound quality scores greater than 80%. Peak alignment and identification were carried out using Mass Profiler Professional (Agilent, v.15.1) with default parameters: mass tolerance of 2 mDa or 20 ppm and retention time tolerance of 0.2 min or 2%. Extracted and aligned features were exported as .csv file for further data analysis

Statistical analyses

Statistical analyses was performed using Graphpad Prism v.8.4.3 graphical software and in R v.3.6.2. Significant differences in F₁ offspring body weight and phenotypes were determined with ‘nested’ analyses, whereby the replicate number is dictated by the number of uniquely exposed fathers (N), not by the number of offspring (n), which is greater. Thus, although we typically generate n > 150 offspring per experiment to robustly capture effect size and partially penetrant phenotypes, statistically speaking each individual derived from the same father is hierarchically nested into a single N value. For example, a nested t-test compares the means of two unmatched groups (all F₁ offspring from control or dysbiotic fathers), for which there is a nested factor in those treatment groups (shared father amongst each litter). This is necessary as using individual offspring as independent variables in intergenerational studies will lead to inflated alpha error rates and spurious significance. Testes-to-body weight ratio, fetoplacental ratio and labyrinth zone were analysed by two-tailed unpaired t-test. Odds ratios (ORs) and 95% confidence intervals (CIs) were computed using with Baptista–Pike method and the statistical significance of the ORs determined using chi-squared test. Kaplan–Meier method was applied to generate survival analysis curves compared by the log-rank (Mantel–Cox) test. ELISA calibration curves were interpolated with a Hyperbola (X is concentration) nonlinear regression model fit (R² > 0.99 was acceptable curve fit).

Bioinformatics analyses

Single-cell RNA-seq alignment and mapping

Raw reads were aligned and mapped using the count module in 10x Genomics Cellranger 6.1.2 (ref. ⁷⁵) to the mm10 transcriptome assembly (2020-A) with default parameters.

Single-cell RNA-seq quality control

All subsequent steps were performed in R (v.4.1.2) using the Seurat package⁷⁶. First, cells were filtered on the basis of three parameters—number of unique molecular identifiers (nCount_RNA; 1,000:5,000), number of unique genes detected (nFeature_RNA; 500:5,000) and mitochondrial rate (percent.mito) and ribosomal rate (percent.ribo; 0:20)—as described below. The nABX and CON samples were then clustered separately using the default Seurat clustering approach. In both conditions, clusters having no uniquely expressed genes (using the FindMarkers function with logfc.threshold = 0.5 and min.pct = 0.5) were discarded.

Single-cell RNA-seq integration

The nABX and CON samples were then integrated using the canonical correlation analysis approach at default parameter settings as recommended by the Seurat package.

Single-cell RNA-seq cell-type annotation

The integrated dataset was then clustered and annotated using uniquely expressed marker genes. As above, clusters showing no uniquely expressed genes were ignored. To define more fine-grained annotations, the somatic and germ cells were split into different Seurat objects and clustered separately.

Single-cell RNA-seq cell-type-specific differential expression

For each cell type identified in the dataset, the FindMarkers function was used to identify DEGs between nABX and CON cells using logfc.threshold = 0.25 and p_val_adj < 0.01.

RNA-seq quality control and data preparation

Raw reads were quality trimmed using Trim Galore (0.4.3.1, -phred33 --quality 20 --stringency 1 -e 0.1 --length 20). These were mapped to the mouse mm10 (GRCm38) genome assembly using RNA Star (2.5.2b-0, default parameters except for --outFilterMultimapNmax 1000) and reads with a MAPQ score less than 20 were discarded to ensure that only unique-mapping high-quality alignments were used for analysis of gene expression. The data were quantified using the RNA-seq quantification pipeline for directional libraries in seqmonk software to generate log₂ reads per million (RPM) or gene-length-adjusted (RPKM) gene expression values.

RNA-seq differential analysis

DEGs were determined using the DESeq2 package (v.1.24.0), inputting raw mapping counts and applying a multiple-testing adjusted P value (false discovery rate (FDR)) < 0.05 significance threshold. An extra foldchange filter of more than 2 was applied to generate final DEGs.

RNA-seq principal component analysis

Principle component analyses of transcriptomes were computed in seqmonk and R statistical software using all expressed genes as input. These were defined as having an RPKM > 0.1 in at least two replicates across all assayed samples.

RNA-seq gene ontology class enrichment

DEGs or gene lists of interest, were inputed into the STRING v.11.0 database⁷⁷ and extracting enrichment analysis related to Reactome and KEGG pathways, filtering by FDR rank.

Small RNA-seq quality control and data preparation

Adaptors were removed using fastx-clipper; fastq files were converted to fasta using a custom perl script and reads of 18–33 nucleotides in length were retained using a custom perl script. Reads were aligned to the mouse genome version mm10 using bowtie, reporting only the best alignment and requiring 0 mismatches (parameters -v 0 -k 1 -- best --sam). Alignment sam files were converted to bam files using samtools v.1.9 and bam files were converted to bed files using bedtools v.2. To quantitate miRNAs intersectBed -c was used to count the number of reads overlapping the positions of the known M. musculus miRNAs (from miRbase www.miRBase.org). The tRNAs were obtained by using intersectBed -c to count the number of reads overlapping a bed file documenting predicted mouse tRNA coordinates downloaded from the tRNA scan database (http://gtrnadb.ucsc.edu/genomes/eukaryota/Mmusc10/). The piwi-interacting RNA coordinates were taken from ref. ⁷⁸ and converted to mm10 using liftover. The piRNAs were quantitated by selecting RNAs between 26 and 32 nucleotides long with a U as the first nucleotide and intersecting these RNAs with the piRNA coordinates using intersectBed -c.

Small RNA-seq differential analysis

To investigate significant differences, data were processed using DESeq2 and the negative binomial test used to identify significant differences after Benjamani–Hochberg multiple test correction.

Whole-genome bisulfite-seq

Raw fastq sequences were quality- and adaptor-trimmed using Trim Galore (0.4.3.1) and reads aligned to mm10 using Bismark (0.20.0), discarding the first 8 bp from the 5′ end and the last 2 bp from the 3′ of a single-end reads. Cytosine methylation status was extracted from mapped reads using the Bismark methylation extractor tool. Genome-wide methylation calls were analysed using Seqmonk software (1.44.0) with five biological independent replicate datasets for each condition. To identify DMRs, the genome was first binned into sliding tiles containing 50 consecutive CpGs and their methylation status determined using the DNA methylation pipeline. DMRs were identified by running read-depth sensitive logistic regression (P(adj) < 0.05), with minimum of ten reads and applying a threshold for an absolute change in DNA methylation of 20%. The methylation level at specific genomic features (for example, imprints) was calculated using the DNA methylation pipeline in Seqmonk over target features.

De novo gDNA-seq

Alignment reads were trimmed using Trim Galore (v.0.6.3)⁷⁹, then the first ten 5′ bases of both reads removed with Cutadapt (v.2.3)⁸⁰. Reads were aligned to M. musculus refence genome (mm10) using bwa mem (BWA-0.7.17)⁸¹ before filtering with samtools (v.1.10)⁸² view with the flags ‘-h -F 256 -f 2 -q 30’ and deduplicated with Picard toolkit (v.2.9.0) MarkDuplicates⁸³. SNPs and small INDELs were called using GATK (v.4.1.6.0)⁸⁴ HaplotypeCaller. Variants were filtered to remove those with a PHRED-called site quality (QUAL) of less than 30, an allele frequency of less than 0.2 or low site coverage (sliding scale) in more than one individual. Variant functional region was annotated using ANNOVAR (2020Jne07) annotate_variation.pl⁸⁵. Structural variants were called with Delly2 call (v.0.8.7)⁸⁶.

Metabolomics quality control and data preparation

All statistical analyses and plotting were performed in R v.3.6.2 (ref. ⁸⁷). To exclude bad sample injection from downstream analyses, the sum of all extracted metabolite features (TIC) was compared between samples (mean = 1.600 × 10¹⁰, range = [1.419 × 10¹⁰; 1.808 × 10¹⁰]) and samples were excluded, if their TIC was not within 3 s.d. from the mean value (0 excluded samples). Missing data were imputed to a fixed threshold, set at 5,000 counts. Exact duplicated features or features falling within (1) a 0.002 amu (absolute threshold) or 20 ppm (relative threshold) window and (2) a 0.15 min (absolute threshold) or 2% (relative threshold) retention time window, were considered to be split peaks and therefore collapsed together. Correlation between testis weight and animal body weight was computed, to verify that there was no significant correlation between the two values (cor = 0.1870, P = 0.2477). Therefore, testis weight Z-scores were used to normalize area-under-the-curve intensity values, to consider variation in signal intensity derived from testis size. Features (1) at zero variance; (2) being singletons; and (3) present in less than 75% of the samples for each class, were removed. Finally, feature tables derived from positive and negative mode were collapsed together after checking for exact duplicated features or if feature falling in a 0.002 Da or 20 ppm window existed; if so, the feature with higher average intensity was retained and both annotations were retained. The resulting table included all 30 samples and a total of 10,582 metabolic features.

Metabolomics metabolite annotation

Annotation was retrieved from the Human Metabolome Database (HMDB, https://hmdb.ca/)^88,89, by searching for metabolites with the exact same mass or falling in a 0.002 amu or 20 ppm window from the mono-isotopic mass of a metabolite. When present, several annotations were retrieved. Only features being annotated were retained for the downstream analysis. Moreover, for each feature, metabolite class and superclass were retrieved, when present, from the HMDB.

Metabolomics principal component analysis

PCA was computed⁹⁰ both for the complete dataset and after stratifying it by sampling week for all metabolic features and annotated features only.

Metabolomics differential analysis

Statistical significance of the feature intensity differences was assessed using a two-sided t-test (stats::t.test function in R) of log-scaled data and P values were FDR-corrected for several hypotheses testing using the Benjamini–Hochberg procedure (stats::p.adjust function in R with BH parameter). Adjusted P values and foldchanges were visualized using Volcano plots. Mass, retention time, HMDB annotation, class, superclass and composite spectrum were retrieved for all features showing significantly different intensity between the treated and the control group and an absolute foldchange greater than 2.

Metabolomics metabolite class enrichment analysis

Metabolite class and superclass enrichment analysis between the treated and control group was calculated using a Fisher’s exact test. All P values were FDR-corrected for several hypotheses testing using the Benjamini–Hochberg procedure (stats::p.adjust function in R with BH parameter).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.