The human gut microbiome in early-onset type 1 diabetes from the TEDDY study

Type 1 diabetes (T1D) is an autoimmune disease that targets pancreatic islet beta cells and incorporates genetic and environmental factors1, including complex genetic elements2, patient exposures3 and the gut microbiome4. Viral infections5 and broader gut dysbioses6 have been identified as potential causes or contributing factors; however, human studies have not yet identified microbial compositional or functional triggers that are predictive of islet autoimmunity or T1D. Here we analyse 10,913 metagenomes in stool samples from 783 mostly white, non-Hispanic children. The samples were collected monthly from three months of age until the clinical end point (islet autoimmunity or T1D) in the The Environmental Determinants of Diabetes in the Young (TEDDY) study, to characterize the natural history of the early gut microbiome in connection to islet autoimmunity, T1D diagnosis, and other common early life events such as antibiotic treatments and probiotics. The microbiomes of control children contained more genes that were related to fermentation and the biosynthesis of short-chain fatty acids, but these were not consistently associated with particular taxa across geographically diverse clinical centres, suggesting that microbial factors associated with T1D are taxonomically diffuse but functionally more coherent. When we investigated the broader establishment and development of the infant microbiome, both taxonomic and functional profiles were dynamic and highly individualized, and dominated in the first year of life by one of three largely exclusive Bifidobacterium species (B. bifidum, B. breve or B. longum) or by the phylum Proteobacteria. In particular, the strain-specific carriage of genes for the utilization of human milk oligosaccharide within a subset of B. longum was present specifically in breast-fed infants. These analyses of TEDDY gut metagenomes provide, to our knowledge, the largest and most detailed longitudinal functional profile of the developing gut microbiome in relation to islet autoimmunity, T1D and other early childhood events. Together with existing evidence from human cohorts7,8 and a T1D mouse model9, these data support the protective effects of short-chain fatty acids in early-onset human T1D.

