Paternal microbiome perturbations impact offspring fitness

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. 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. 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. 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. 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. 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. 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. 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. 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.