The gut microbiota contributes to the pathogenesis of anorexia nervosa in humans and mice

Study participants

The AN patients were all diagnosed by an experienced consultant psychiatrist and they met the DSM-5 criteria for the restricting or binge/purging subtypes of AN¹. Since 90–95% of individuals with diagnosed AN are females, we decided to include only women in the present study and since ethnicity may influence gut microbiota, we only included Danish Caucasian women with AN cases, recruited from three specialized centres in Denmark from 1 September 12014 to 31 July 2016. The centres were: Center for Eating Disorders (Odense University Hospital), Child and Adolescent Psychiatric Unit (Aarhus University Hospital) and Unit for Psychiatric Research (Aalborg University Hospital). Exclusion criteria comprised antibiotic or antifungal treatment within the previous 3 months, any acute or chronic somatic diseases or infections. All the included patients were treated in specialized centres, and they were interviewed by an experienced and specialized psychologist or psychiatrist at the start of their treatment. The validated Eating Disorder Inventory (EDI, details given in Supplementary Note 2) was used for the interview and as a questionnaire filled out by trained health professional specialists²¹. The exclusion criteria for the age-matched healthy control women were BMI below 18.5 or above 25 kg m⁻², regular medication of any kind apart from birth control pills, and antibiotics within the last 3 months. The control participants were recruited via public advertisement and via direct contact to health staff, medical students and their relatives. BMI and other clinical characteristics of the healthy controls are listed in Supplementary Table 2.

The study protocol was registered at ClinicalTrials.gov (NCT02217384) and the study was conducted in accordance with the Helsinki declaration and approved by the Regional Scientific Ethical Committee for Southern Denmark (file no 42053 S-20140040). All participants involved in this study provided written informed consent.

In-patients. All the below-described measurements were conducted during routine treatment in the specialized unit for somatic and psychological stabilization of patients with severe AN. Safe and effective weight restoration of 2.0–3.0% per week is the goal of the treatment in the inpatient unit. Care was given by a multidisciplinary team and is in accordance with international guidelines. Individually customized meals were given under the supervision of trained nurses or dietitians at scheduled times. If the meals could not be consumed within the predefined timeframes (15 min for a snack and 30 min for a main meal), supplemental nutrition drinks were added either orally or via a duodenal tube. To account for individual preferences, a choice of three different meals was offered. The macronutrient content was consistent and within recommended energy percent ranges: 40–50% carbohydrate (maximum of 10% sugar), 30–40% fat and 20–25% protein. The daily energy intake was individualized according to the weight course. If a patient failed to reach 2% of weekly weight gain, the energy content of the daily menu was increased. All meals were followed by a supervised rest varying from 30 to 60 min in a seated position. Between the rests, light physical activity such as a walk was allowed; however, no forms of exercise training were allowed. For patients with an urge for excessive exercise or otherwise a lack of compliance with behaviour rules, behaviour supervision was extended and if needed, to 24 h a day.

Out-patients. AN patients paid regular visits to outpatient units where care was given by a multidisciplinary team involving medical doctors, nurses, psychologists and health behaviour educators.

They were given an individual diet plan, which had the same macronutrient composition as mentioned above for in-patients. As with the in-patients, all out-patients were also enrolled in cognitive behavioural therapeutic treatment courses.

Height was measured on a wall-mounted stadiometer and weight was measured in the morning before breakfast on a calibrated platform scale. BMI was calculated as weight divided by the square of height (=kg/m²).

Biochemical analyses

Blood samples were taken in the morning after an overnight fast. The blood samples were collected on ice and processed to obtain serum and plasma, and subsequently stored at −80 °C. Serum concentrations of sodium, potassium, albumin and creatinine were measured by enzymatic assays on a Roche/Hitachi cobas c system. Serum concentrations of total cholesterol, high-density lipoprotein cholesterol and low-density lipoprotein cholesterol were determined using the phosphotungstic acid magnesium chloride precipitation method. Serum immunoreactive insulin levels were measured using an enzyme-linked immunosorbent assay, while serum concentration of alanine aminotransferase was analysed with an enzymatic coulometric method including pyridoxal phosphate activation.

Analysis of plasma ClpB was performed as previously described⁵⁵.

Faecal sample collection and DNA extraction