species, and deficiency of bacteria that produce short-chain fatty acids (SCFAs) 7,8 in cases of T1D or islet autoimmunity (IA) 8,11,15,18 . Corroborating these findings, decreased levels of SCFA-producing bacteria were found in adults with type 2 diabetes (T2D) 19 . In addition, increased intestinal permeability 14 and decreased microbial diversity 12 after IA but before T1D diagnosis have been reported. Studies using the nonobese diabetic (NOD) mouse model have determined immune mechanisms that mediate the protective effects of SCFAs 9 and the microbiome-linked sex bias in autoimmunity 20 . NOD mice fed specialized diets resulting in high bacterial release of the SCFAs acetate and butyrate were almost completely protected from T1D 9 . A study in a streptozotocin-induced T1D mouse model demonstrated that bacterial products recognized in pancreatic lymph nodes contribute to pathogenesis 21 .
Even in the absence of immune perturbation, the first few weeks, months and years of life represent a unique human microbial environment that has only recently been detailed 22,23 . Infants have a markedly different gut microbial profile from adults, characterized by a distinct taxonomic profile, greater proportion of aerobic energy harvest metabo lism, and more extreme dynamic change 24 . These differences gradually fade over the first few years of life, particularly in response to the introduction of solid food, and individual microbial developmental trajectories are influenced by environment, delivery mode, breast (versus formula) feeding, and antibiotics [25][26][27] . Most studies that address the development of the gut microbiome, both generally and in association with T1D, have used gene analysis of 16S rRNA, which leaves open the question of functional and strain-specific differences that are not easily detected by this technology that might contribute to disease pathogenesis 12 .
Bridging this gap is one goal of the The Environmental Determinants of Diabetes in the Young (TEDDY) study, a prospective study that aims to identify environmental causes of T1D 28 . It includes six clinical research centres in the United States (Colorado, Georgia/Florida and Washington) and Europe (Finland, Germany and Sweden), which together have recruited several thousand newborns with a genetic predisposition for T1D or first-degree relative(s) with T1D. This has enabled the TEDDY study to collect a range of biospecimens, including monthly stool samples starting at three months of age, coupled with extensive clinical and personal data such as diet, illnesses, medications and other life experiences. To characterize microbial, environmental, genetic, immunological and additional contributors to the development Letter reSeArCH of T1D, the TEDDY study group further assembled nested case-control studies for IA (n = 418 case-control pairs) and T1D (n = 114) 29 . Casecontrol pairs were matched by clinical centre, sex and family history of T1D, which are all known confounding factors for T1D susceptibility and microbiome composition.
Here, we assessed 783 children followed from three months to up to five years of age from six clinical centres in four countries (Finland, Germany, Sweden and the United States) who either progressed to persistent IA or T1D or were matched as controls (Fig. 1a, b, Extended Data Table 1). Stool samples were collected, on average, monthly starting at three months of age and continuing until the clinical end point (IA or T1D). This study focused solely on analysing metagenomic sequencing data (n = 10,903 samples, n = 783 subjects), while a companion paper by Stewart et al. 30 interrogated corresponding 16S rRNA amplicon sequencing information.
We first investigated the taxonomic composition of early gut metagenomes at the species level. Principal coordinate analysis ordination of Bray-Curtis beta diversities showed a strong longitudinal gradient and marked heterogeneity among the earliest samples (Fig. 2a, Extended Data Fig. 1a-k, Supplementary Note 1). Permutational analysis of variance (ANOVA) of Bray-Curtis beta diversities indicated that inter-subject differences explained 35% of microbial taxonomic variation (permutation test, P < 0.001, 1,000 permutations), followed by age at stool sampling at roughly 4% of variance (P < 0.001). Using cross-sectional analysis to test for associations between taxonomic beta diversities and other collected metadata, we found that in addition to subject ID and age, geographical location and breastfeeding had strong and systematic effects on the composition of the microbial community (Supplementary Table 1, Extended Data Fig. 2a-d, Supplementary Note 1). To investigate the stability and individuality of the microbial profiles further, we compared intra-and intersubject Bray-Curtis beta diversities. The gap between individual stability and similarity within or across clinical centres was largest at the beginning of the sampling period, indicating that the children had particularly dissimilar micro biota during these early months (Fig. 2b, Supplementary Note 1). Finally, we tested microbial alpha diversity (Shannon's diversity index) of taxonomic profiles for associations with collected metadata, and found that the cessation of breastfeeding had the largest effect (ANOVA, partial η 2 = 0.053) in the accrual of alpha diversity in early life (Supplementary Table 2 Table 2, Supplementary Note 2). Previous studies have found Bifidobacterium species to be especially vulnerable to antibiotics 31,32 , leading us to investigate how antibiotic perturbations influenced these common dominant members of the early gut. Comparing microbial relative abundances before and after antibiotics (assuming that the given species was present in the preceding sample), we saw a decrease in the abundances of the Bifidobacterium members B. bifidum, B. pseudocatenulatum, B. adolescentis, B. dentium and B. catenulatum, whereas B. longum and B. breve did not systematically decline owing to antibiotics (Fig. 2d), suggesting that certain Bifidobacterium species are particularly susceptible to out-competition by other community members after depletion by antibiotics. Given their dominance in the typical developing gut microbiota and finely tuned balance of metabolic interactions with breast milk, this finding underscores the importance of approaching antibiotic prescriptions in early childhood with care, especially during breastfeeding.
Accompanying our taxonomic profiling, functional profiling of these metagenomes suggested the development of a consistent microbial functional core during infancy, with a smaller subject-specific variable functional pool (Extended Data Fig. 5a, b, Supplementary Note 3). As in most microbial community studies 33 , microbial gene families of uncharacterized function made up a substantial fraction of these profiles, averaging roughly 50% based on Gene Ontology 34 annotations (Extended Data Fig. 5c) and more than 90% based on more functionally specific MetaCyc pathways (Extended Data Fig. 5d). We observed an increasing longitudinal trend in the proportion of unmapped reads (Extended Data Fig. 5e, Pearson's r = 0.318, P < 2.2 × 10 −16 ). However, within the reads that mapped to either microbial pangenomes or known protein sequences (the proportion of which decreased with age), we saw an increase in the proportion of reads with MetaCyc annotation, mainly during the first year (Extended Data Fig. 5f, Pearson r = 0.391, P < 2.2 × 10 −16 ). This suggests that although the early life microbiome is relatively well-covered by current microbial reference genomes, less functional and biochemical characterization has been carried out on gene families within these microorganisms, which will thus particularly benefit from future work. In addition to broadly conserved and subject-specific functions, we identified a range of microbial metabolic enzymes that consistently increased or decreased in abundance over the first year of life, paralleling shifts in community structure and infant diet (Fig. 3, Supplementary Note 3, Supplementary Table 3). For example, the enzyme l-lactate dehydrogenase (1.1.1.27), which is well-characterized in Bifidobacteria for its role in milk fermentation 35 , was among the most consistently declining enzymes over this period, notably coinciding with the cessation of breastfeeding in many infants (from 73% breastfed at month 3 to 28% at year 1). Conversely, the enzyme transketolase (2.2.1.1), which has been implicated previously 36 in the metabolism of fibre, was among the most consistently increasing enzymes, which also coincided with increased incorporation of solid food (a component of 53% of infants' diets at month 3 versus 100% at year 1). Hence, these notable changes in community functional potential highlight the unique metabolic environment of the early infant gut, and the subsequent transition to a more adult-like gut microbiome that is adapted to variable, fermentative energy sources.
Combining taxonomic and functional profiles to test for differences between cases and controls, we used linear mixed-effects modelling and identified a relatively small number of individual taxonomic and functional features that were associated with case-control outcome (Supplementary Table 4), most with borderline statistical significance (false discovery rate (FDR) corrected q-values indicated below). We confirmed separation between cases and controls by random forest classifiers (Extended Data Fig. 6a, b, Supplementary Note 4). In the IA case-control cohort, healthy controls contained higher levels of Lactobacillus rhamnosus (q = 0.055), supporting protection against IA by early probiotic supplementation 37 (Extended Data Fig. 6c, d, Supplementary Note 5). IA controls also had more Bifidobacterium dentium (q = 0.054), whereas IA cases had on average higher abundance of Streptococcus group mitis/oralis/pneumoniae species (q = 0.11). In T1D case-control comparisons, controls had higher levels of Streptococcus thermophilus (q = 0.078) and Lactococcus lactis (q = 0.094) species, both common in dairy products, whereas cases contained higher levels of species such as Bifidobacterium pseudocatenulatum (q = 0.078), Roseburia hominis (q = 0.11) and Alistipes shahii (q = 0.14). Even though our modelling approach controlled for regional differences in clinical centres, we found additional but often weak associations with outcome in some clinical centres when tested separately (Supplementary Table 4). Finnish IA cases had more Streptococcus group mitis/oralis/pneumoniae species (q = 0.0008), IA controls from Colorado had more Streptococcus thermophilus (q = 0.0059), and Swedish IA cases contained more Bacteroides vulgatus (q = 0.090).
Pathways with the highest statistical significance in case-control comparisons were related to bacterial fermentation (Supplementary Table 4). The superpathway of fermentation (MetaCyc identifier PWY4LZ-257) was increased in controls in the T1D cohort (q = 0.019) and Finnish IA cohort (q = 0.049). SCFAs such as butyrate, acetate and propionate are common by-products of bacterial fermentation, and butyrate and acetate protected NOD mice against T1D 9 . Consistently, we observed that several bacterial pathways that contribute to the biosynthesis of short-chain fatty acids were increased in healthy controls.
Among pathways involved in butyrate production, the degradation of l-arginine, putrescine and 4-aminobutanoate (ARGDEG-PWY) superpathway was increased in T1D controls cohort-wide (q = 0.043), whereas the fermentation of acetyl coenzyme A to butanoate (PWY-5676) was more abundant in the Finnish T1D controls (q = 0.053). The degradation of acetylene (P161-PWY), which contributes to acetate production, was increased in T1D controls cohort-wide (q = 0.14), and the degradation of l-1,2-propanediol (PWY-7013), which is involved in propionate biosynthesis, was higher in the German T1D controls (q = 0.019). These findings support existing evidence for the protective effects of SCFAs in human T1D 7,8 and T2D 19 cohorts and the NOD mouse model 9 .
As reflected by the community-level analyses, human milk with its pro-and prebiotic functions is one of the main factors that determine the community composition of the infant gut microbiome. Bifidobacterium longum subsp. infantis is a particularly versatile degrader of human milk oligosaccharide (HMO) that is often found in stool samples collected during breastfeeding 38 . By following the families representing genes in the B. longum subsp. infantis HMO gene cluster 39,40 in our data, we found that an additional 30 bacterial species  Gut microbial communities were highly individual. c, Influence of antibiotic (Abx) courses on microbial stability, measured by Bray-Curtis dissimilarity over consecutive stool samples (<50 days apart) from the same individual during the first three years of life, and stratified by whether antibiotics were given between the two samples (n = 654 observations with antibiotics, n = 6,734 observations without antibiotics). Curves show locally weighted scatterplot smoothing (LOESS) for the data per category. Shaded areas show permutation-based 95% confidence intervals for the fit. d, Decreases in the most common Bifidobacterium species in connection to oral antibiotic treatments. Fold change was measured between consecutive samples with an antibiotic course between them, given that the species in question was present in the first of the two samples. Sample size per species (n) indicates the number of sample pairs in which the species in question was present in the sample before the antibiotic treatment. Bars show bootstrapped mean log 2 (fold change) (that is, decrease), and error bars denote s.d. (n = 1,000 bootstrap samples).
Letter reSeArCH carried at least one homologue with more than 50% sequence identity to one or more HMO utilization genes (Supplementary Table 5). As expected, many Bifidobacteria carried several homologues, but surprisingly three Enterococcus species (E. casseliflavus, E. faecalis and E. faecium) also carried seven or more homologues (Supplementary Table 5).
To identify strain-level adaptation similar to B. longum subsp. infantis, we further examined whether any of these genes showed contrasting prevalence between samples collected during breastfeeding and after weaning, given that the carrier species itself was present. In total, 41 gene families were observed more often during breastfeeding (Supplementary Table 5, test of proportions, adjusted P < 0.001); most (37 out of 41) were carried by B. longum (Fig. 4), and B. pseudocatenulatum contained four such gene families (Extended Data Fig. 7, Supplementary Table 5). In samples with B. longum, this implicated a clear strain shift after weaning, when fewer B. longum strains carried these genes (Fig. 4). In samples with B. pseudocatenulatum, four gene families showed a similar but less contrasting pattern (Extended Data Fig. 7). Overall, these observations identify new candidate species that contribute to HMO processing or exploitation, and link strain composition to specific driving molecular functions that potentially explain selective sweeps during microbiome development, in this case specifically related to breastfeeding.
Despite ample sample size, scrutiny of the study design, and thorough statistical analyses, most of the taxonomic and functional signals we detected in case-control comparisons were modest in effect size and statistical significance. This could be due to several reasonsdifferences between T1D endotypes, temporally diffuse signals, geographical heterogeneity, or lack of stool samples for the first two months of life -and these should be considered in future investigations (Supplementary Note 6). Furthermore, the data used in these investigations was composed of samples from the genetically predisposed and mostly white, non-Hispanic case-control groups designed into the TEDDY study. Results cannot be guaranteed to reflect the whole TEDDY cohort or child populations in the respective countries.
Future targeted approaches to identify subject-specific connections between the gut microbiota and T1D pathogenesis may be beneficial, particularly given the apparent population-level heterogeneity revealed here. For example, laboratory experiments involving dietary factors that have been associated with the onset of T1D 3 may reveal biochemically specific signals that are mediated by the microbiome. Different endotypes of disease, such as differences in the first appearing autoantibody (IAA versus GADA), the number of appearing autoantibodies, the time from seroconversion to T1D diagnosis, genetic host risk alleles and ethnic backgrounds, may be characterized by distinct microbial configurations (Supplementary Note 6). Finally, components of the microbiome that were poorly measured in these data may also have crucial roles: viruses, fungi, microbial transcription or small-molecule biochemistry. By surveying these additional molecular activities by cross-sectional analysis and in more detailed longitudinal populations, this study lays

Competing interests
The authors declare no competing interests. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http:// creativecommons.org/licenses/by/4.0/.

MEthOdS
Cohort and study design. TEDDY is a prospective cohort study funded by the National Institutes of Health with the primary goal to identify environmental causes of T1D. It includes six clinical research centres-three in the United States (Colorado, Georgia/Florida, Washington) and three in Europe (Finland, Germany and Sweden). Detailed study design and methods have been previously published 28,41,42 . Written informed consents were obtained for all study participants from a parent or primary caretaker, separately, for genetic screening and participation in a prospective follow-up. The This analysis used stool samples and clinical metadata from two nested casecontrol studies (persistent, confirmed IA or T1D) using risk set sampling 29 . The data used here were collected as of 31 May 2012, as a 1:1 match in which one control per case of persistent confirmed IA or T1D was selected from the full TEDDY cohort. A control was a participant who had not developed persistent, confirmed IA or T1D by the time the case to which it was matched had developed IA or T1D, within ±45 days of the event time. Matching factors were clinical centre, sex and family history of T1D to control for differences in geographical area, genetic background and in sample or data handling between clinical centres. In all case-control comparisons, we removed all case-control pairs in which the control later progressed to case status (that is, progressed to IA or T1D). In additional, 17 subjects with missing information about breastfeeding together with their matched pairs were excluded from the case-control comparisons to avoid confounding effects from unknown breastfeeding status.
The development of persistent, confirmed IA was assessed every three months. Persistent autoimmunity was defined by the presence of confirmed islet autoantibody on two or more consecutive visits. The date of persistent autoimmunity was defined as the draw date of the first sample of the two consecutive samples that deemed the child persistently positive for a specific autoantibody (or any autoantibody). T1D was defined according to American Diabetes Association criteria for diagnosis 43 .
Stool samples were collected monthly starting at three months of age and continuing up until 48 months of age, then every three months until 10 years of age and then biannually thereafter, into the three plastic stool containers provided by the clinical centre. Children who were antibody negative after 4 years of age were encouraged to submit four times a year even though after 4 years their visits schedule switched to biannual. Parents sent the stool containers at either ambient or +4 °C temperature with guaranteed delivery within 24 h in the appropriate shipping box to the NIDDK repository if living in the United States or their affiliated clinical centre if living in Europe. The European clinical centres stored the stool samples and sent monthly bulk shipments of frozen stool to the NIDDK repository. The TEDDY Manual of Operations, including the stool sample collection protocol, can be accessed online at https://repository.niddk.nih.gov/static/studies/ teddy/teddy_moop.pdf.
A priori power calculations using discrete Cox's proportional hazards regression 44 for the matched IA case-control study estimated 80% power, α = 0.01, twosided test to detect an odds ratio > 3 for an exposure with 5% prevalence, to an odds ratio > 1.8 for an exposure with 20% prevalence. The experiments were not randomized, and investigators were not blinded to allocation during experiments and outcome assessment. Metagenomic sequencing and initial bioinformatics. Samples were metagenomically sequenced as one library each multiplexed through Illumina HiSeq machines using the 2 × 100-bp paired-end read protocol. Samples with limited DNA quantity and/or too few high-quality reads were filtered out, resulting in a discrepancy of sample frequencies between the metagenomic data and the 16S rRNA amplicon sequencing data analysed in the companion paper 30 . Casava v1.8.2 (Illumina) output initial FASTQ files from the resulting data were processed using cutadapt v1.9dev2 for adaptor removal, Trim Galore v0.2.8 (Babraham Bioinformatics) for removing low-quality bases and PRINSEQ v0.20.3 45 for sample demultiplexing. Bowtie2 v2.2.3 was used to map reads to the human genome for decontamination before subsequent analysis.

Taxonomic and functional profiling by MetaPhlAn2 and HUMAnN2.
Taxonomic profiling of the metagenomic samples was performed using MetaPhlAn2 46 v2.6.0, which uses a library of clade-specific markers to provide pan-microbial (bacterial, archaeal, viral and eukaryotic) quantification at the species level. MetaPhlAn2 was run using default settings.
Functional profiling was performed with HUMAnN2 47 v0.9.4. For an input metagenome, HUMAnN2 constructs a sample-specific reference database by concatenating and indexing the pangenomes of species detected in the sample by MetaPhlAn2 (pangenomes are pre-clustered, pre-annotated catalogues of open reading frames found across isolate genomes from a given species 48 ). HUMAnN2 then maps sample reads against this database to quantify gene presence and abundance in a species-stratified manner, with unmapped reads further used in a translated search against UniRef90 49 to include taxonomically unclassified but functionally distinct gene family abundances. Finally, for community-total, species-stratified, and unclassified gene family abundance, HUMAnN2 reconstructs metabolic pathway abundance based on the subset of gene families annotated to metabolic reactions (based on reaction and pathway definitions from MetaCyc 50 ). Enzyme (level-4 Enzyme Commission (EC) categories) abundances were further computed by summing the abundances of individual gene families annotated to each EC number based on UniRef90-EC annotations from UniProt 51 . Phenotype and covariate analysis. This study includes extensive collection of clinical covariates that cover several aspects of common and rare life events in early childhood from infancy up to five years of age. In these analyses, we used information that is, according to the literature, of high relevance in terms of gut microbiome development. Information about mothers, pregnancy and birth was collected during the three-month clinic visit by questionnaire and included the mode of birth (vaginal birth versus caesarean section), gestational age, infant's 5-min Apgar score, information about maternal diabetes (T1D, T2D or gestational diabetes) and maternal insulin and medication use (antibiotics, angiotensinconverting enzyme inhibitors, metformin, glyburide, antihypertensives) during pregnancy. Dietary information used in these analyses includes the start (and end) date for the following dietary compounds: breastfeeding, baby formula, cow's milk, gluten, cereals, meat, vegetables and fruits. The start of solid food (anything other than breast milk or cow's milk) was also analysed separately. The T1D-associated autoantibodies IAA, GADA and IA2A were analysed from serum samples collected at each clinic visit. In addition to IA, defined as persistent, confirmed autoantibody seropositivity, we analysed the data in terms of the persistency and cumulative frequency of autoantibodies (single or multiple autoantibodies). In TEDDY, all prescribed antibiotic courses are recorded. We further stratified these data by the type of antibiotic in five categories: amoxicillin, penicillin, cephalosporins, macrolide and other antibiotics. Information about probiotics covered the dates for starting and stopping probiotic supplementation, but not the specific types of probiotics used. In addition, sex, information about whether first degree relatives in family had T1D, and HLA haplotypes of the subjects were used in these analyses. Subjects screened from the general population were identified with high-risk alleles (89%) including: DRB1*04-DQA1*03-DQB1*03:02/DRB1*03-DQA1*05-DQB1*02:01 (DR3/4), DRB1*04-DQA1*03-DQB1*03:02/DRB1*04-DQA1*03-DQB1*03:02 (DR4/4), DRB1*04-DQA1*03-DQB1*03:02/DRB1*08-DQA1*04-DQB1*04:02 (DR4/8) and DRB1*03-DQA1*05-DQB1*02:01/DRB1*03-DQA1*05-DQB1*02:01 (DR3/3), plus six genotypes specific to first-degree relatives 28 .
Principal coordinate analysis (PCoA) ordination was generated using t-distributed stochastic neighbour embedding (t-SNE) as implemented in Rtsne package in R with Bray-Curtis dissimilarity as the distance measure and perplexity (a free parameter) equal to 50. Statistical significance of the trends between early clusters and metadata were tested using mixed-effect logistic regression and samples collected during the first year of life as follows. The target variable used was a binary indicator of whether the relative abundance of the taxon of interest (three different Bifidobacterium species or phylum Proteobacteria) was greater than 0.5 (definition of the cluster). The age of sample collection, mode of delivery, clinical centre, breastfeeding status (ongoing/stopped), solid food status (binary variable indicating whether solid food was introduced in the diet) and antibiotics status (binary variable indicating whether the subject received antibiotics during the last 30 days) were used as fixed effects, and the subject ID was used as a random effect.
Associations between microbial feature abundances and clinical outcome were determined using MaAsLin 52 . In brief, this multivariate linear modelling system for microbial data selects from among a set of (potentially high-dimensional) covariates to associate with microbial taxon or pathway abundances. Mixed-effects linear models using a variance-stabilizing arcsin square root transform on relative abundances are then used to determine the significance of putative associations from among this reduced set. In the models, subject ID was used as a random effect, and the age of sample collection, mode of delivery, clinical centre (for cohort-wide comparisons), breastfeeding status (ongoing or stopped), solid food status (binary variable indicating whether solid food was introduced in the diet), number of sequencing reads and case-control outcome were used as fixed effects. Nominal Letter reSeArCH P values were adjusted using the Benjamini-Hochberg FDR method. Here, microbial features with corrected q < 0.25 were reported. For metabolic pathways, pseudocount 2 6 was added to CPM values to stabilize the variation in lowly abundant and/or prevalent but highly variable categories, and data were log 2 -transformed.
As previously described 40 , to associate microbial diversity with covariates while accounting for nonlinear, age-dependent effects, we first fitted a sigmoid function (nls function in R) to account for the longitudinal trend. Residuals of this model were then used as inputs for a mixed-effect model (glmmPQL function in the MASS R package), with subject IDs as random effects to account for repeated measurements in the data. Other factors were included in the model as fixed effects, and their significance levels were evaluated using P values reported by the model (Supplementary Table 2).
The association between T1D case-control outcome and microbial alpha diversity in individual clinical centres was tested using a linear mixed-effects model (glmmPQL function in MASS R package) on samples 730 days or less before T1D diagnosis. In the model, the age at stool sample collection and T1D case-control outcome were used as fixed effects, and subject ID was used as a random effect. Microbial variance explained by clinical and other covariates. Variance analysis was conducted using the adonis function in the vegan R package given a Bray-Curtis dissimilarity matrix of the taxonomic profiles and all TEDDY clinical metadata listed above. In brief, adonis conducts multivariate ANOVA using the dissimilarity matrix (that is, partitions the sums of squares) given the metadata as covariates. Statistical significance of the fit was assessed using permutation tests. HMO gene homology. The HMO gene cluster homologues between B. longum subsp. infantis and multiple taxa were analysed as follows. UniRef90 gene families corresponding to the protein sequences in the B. longum subsp. infantis HMO gene cluster 39 (protein sequences Blon_2331-Blon_2361 in NCBI protein sequence database) were identified by translated BLAST search against ChocoPhlAn pangenome collection 48 used by HUMAnN2. Identified hits were further filtered by requiring ≥50% alignment identity and ≥80% mutual coverage. Combining this information with HUMAnN2 species-stratified UniRef90 gene family quantification enabled calling these genes present given that they had sufficient read coverage, here defined as log 10 (counts per million) > 0.1 in at least 50 samples collected during breastfeeding. Differential gene prevalence during breastfeeding was tested using the samples in which the carrier species had >1% relative abundance. Testing was conducted using the test of equal or given proportions (prop.test function in R) and by comparing the prevalence (proportion of the samples for which the species in question harboured the gene according to the metagenomic data) of the gene in samples collected during breastfeeding with the samples collected after weaning. P values were adjusted for multiple testing by Benjamini-Hochberg method (p.adjust function in R). All homologues together with their BLAST search metrics, prevalence in the metagenomic data and corresponding B. infantis HMO gene are reported in Supplementary Table 5. Bacterial growth assays. Bifidobacterium bifidum strain RJX-1201, Bifidobacterium breve RJX-1202 and Bifidobacterium longum RJX-1203 were streaked on brain heart infusion agar (BD) supplemented with 1% vitamin K/hemin solution (BD; sBHI), and incubated for 48 h in a vinyl anaerobic chamber (Coy Laboratory Products) containing 5% CO 2 , 5% H 2 and 90% N 2 and maintained at 37 °C. Cells were transferred to sBHI liquid medium (BHI broth, BD, supplemented as above) and grown for 24 h in anaerobic conditions. Cultures were washed twice with PBS and optical density at 600 nm (OD 600 ) was measured using a BioTek PowerWave 340 plate reader. OD 600 was normalized to 0.2 for all strains and 5 μl bacteria inoculum was added to a final volume of 200 μl containing 10% sBHI and 125 mM carbon source (glucose, fucose, galactose or lactose) in a 96-well plate. OD 600 was measured in the plate reader every hour for 48 h with 5 s of medium shaking before each measurement. All of the measurements were normalized to a medium-only blank. Experiment was repeated three times (n = 3) in triplicate and one representative experiment is shown. Error bars are s.d. of three technical replicates. Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this paper. Code availability. Code for Random Forest case-control comparisons and cohort wide MaAsLin association analyses in Supplementary Table 4 has been made publicly available at https://github.com/tvatanen/broad_teddy_microbiome_analyses. Other analysis software including quality control, taxonomic, and functional profilers is publicly available and referenced as appropriate.

Data availability
TEDDY microbiome 16S and whole-genome sequencing data that support the findings of this study are available in the NCBI database of Genotypes and Phenotypes (dbGaP) with the primary accession code phs001443.v1.p1, in accordance with the dbGaP controlled-access authorization process. Clinical metadata analysed during the current study are available in the NIDDK Central Repository at https://www.niddkrepository.org/studies/teddy. The curve shows the generalized additive model (GAM) fit for the data, and the shaded area shows the 95% confidence interval for each fit, as implemented in geom_smooth function in ggplot2 R package. b, Shannon's diversity for the samples in the IA case-control cohort (n = 7,051) with respect to the time to the appearance of first autoantibody (seroconversion). The curves show LOESS fits for cases and controls separately, and the shaded area shows 95% confidence intervals for each fit. c, Shannon's diversity for the samples in the T1D case-control cohort (n = 3,309) with respect to the time to T1D diagnosis. The curves and shaded areas are as in b. d, As in c, but only for data (n = 983 samples) for subjects in Finland. No difference between cases and controls. e, As in c, but only for data (n = 142 samples, n = 6 subjects) for subjects in Georgia, USA. Cases show a drop in alpha diversity before the diagnosis of T1D (linear mixed-effects model, P = 0.0033).    Although the baseline level of functional similarity is significantly greater than that of taxa (see Fig. 2b), functional states and development trajectories also both retain a level of personalization. The stability of the functional profiles was evaluated in three-month time windows, over two-month increments. Lines show the median per time window, and shaded area denotes the 99% confidence interval estimated using binomial distribution. c, d, Proportion of metagenomic gene abundance with functional annotation through Gene Ontology (c) and MetaCyc (d) databases. The metagenomic reads were divided into the following categories: reads that could be mapped to genes with functional assignment in the database in question (annotated), and reads with no annotation but alignment to species pangenomes or UniProt proteins (unannotated). The proportion of the unknown genes (unmapped) was estimated using the number of reads with unknown origin. e, The proportion of unmapped reads, reflecting the relative abundances of reads not mappable to any microbial pangenomes in the available reference set or to UniProt. An increasing trend of unmapped reads with respect to the age at sample collection continued through approximately two years of age. f, The proportion of reads with confident functional annotation in MetaCyc within the genes that mapped to species pangenomes or UniProt proteins. The data again showed an increasing longitudinal trend, implicating a deficit of functional and biochemical annotations within microorganisms that are abundant during the first year of life.  Fig. 6 | Differences between cases and controls. a, The gut microbiome functional (left) and taxonomic (right) profiles were classified between cases and controls using leave-one-out cross-validation (n = 3,366 samples), in which one case-control pair was held-out in turn. Data show error rates for classifying these held-out samples per fold (a data point per fold, n = 100 folds). This suggests weak but betterthan-random classification between cases and controls. Notched box plots are as in Extended Data Fig. 4. b, Average longitudinal abundance of Ruminococcus gnavus in Finland (n = 2,630 samples) stratified by the number of observed persistent autoantibodies (AABs); no autoantibodies (that is, healthy control), a single autoantibody, or multiple (two or more) autoantibodies. c, Average longitudinal abundance of Lactobacillus rhamnosus in IA cases and controls (n = 7,017 samples). L. rhamnosus is more abundant in controls (q = 0.055). The curves in b and c show LOESS fit per group, and shaded areas show 95% confidence interval for each fit, as implemented in geom_smooth function in ggplot2 R package. d, Abundance (left) and prevalence (right) of Lactobacillus reuteri and L. rhamnosus in the first stool sample of each individual (collected at approximately three months of age) in association with early probiotic supplementation. 'No probiotic' indicates no probiotics given before the first stool sample (n = 583); 'later probiotic' refers to probiotics given later than the first four weeks but before the first stool sample (n = 45); 'early probiotic' refers to probiotics given during the first four weeks of life (n = 84). Numbers (n) per clinical centre are given in Extended Data Table 2. L. reuteri and L. rhamnosus were more abundant and prevalent in groups with probiotics supplementation. Visual jitter was added to make data equal to zero distinguishable, and boxes denote the IQR, when applicable. The shown P values were obtained by applying Fisher's exact test (two-sided) to presence or absence count data (counting samples in which the species were present). Data on subjects' ethnic background were not systematically collected in European clinical centres but these study populations were predominantly white, non-Hispanic. Reported antihypertensive drugs were atenolol (n = 2), bisoprolol (n = 1), labetalol (n = 6), methyldopa (n = 1), methyldopa plus methyldopate (n = 3), metoprolol (n = 4) and nifedipine (n = 5). No use of angiotensin-converting enzyme (ACE) inhibitors was reported. Numbers indicate the number of subjects (n) if not specified otherwise.

Letter reSeArCH
Extended data

Statistical parameters
When statistical analyses are reported, confirm that the following items are present in the relevant location (e.g. figure legend, table legend, main text, or Methods section).

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement An indication of whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistics including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers upon request. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability TEDDY Microbiome 16S and WGS data that support the findings of this study are available in NCBI's database of Genotypes and Phenotypes (dbGaP) with the