Milk samples from 518 healthy mothers (one individual sample per mother) were obtained in 16 different countries, including cohorts from Equatorial Guinea (n = 33), Kenya (n = 30), Senegal (n = 60), South Sudan (n = 53), Tanzania (n = 4), Ecuador (n = 35), El Salvador (n = 6), Mexico (n = 18), Peru (n = 8), Puerto Rico (n = 5), mainland United States (n = 85), Austria (n = 38), Germany (n = 32), The Netherlands (n = 15), Norway (n = 20), and Spain (n = 76). The subject’s age, body weight, sampling time, and birth mode of the baby were listed in Table 1.

The study procedures related to the samples obtained from Equatorial Guinea, Kenya, Senegal, South Sudan, Ecuador, El Salvador, Mexico, Peru, Austria, Germany, The Netherlands, Norway, and Spain were approved by the overarching Institutional Review Board of the European Commission in the frame of the EU project “Variations in biochemical and microbiological milk composition among highly diverse human populations and their impact on infant gut ecosystem” (call FP-7-PEOPLE-2013-IEF) (protocol #624773, approved on 14 February 2014) and at each study location, and consent was obtained from each participating woman. Human milk samples from Puerto Rico were approved by the University of Puerto Rico, Medical Sciences Campus (IRB ProB2310120, approved on 24 March 2021). The United States samples were approved by the Rutgers University-New Brunswick Health Sciences Institutional Review Board (Protocol #Pro2018002781, approved on 17 September 2021; Protocol #Pro2020002169, approved on 9 June 2021), and consent was obtained from each participating woman. The samples from Tanzania are from a previous study approved by the IRB reference number 164-12-21052012 (46).

Collection of human milk was performed by the mothers as described by McGuire et al. (47); briefly, the aureola skin was wiped with antiseptic wipes containing chlorhexidine digluconate (bactiseptic wipes, Vesismin Health, Barcelona, Spain) using gloved hands, and milk was manually expressed into disposable sterile containers. Samples were shipped on dry ice to the Complutense University of Madrid (Spain), where they were stored at −80°C until an aliquot was sent on dry ice to Rutgers University (USA) for analyses. Human milk samples from Puerto Rico, the U.S. mainland, and Tanzania were collected either by hand expression or using the mother’s own sterile pump. Samples were frozen at the participants’ homes and transported on ice to the laboratory, where they were stored at −80°C until analysis.

The samples (1 mL) were centrifuged (15,000×g for 10 min at 4°C), and the fat layer was removed using a sterile swab. This step was repeated twice more to remove all fat. Then, the pellets together with a 200 μL fraction of the supernatants, were used for total DNA extraction employing the Dneasy PowerSoil Pro Kit (QIAGEN, Hilden, Germany), following the manufacturer’s instructions.

Primers 515 F (5’-GTGYCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACNVGGGTWTCTAAT-3′) were used to amplify the V4 hypervariable region of the bacterial 16S rRNA gene following the protocols for the Earth Microbiome Project.¹ The concentration of the pooled, purified and barcoded DNA amplicons was determined using the Qubit dsDNA HS assay kit (Thermo Fisher, Waltham, MA, USA). Amplicons were sequenced at Genewiz, LLC. (South Plainfield, NJ, USA) using the Illumina MiSeq platform (Illumina, CA, USA) with the Illumina MiSeq 2 × 150 bp paired-end protocol (Illumina Inc., San Diego, CA, USA) using the Illumina MiSeq platform.

For fungi, the Internal Transcribed Spacer 2 (ITS2) region was amplified from DNA obtained from extracted milk samples with a single ITS3 forward primer (5′-AATGATACGGCGACCACCGAG ATCTACACTATGGTAATTGTGCATCGATGAAGAACGCAGC-3′) and a barcoded ITS4 reverse primer (5′-CAAGCAGAAGACGGCATA CGAGATTCCCTTGTCTCCAGTCAGTCAGCCTCCTCCGCTTAT TGATATGC-3′) (38). Primers included Illumina adapters and the reverse primer included a unique 12 base Golay barcode (XX) for pooled demultiplexing (5′-CAAGCAGAAGACGGCATACGAG ATXXXXXXXXXXXXCGGCTGCGTTCTTCATCGATGC-3′). To determine PCR cycle counts appropriate for the fraction of fungal material, qPCR was performed on 518 extracted human milk samples. For samples with a Ct value less than 30, 25 cycles of PCR was performed. For samples with a Ct value between 30–31, 33 cycles of PCR was performed; the higher cycle number was needed to ensure there was enough fungal material for sequencing. Samples with a Ct value above 33 were considered to contain no fungal material (matching a no-template control). Samples were purified by SPRI, pooled, and sequenced on the Illumina MiSeq platform using the Illumina MiSeq Reagent v3 600-cycle (2 × 300 bp).

For the detection of Methanogen from the samples, we performed qPCR to detect the copy number of Methanobacteriales with specific primers (Fwd: 5’-AGGAATTGGCGGGGGAGCAC-3′, Rev.: 5′-TGGGTCTCGCTCGTTG-3′) targeting the 16S rRNA gene fragment between positions 915 and 1,100 (48, 49). First, we use a universal primer (Fwd: 5′-ACTCCTACGGGAGGCAGCAG-3′, Rev.: 5′-ATTACCGCGGCTGCTGG-3′) targeting the 16S rRNA gene between position 314 and 540 (50) to detect the number of total bacteria and archaea. Part of the E. coli 16S gene fragment (NR_112558) was synthesized and diluted to 10⁷, 10⁶, 10⁵, 10⁴, 10³, 10², and 10 copies as standards. The PCR program for total bacteria is initial at 95°C for 5 min, followed by 45 cycles of 10 s at 95°C, 10 s at 60°C, and 10 s at 72°C. Only the samples that had total bacterial copy numbers > 10⁵ were used to detect Methanogens. The standards for Methanogen were synthesized from part of the Methanobacterium espanolae 16S rRNA gene (NR_104983.1) and diluted to range from 10⁷ to 10 copies. The PCR program for Methanogen detection started with an initial step at 95°C for 10 min, followed by 40 cycles of 10 s at 95°C, 30 s at 57°C, and 20 s at 72°C. Both PCR was performed using a Quantstudio 3 system (Thermo Fisher, Waltham, MA, USA) with a Quantinova SYBR green PCR kit (Qiagen, Hilden, Germany) and 1uL of extracted DNA was used as the template. To avoid extrapolating beyond the standard curve, we set the detection limit at 10 copies for both PCR experiments.

Raw reads were demultiplexed and quality-filtered using the QIIME2 pipeline (v2022.2) (51). The quality-filtered reads were denoised and concatenated with DADA2 (52). Taxonomy was assigned to Amplicon sequence variants (ASVs) against the SILVA database v138.1 (53), using in QIIME 2 (q2-feature-classifier). The phylogenetic tree was generated by the FastTree algorithm (54). Negative control samples were also sequenced, and ASVs identified as contaminated by decontam (55) were removed from further analysis. The reads count was rarefied to 2,857 reads per sample (21 samples were excluded, Supplementary Table 1). The ASVs determined as contamination are listed in Supplementary Table 2.

Faith’s Phylogenetic diversity, observed ASV, and Shannon Index were calculated as alpha diversity metrics. The difference in alpha diversity between countries or continents was tested using the Kruskal-Wallis group test and with FDR p-value adjustment.

Jaccard, Bray Curtis, unweighted and weighted Unifrac distances were calculated to obtain pairwise beta-diversity, and dimensionality reduction on the distances was performed using Principal Coordinates Analysis (PCoA) method (ape v5.8). All alpha and beta diversity metrics are calculated on ASV level abundance. Permutational multivariate analysis of variance (PERMANOVA) was used to test the significance of different groups with 999 permutations (vegan v2.6-8).

Differentiated taxa were detected with ANCOMBC (2.6.0) (56) with Holm–Bonferroni correction, with a default prevalence cutoff of 10%. The shared taxa between countries were shown as an upset plot by UpSetR (1.4.0). The co-abundance network was performed on the genus level with taxa with abundance > 0.01 and prevalence > 10% using the SparCC method (57) (SpiecEasi v1.1.3), and the correlation cutoff was set to ± 0.3.

The 16S rRNA gene (V4 region) sequencing analysis of the 497 milk samples included in this work (Supplementary Table 3) yielded 10,767,189 high-quality filtered sequences in total, ranging from 18 to 76,462 reads per sample [mean = 20,786 reads per sample; median (IQR) = 18,069 (9,542–28,173) sequences per sample]. The samples were rarefied to 2,857 sequences per sample (Supplementary Figure 1).

First, we applied a linear model to analyze the effect of country as well as some potential cofounders like maternal age, maternal body weight, birth mode, and the postpartum days of milk collection on the alpha diversity. Only the country showed a significant effect on the milk alpha diversity, and the effect size is bigger than the other variables (Supplementary Table 4). Next, we only focused on the difference between countries. As assessed by using the Faith PD diversity index, the milk alpha diversity oscillated between Equatorial Guinea and Mexico (Figure 1A); the same range between Equatorial Guinea (lowest) and Mexico (highest) in Observed ASVs and Shannon index (Supplementary Figure 2). Overall, there was a significant difference in alpha diversity between different countries. Among all the pairwise comparisons, Ecuador and Mexico had significantly higher Faith PD diversity than Germany, Norway, Spain, the United States, and Equatorial Guinea, Equatorial Guinea had significantly lower alpha diversity than Peru, Senegal, South Sudan, Austria, Mexico, and Ecuador. Most other pairwise comparisons were not significant after adjustment for multiple comparisons. Similar trends are shown in Observed ASVs, but not in Pielou evenness or Shannon index (Supplementary Table 5). When aggregated into the Continent level, America showed significantly higher alpha diversity in Faith PD and richness (observed ASVs), but metrics accounting for evenness (Pielou and Shannon index) showed higher evenness in Africa (Supplementary Figure 3).

For beta diversity, we also included all potential cofounders in the first PERMANOVA analysis and found that only country showed a significant effect on beta diversity (Supplementary Table 6). Therefore, in the following analysis, we only included the country or continent. We detected a significant effect by country with the omega square effect size is 0.103 (meaning about 10.3% of the variation is from the country effect; p-value = 0.001). Mexico, Equatorial Guinea, and Norway showed the most segregated centroids in PCoA analysis based on Bray Curtis distance (Figure 1B). Some European and African countries were close to each other in the PCoA plot, including Spain, Germany, Austria, Kenya, Senegal, and Tanzania, indicating similar bacterial structures. The countries from different continents also differ in PC1, the countries from America all had negative PC1 values while most of the rest countries had positive PC1 values. Jaccard and Unifrac distance showed a similar effect, but the effect size of weighted Unifrac (0.120) is bigger than Jaccard (0.068) or unweighted Unifrac (0.064), suggesting the difference between countries was due to abundant taxa (Supplementary Figure 4). The effect of the continent on beta diversity was also significant but had a smaller effect size (omega square range from 0.022 to 0.043), and based on Bray Curtis and Jaccard, the distance of the center of samples from Europe and Africa was closer than the distance to America (Supplementary Figure 5).

In total 14,581 clean ASVs were generated from 16S sequencing, belonging to 41 phyla, 1,333 genera and 2,774 species. Most of them corresponded to seven major phyla: Bacillota, Pseudomonadota, Actinomycetota, Bacteroidota, Spirochaetota, Verrucomicrobiota, and Deinococcota, (Figure 1C). Among the 1,333 genera detected in this work, only 38 of them were found in samples from all the countries (Figure 2), while 18 were shared by all countries except for those from Tanzania and 7 were shared by all countries except for those from El Salvador. The genera with the top three highest frequency of detection were Staphylococcus (99.2% of samples; present in samples from all countries), Streptococcus (97.4% of samples; present in all countries), and Bifidobacterium (93.2% of samples; present in all countries). Overall, the top three most abundant genera were the same as the top three prevalent genera, Streptococcus (mean abundance 23%), Staphylococcus (mean abundance 16.8%), and Bifidobacterium (mean abundance 5.9%). Other highly abundant (> 1%) and universal genera include Pseudomonas, Corynebacterium, Acinetobacter, Rothia, Lactobacillus, Prevotella, Bacteroides, and unclassified genera from Muribaculaceae, Enterobacteriaceae, and Lachnospiraceae (Supplementary Table 7).

Samples from mainland USA are enriched in Pseudomonas, and depleted in Streptococcus, and Lactobacillus. Mexico and Ecuador are enriched in Bifidobacterium, Bacteroides, and Prevotella. Many countries from the European cohort are depleted in Bacillus and Alistipes (except Germany), enriched in Staphylococcus. Many countries from the African cohort are enriched in Alistipes and Bacillus and depleted in Bacteroides and Gemella. Kenya enriched in Acinetobacter, Alistipes, Lactobacillus, and Muribaculaceae (Figure 3 and Supplementary Figure 6).

Co-abundance network analysis of human milk bacteria only showed positive associations, particularly Bacteroides, Prevotella, and Bifidobacterium formed positive connections between each other. Other positive connections are between Staphylococcus and Corynebacterium; Streptococcus and Rothia; Pseudomonas and Acinetobacter (Figure 4A). That different bacterial genera occupied central positions, indicates country-specific structures and ecological roles within the microbial communities. For example, in Austria, Blautia and Rothia served as central nodes with strong associations to multiple genera. In South Sudan, Lactobacillus, Corynebacterium, and Kocuria were the most highly connected, while in Ecuador, Streptococcus, Veillonella, Staphylococcus, Blautia, Treponema, Bifidobacterium, and Prevotella showed the highest connectivity. These variations likely reflect underlying differences in genus abundance. Indeed, genera such as Bifidobacterium, Staphylococcus, Streptococcus, Lactobacillus, Corynebacterium, Rothia, Bacteroides, Prevotella, and Pseudomonas were identified as differentially abundant and consistently appeared as key hubs in their respective networks. Notably, the global network combining all countries’ data was the most robust and densely interconnected, highlighting overarching co-abundance patterns across diverse populations.

Despite these differences, we also observed some recurring correlation patterns across countries, similar to those in the global network. For instance, connections among Bacteroides, Prevotella, and Bifidobacterium were seen in Austria, Ecuador, South Sudan, and the United States. However, specific patterns varied: in Austria, Blautia and Muribaculaceae were also connected (Figure 4B); in Ecuador, Blautia, Succinivibrio, Treponema, and Bacteroides formed key links (Figure 4C); in South Sudan, Muribaculaceae showed strong associations (Figure 4F); and in the United States, Streptococcus emerged as a connected node (Figure 4H). Networks in Austria, Ecuador, Equatorial Guinea, and South Sudan exhibited greater complexity and connectivity compared to others.

Based on 16S sequencing, we were able to detect three genera from Methanobacteriales (Methanobrevibacter (0–2.6%), Methanosphaera (0–16.7%), and Methanobacterium (0–33.3%) Supplementary Table 8).

We also applied a standard curve qPCR method to detect Methanobacteriales. Based on our 16S rRNA sequencing results, which showed that the methanogen relative abundance ranged from 0.03 to 1% in samples where methanogen reads were detected. Given these proportions, if a sample contained fewer than 10⁵ total bacterial cells, the number of methanogen cells could be as low as 30 copies -close to the detection limit-. Therefore, we performed methanogen qPCR only on samples with a total bacterial cell count exceeding 10⁵. When examining a subset of the samples (268 samples from 7 countries) with qPCR, we found that only 21 (8%) had total bacteria/archaea 16S gene copy number >10⁵ per 1 μL DNA, evidencing the low bacterial density in these human milk samples from healthy women. Furthermore, this threshold of detection of 10⁵ gene copies found in our study likely overestimates the number of bacterial cells in the samples, since bacteria can have multiple 16S gene copies. Among the 21 samples, only two, one from South Sudan and the other from Equatorial Guinea had detectable methanogenic archaea (Supplementary Table 9).

The presence of fungi was evaluated in the whole collection of 518 samples (Supplementary Table 10). Only 6.6% of samples (n = 34) were positive based on ITS2 amplification and they belonged to the following cohorts: Equatorial Guinea (n = 12, 36.4%), Senegal (n = 1, 1.7%), South Sudan (n = 3, 5.7%), Tanzania (n = 3, 75%), Mexico (n = 2, 11.1%), mainland US (n = 7, 8.2%), Puerto Rico (n = 3, 60%), and Spain (n = 3, 3.9%). Among the fungi-positive samples, three fungal phyla (Ascomycota, Basidiomycota, and Mucoromycota) and 12 major genera were detected, dominated by Candida (Figure 4). Some of the fungal genera, such as Candida, and Malassezia, are frequent inhabitants of the human skin while others, such as Dichostereum or Blakeslea, are typically associated other to air, soil, and decaying plant matter, suggesting environmental or laboratory contamination of milk samples (Figure 5, Supplementary Table 11).

The interkingdom co-abundance network showed that Prevotella, Corynebacterium and Clavispora were the centers of the network (Figure 6A). Besides Clavispora, there were three fungal genera involved in the network, including Candida, Pichia, and Yarrowia. Most bacterial genera involved were the common genera found in the bacteria-only networks, including Bacteroides, Prevotella, Bifidobacterium, Staphylococcus, Corynebacterium, Streptococcus, Pseudomonas, Lactobacillus, and Rothia. Compared with the bacteria-only network with all samples (Figure 4A) or the same sample size as the interkingdom network (Figure 6B), the bacteria correlations remained similar, and the fungal genera acted as “bridges” connecting the isolated parts of the bacteria-only network.