Stools were collected according to International Human Microbiome Standards (IHMS) guidelines (SOP 03 V1) in kits by AN cases and HC at home and immediately frozen at −20 °C until they were transported on dry ice and frozen 4–24 h later at −80 °C in plastic tubes at the biobanks. DNA extraction from aliquots of faecal samples was performed following IHMS SOP P7 V2⁵⁶.

Bacterial cell counting

For bacterial cell counting (Supplementary Fig. 7), 0.08–0.12 g of frozen (−80 °C) faecal samples were diluted 15 times in pH 7.2 Dulbecco’s phosphate-buffered saline (DPBS) (Sigma-Aldrich), mechanically homogenized using tissue lyser (40 min, 12.5 agitations per second; QIAGEN) and fixed with 2% paraformaldehyde (10 min, room temperature; Biotum). Then the samples were diluted 120 times in filtered staining buffer (1 mM EDTA, 0.01% Tween20, pH 7.2 DPBS, 1% BSA (Sigma-Aldrich)). To minimize clumps, the samples were filtered through a cell strainer (pore size 5 μm; pluriSelect), pre-wet in the staining buffer. Next, the bacterial cell suspension was stained with SYBR Green I (1:200,000; Thermo Fisher) in DMSO (Sigma-Aldrich) and incubated in the dark for 30 min. For accurate determination of bacterial cell counts, a known concentration of 123count eBeads (Invitrogen) was added to the samples before the analysis. Measurements were performed using a BD Fortessa LSRII flow cytometer (BD Biosciences) and data were acquired using BD FACSDiVaTM software. A threshold value of 200 was applied on the FITC (530/30 nm) channel. Fluorescence intensity at green (530/30 nm, FITC), blue (450/50 nm, Pacific Blue), yellow (575/26 nm, PE) and red (695/40 nm, PerCP-Cy5-5) fluorescence channels as well as forward- and side-scattered (FSC and SSC) light intensities were collected. Measurements were performed at a pre-set flow rate of 0.5 μl s⁻¹. Data were processed in R using flowcore package (v1.11.20)⁵⁷ in R (v4.1.2). Fixed gating strategy separated the microbial fluorescent events from the faecal sample background.

Shotgun sequencing

DNA was quantitated using Qubit fluorometric quantitation (Thermo Fisher) and qualified using DNA size profiling on a fragment analyser (Agilent). High molecular weight DNA (>10 kbp; 3 µg) was used to build the library. Shearing of DNA into fragments of approximately 150 bp was performed using an ultrasonicator (Covaris) and DNA fragment library construction was performed using the Ion Plus Fragment Library and Ion Xpress Barcode Adapters kits (Thermo Fisher). Purified and amplified DNA fragment libraries were sequenced using the Ion Proton Sequencer (Thermo Fisher), with a minimum of 20 million high-quality reads of 150 bp (on average) generated per library.

Construction of a gene count table

To construct the gene count table, METEOR software was used⁵⁸: first, reads were filtered for low quality by AlienTrimmer⁵⁹. After removal of low-quality reads and human DNA reads, 75.7% ± 2.7% high-quality metagenomics sequencing reads of faecal DNA were mapped onto the Integrated Gut Catalog 2 (IGC2)⁶⁰, comprising 10.4 million genes, using Bowtie2 (ref. ⁶¹). Reads mapped to a unique gene in the catalogue were attributed to their corresponding genes. Then, reads that mapped with the same alignment score to multiple genes in the catalogue were attributed according to the ratio of their unique mapping counts to the captured genes. The resulting count table was further processed using the R package MetaOMineR v1.31⁶². It was downsized at 14 million mapped reads to take into account differences in sequencing depth and mapping rate across samples. Then the downsized matrix was normalized for gene length and transformed into a frequency matrix by fragments per kilobase of transcript per million fragments mapped normalization. Gene count was computed as the number of genes present (abundance strictly positive) in the frequency matrix.

Profiling and annotation of MSPs and gut enterotypes

IGC2 was previously organized into 1,990 MSPs with MSPminer⁶³ using a publicly available updated MSP dataset⁶⁴. Relative abundance of MSP was computed as the mean abundance of its 100 ‘marker’ genes (that is, the genes that correlate the most altogether). If less than 10% of ‘marker’ genes were seen in a sample, the abundance of the MSP was set to 0. This approach was used in the MetaHIT⁶² and Metacardis⁶⁵ consortia. For the 4 MSPs with less than 100 core genes, all available core genes were used.

Abundances at higher taxonomical ranks were computed as the sum of the MSP that belong to a given taxa. MSP count was assessed as the number of MSPs present in a sample (that is, whose abundance is strictly positive). Enterotypes profiling was performed as previously demonstrated⁶⁶.

Estimating functional modules of gut bacteriome

Genes from the IGC2 catalogue were mapped with diamond⁶⁷ onto KEGG orthologues (KO) from the KEGG database⁶⁸ (v8.9). Each gene was assigned to the best-ranked KO among hits with e-value < 10 × 10⁻⁵ and bit score >60. Then we assessed presence and abundance of GMMs²⁸ and GBMs²⁹ in a metagenomic sample by the pipeline implemented in the R package omixerRpm (v0.3.2) as previously described^28,29.

Estimation of dynamic growth rate of bacteria from metagenomic samples

We used the computational pipeline to infer gut bacterial growth dynamics from metagenomic samples as previously described²⁷. Sequencing reads were mapped to a database that contains complete genomic references of 2,991 strains belonging to 1,509 microbial species. For each bacterial species, a reference strain with prevalence of 100% across the samples was selected. A coverage map was then assembled on the basis of aligned reads to the reference genome. Genomic segments were binned into 10-kbp regions and coverage of the resulting bins was calculated and smoothed. The location of the origin and terminus of replication was predicted by fits of the same strain across multiple samples. Lastly, PTRs were calculated for each bacterial species in every sample as the smoothed sequencing coverage of the representative strain at the predicted peak location, divided by that at the predicted trough location.

Studies of bacterial structural variations

Before the SVs classification, the pipeline with iterative coverage-based read assignment algorithm was performed to reassign the ambiguous reads to the most likely reference with high accuracy^31,32. The reference genomes provided in the proGenomes database (http://progenomes1.embl.de/) were concatenated and then divided into genomic 1-kbp bins and applied for the detection of highly variable genomic segments. The SGV-Finder pipeline³¹ was used to detect the SVs that are either (1) with deletion percentage of the genomic segment across the population of <25% (variable SVs, vSVs), (2) with deletion percentage between 25% and 75% (deletion SVs, dSVs; the absence or presence of this particular genomic segment was kept) or (3) with deletion percentage of >75% (this genomic segment was excluded from the analysis). All bacterial species with SV calling were present in at least 10% of the total samples and were used for subsequent analysis.

Mediation analysis

The R package ‘mediation’⁶⁹ (v4.5.0) was used to infer causal relationships between gut microbial features, polar and microbiota-related metabolites, and metabolic traits and eating disorder scores. To reduce the testing numbers, we only kept the candidate groups consisting of variables that were strongly associated with each other; that is, for a candidate microbial feature-metabolite-phenotypic variable causal group, the association between gut microbial feature and serum metabolite was significant (P_adj < 0.1); the association between metabolite and phenotypic variable was significant (P_adj < 0.1); and the association between microbial feature and phenotypic variable was significant (P_adj < 0.1). After performing mediation analysis, only candidate groups with significance in direction 1 were kept for Sankey diagram visualization.

Profiling and analysis of viral gut microbiota

We profiled the viral gut microbiota using MiCoP⁷⁰, as this method is optimized to call viruses directly from bulk metagenomics sequencing reads and compute relative abundance within the virome dataset. As a reference dataset, MiCoP draws upon the NCBI’s RefSeq Viral database⁷¹. We identified a total of 209 viral species with prevalence of > and =10% and relative abundance of > and =0.01% for 147 (77 AN versus 70 HC) individuals included in the dataset. Richness, alpha and beta diversity were calculated with the R package ‘fossil’⁷² and ‘vegan’⁷³. Two-tailed Wilcoxon’s rank-sum test was used to determine statistically significant differences in richness and alpha diversity indices between groups. Permutational multivariate analysis of variance (PERMANOVA) at n = 999 was performed for Canberra distance. The viral–bacterial interactions in both AN-RS and AN-BP microbiome data were computed using the Sparse Correlations for Compositional (SparCC)⁷⁴ algorithm. Before the SparCC analysis, the AN bacterial and viral microbiota datasets were subset to AN-RS and AN-BP datasets, which were then separately submitted for SparCC analysis.

Analysis of serum polar metabolites by gas chromatography–time-of-flight mass spectrometry

The metabolites listed as gut microbiota-related metabolites were based on literature mining^75,76. Serum samples were randomized and sample preparation was done as described previously^43,77. Briefly, 400 μl of methanol (MeOH) containing internal standards (heptadecanoic acid, deuterium-labelled dl-valine, deuterium-labelled succinic acid and deuterium-labelled glutamic acid, 1 µg ml⁻¹) was added to 30 µl of the serum samples, which were then vortex mixed and incubated on ice for 30 min. Samples were then centrifuged (9,400 × g, 3 min) and 350 μl of the supernatant was collected after centrifugation. The solvent was evaporated to dryness, 25 μl of MOX reagent was added and the sample was incubated for 60 min at 45 °C. N-Methyl-N-(trimethylsilyl)trifluoroacetamide (25 μl) was added and after 60 min incubation at 45 °C, 25 μl of the retention index standard mixture (n-alkanes, 10 µg ml⁻¹) was added.

The analyses were done using an Agilent 7890B gas chromatograph coupled to an Agilent 7200 quadrupole time-of-flight mass spectrometer. The following parameters were used: injection volume was 1 µl with 100:1 split on PTV at 70 °C, heating to 300 °C at 120 °C min⁻¹; column: Zebron ZB-SemiVolatiles with length of 20 m, inner diameter of 0.18 mm, film thickness of 0.18 µm, with initial helium flow of 1.2 ml min⁻¹, increasing to 2.4 ml min⁻¹ after 16 min. Oven temperature programme: 50 °C (5 min), then to 270 °C at 20 °C min⁻¹ and then to 300 at 40 °C min⁻¹ (5 min). EI source: 250 °C, 70 eV electron energy, 35 µA emission, solvent delay 3 min. Mass range 55 to 650 amu, acquisition rate 5 spectra per second, acquisition time 200 ms per spectrum. Quad at 150 °C, 1.5 ml min⁻¹ N₂ collision flow, aux-2 temperature 280 °C.

Calibration curves were constructed using alanine, citric acid, fumaric acid, glutamic acid, glycine, lactic acid, malic acid, 2-hydroxybutyric acid, 3-hydroxybutyric acid, linoleic acid, oleic acid, palmitic acid, stearic acid, cholesterol, fructose, glutamine, indole-3-propionic acid, isoleucine, leucine, proline, succinic acid, valine, asparagine, aspartic acid, arachidonic acid, glycerol-3-phosphate, lysine, methionine, ornithine, phenylalanine, serine and threonine purchased from Sigma-Aldrich at a concentration range of 0.1–80 μg ml⁻¹. An aliquot of each sample was collected, pooled and used as quality-control samples, together with National Institute of Standards and Technology (NIST) CRM1950 serum sample, an in-house pooled serum sample. The relative standard deviation of the concentrations was on average 16% for the pooled quality-control samples and 10% for the NIST samples.

Analysis of serum bile acids and serum semipolar metabolites

The sample preparation procedure was performed as described previously⁷⁸. The plate was preconditioned with 450 µl acetonitrile before the addition of 100 µl of sample and 10 µl of polyfluoroalkyl substances (PFAS) and bile acids (BAs) internal standard mixture (200 ng ml⁻¹ and 1,000 ng ml⁻¹, respectively). Thereafter, 450 µl of acetonitrile containing 1% formic acid were added to each well and the samples extracted using a 10″ vacuum manifold. The eluate was evaporated to dryness under nitrogen gas flow and reconstituted to 80 µl of MeOH/2 mM aqueous ammonium ethanoate.

Chromatographic separation was carried out using an Acquity UPLC BEH C18 column (100 mm × 2.1 mm inner diameter, 1.7 µm particle size), fitted with a C18 precolumn (Waters). Mobile phase A consisted of H₂O:MeOH (v/v 70:30) and mobile phase B of MeOH with both phases containing 2 mM ammonium acetate as an ionization agent. The flow rate was set at 0.4 ml min⁻¹ with the elution gradient as follows: 0–1.5 min, mobile phase B was increased from 5% to 30%; 1.5–4.5 min, mobile phase B was increased to 70%; 4.5–7.5 min, mobile phase B was increased to 100% and held for 5.5 min. A post-time of 5 min was used to regain the initial conditions for the next analysis. The total run time per sample was 18 min. The dual electrospray ionization source settings were as follows: capillary voltage was 4.5 kV, nozzle voltage 1,500 V, N₂ pressure in the nebulizer was 21 psi and the N₂ flow rate and temperature as sheath gas were 11 l min⁻¹ and 379 °C, respectively. To obtain accurate mass spectra in the mass spectrometry (MS) scan, the m/z range was set to 100–1,700 in negative ion mode. MassHunter B.06.01 software (Agilent) was used for all data acquisition.

Identification of compounds was done with an in-house spectral library using MS (and retention time), tandem mass spectrometry information. Quantitation was based on a matrix-matched calibration curve spiked with native compounds. The calibration curve consisted of concentrations ranging from 0–1,600 ng ml⁻¹ for BAs. The relative standard deviation for the BAs was on average 17.8% for the quality-control samples and 19.4% for the NIST samples.

Animal experiments

Animal protocols were approved by the Science Ethics Committees of the Capital Region of Copenhagen, Denmark. Female germ-free Swiss Webster mice were bred and maintained in flexible film gnotobiotic isolators until the start of experiments at the Department of Experimental Medicine, University of Copenhagen. The mice were fed autoclaved chow diet (7% simple sugars, 3% fat, 50% polysaccharide, 15% protein (w/w), energy 3.5 kcal g⁻¹) and water ad libitum under a 12 h light/12 h dark cycle (lights on at 7:30 a.m.) and constant temperature (21–22 °C) and humidity (55 ± 5%).

Faecal samples from randomly selected subsets of three patients with AN (females aged 20, 22 and 20 yr with BMI of 10.3, 11.5 and 11.7 kg m⁻², respectively) and three HC (females aged 20, 22 and 21 yr with BMI of 22.6, 21.1 and 21.2 kg m⁻², respectively) who were representatives of cases and controls were used to colonize 6-week-old female GF littermates. Briefly, 250 mg of faecal samples were suspended with 5 ml of LYBHI media (supplemented with 0.05% Cysteine and 0.2% Hemin as reducing agents) diluted in 20% glycerol (20 ml g⁻¹ of faeces) in an anaerobic cryovial; these inoculum samples were then vortexed for 5 min, followed by 5 min standing to precipitate particles. The faecal slurries were then transferred to 1 ml cryovials and immediately frozen at −80 °C. At day one, both groups of mice were housed in autoclaved individually ventilated cages where they received the first dose (200 µl) of faecal slurries. The mice were then given autoclaved chow diet and water ad libitum for 2 d and their food intake was recorded. At day 3, mice were gavaged with a second dose of faecal material from the same matched AN and HC donors as before. Thereafter, mice in both groups were singly housed and subjected to 30% energy-restricted autoclaved chow diet (amount of given food was set at 70% of ad libitum food intake for each mouse) for 3 weeks where water was given ad libitum. Both the anorexia-transplanted (AN-T) and the normal control-transplanted (HC-T) mice were weighed every 5 d after the start of energy-restricted dieting.

At the end of the study, mice were anaesthetized with isoflurane and blood from the vena cava was collected in tubes containing EDTA. Blood samples were centrifuged for 6 min at 4,032 g at 4 °C. Plasma was isolated and stored at −80 °C for subsequent biochemical testing. Inguinal subcutaneous white adipose tissue, the caecum and the hypothalamus of each mouse were precisely dissected and collected for quantitative PCR analysis. No animal or data points were excluded from the present study.

Total RNA was extracted from tissues using Trizol reagent (Invitrogen) according to the manufacturer’s instructions, followed by concentration measurement. One µg of RNA was transcribed to complementary DNA using the Reverse Transcription System (Promega). Real-time PCR was performed using the LC480 detection system (Roche Diagnostics) and SYBR Green I Supermix (Takara). All qPCRs were run on the thermal cycles at 95 °C for 10 min, followed by 45 cycles of 0.01 s at 95 °C and of 20 s at 60 °C. Data were normalized to the housekeeping Rpl36 gene for adipose tissue and Rplp0 gene⁷⁹ for hypothalamus and analysed according to the delta-delta CT method. Sequences of oligonucleotides used in this study are provided in Supplementary Table 8.

DNA extraction, 16S rRNA sequencing and data analyses in mice experiments

Microbial DNAs were isolated and purified from stool samples (~250 mg) of human donors and mouse recipients by using NucleoSpin soil mini kit (MACHEREY‑NAGEL). The DNA was then amplified using Phusion High-Fidelity PCR master mix (New England Biolabs) by PCR targeting the V3-V4 region of the 16S rRNA gene (primer sequences provided in Supplementary Table 8). The following PCR programme was used: 98 °C for 30 s, 25× (98 °C for 10 s, 55 °C for 20 s, 72 °C for 20 s), 72 °C for 5 min. Amplification was verified by running the products on an agarose gel. Indices were added in a subsequent PCR using an Illumina Nextera kit with the following PCR programme: 98 °C for 30 s, 8× (98 °C for 10 s, 55 °C for 20 s, 72 °C for 20 s), 72 °C for 5 min. Attachment of indices was verified by running the products on an agarose gel. Products from the nested PCR were pooled on the basis of band intensity and the resulting library cleaned with magnetic beads. The DNA concentration of pooled libraries was measured fluorometrically. Sequencing was done on an Illumina MiSeq desktop sequencer using the MiSeq reagent kit V3 (Illumina) for 2 × 300 bp paired-end sequencing. Paired-end reads were subsequently trimmed, merged and analysed using the DADA2 (v1.16.0) pipeline⁸⁰.

Statistical analysis

No data were excluded before the statistical analysis in the present study. No allocation and randomization were included as the study is observational. This study includes all available samples (n = 147) of patients with anorexia and healthy individuals. Although no statistical methods were used to predetermine sample sizes, our sample sizes are similar to those reported in previous publications^43,81. Samples were randomly distributed across metagenomics and metabolomics batches. Investigators were blinded to group allocation during data collection in metagenomic, biochemical and metabolomics analyses. All analyses of human samples were performed using R (v4.1.2). Gene expression analyses and body weight comparisons for animal studies were performed using GraphPad Prism (v9.3.0).

1. (1)

   Differential analysis

   We carried out the differential analysis using the metadeconfoundR pipeline implemented in the R package metadeconfoundR (v0.1.8; see https://github.com/TillBirkner/ metadeconfoundR or https://doi.org/10.5281/zenodo.4721078) where we assessed the extent to which the observed differences between AN and HC participants in microbiome or metabolome analyses are confounded by covariates including age, BMI, smoking and medication. This pipeline initially used univariate statistics to find associations between microbiome features and disease status, followed by nested linear model comparison post hoc testing to check for the confounding effects of potential covariates and finally, returning a status label.

2. (2)

   Association analysis

   For the association and mediation analysis between omics features and eating behaviour and psychological traits within the AN group, we first checked the normality of continuous variables with Shapiro-Wilk normality test, finding most variables not to be normally distributed. Therefore, before association analysis, we standardized the continuous variables using empirical normal quantile transformation to follow a standard normal distribution (N ~ (0, 1)). Then we implemented a linear regression model to assess the associations between omics features and eating behaviour and psychological traits using the following formula where confounding factors were added as covariates.

   $$\begin{array}{l}{\mathrm{Eating}}\,{\mathrm{behavior}} – {\mathrm{and}}\,{\mathrm{psychological}}\,{\mathrm{traits}}\sim {\mathrm{omics}}\,{\mathrm{features}} \\ ( {\mathrm{for}}\,{\mathrm{example,}}\, {\mathrm{metabolite/msp/SV}}) \\ + {\mathrm{Age}} + {\mathrm{BMI}} + {\mathrm{Smoking}} + {\mathrm{Medication}}\end{array}$$

   Medication included selective serotonin re-uptake inhibitors, antipsychotics and benzodiazepines.

3. (3)

   For the association analysis between omics features and host metabolic traits, a normality check and data standardization were also performed before the linear regression analysis. In the linear regression model, the above-mentioned confounding factors were included as covariates, except for BMI as the extremely low BMI is the most remarkable phenotypic change for AN patients compared with HC individuals.

Differences in gut microbial diversity (gene richness, species count, taxonomic composition) between AN and HC were calculated using Wilcoxon tests, and Kruskal-Wallis test was used for assessing the significance of differences between multiple groups. Unless otherwise stated, all P values were corrected using the Benjamini-Hochberg method and P_adj < 0.05 was considered statistically significant.

Reporting summary

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