We conducted a cross-sectional analysis of a subset of female participants in the Growth and Obesity Cohort Study (GOCS). The original GOCS, which began in 2006, has been described previously (Corvalán et al., 2009). In brief, 1,196 children aged 2.5 to 4 years from low- and middle- income families in the South East region of the Santiago metropolitan area in Chile and enrolled in preschool at the National Board of Preschool Council Program (Junta Nacional de Jardines Infantiles) were enrolled in the study; approximately half (601) were girls. Participants in the study visit the Institute of Nutrition and Food Technology (INTA) Health Clinic at the Universidad de Chile in Santiago, Chile at least once per year for anthropometric assessments, pubertal (Tanner) evaluation, collection of biospecimens, and to complete 24-hour dietary recall interviews. A limited set of behavioral and demographic information was also collected via questionnaire. Breast composition was measured when the participants were two years post-menarche (2PM). The current study included a subset of 218 girls randomly selected to provide a fecal sample and who had a breast composition measurement at 2PM. The 2PM breast assessment and fecal sample collection occurred between 2018 and 2019. The study protocol and written consent forms were approved by the University of Chile Ethics Committee at INTA.

Breast composition was measured at the clinic visit corresponding to a timepoint of 2PM for each girl using the dual energy X-ray absorptiometry (DXA) breast scanning protocol developed by Shepherd et al. at the University of California, San Francisco (Shepherd et al., 2006). The Prodigy DXA system software (version 13.6, series 200674; GE Healthcare) was used to scan each breast for quantification of adipose fat and fibroglandular (FG) tissue. Stable calibration of the system was continuously performed using a quality control breast phantom. Absolute fibroglandular volume (FGV; cm3) and total breast volume for each breast were derived from a two-compartment model of adipose fat and fibroglandular tissue. The percentage of FGV tissue in the breast was derived by dividing the absolute FGV by total breast volume and multiplying by 100. Percent FGV and absolute FGV for the left and right breast were averaged to obtain two single measures of breast density: percent FGV (%FGV) and absolute FGV (aFGV). DXA is frequently used in studies of bone density in children; exposure to ionizing radiation from the DXA protocol is low (Crabtree et al., 2014). The DXA approach for breast composition assessment has high validity and precision among adolescent girls (Shepherd et al., 2008). Both breast density outcomes were categorized into terciles (‘T1’, ‘T2’, ‘T3’) based on the distribution of the sample for statistical analyses.

Fecal samples were collected after annual visits to the INTA health clinic occurring when the girls were between 13 and 15 years old. During the visit, girls were provided with materials (sealed plastic bag, stool catcher, a plastic, sterile container with a spoon to manipulate the sample and procedure gloves) and instructions for at-home fecal collection. Briefly, after the stool deposit, girls were asked to collect a part of the stool approximately the size of a walnut and place it in the plastic sterile container and label with the date and time of sample collection. Within 15 minutes of collection, the samples were sealed and then stored temporarily in the participant’s freezer, at which time study personnel were contacted to schedule the sample pick up. Samples were retrieved from the participant’s home no later than three days after initial collection, transported on ice to INTA, labeled with a de-identified key, and stored in a -80°C freezer prior to shipment. All stool samples were shipped on dry ice to the National Exposure Assessment Laboratory at Emory University.

Fecal processing of 279 samples (266 unique samples and 13 duplicates) was performed at the National Exposure Assessment Laboratory at Emory University, a Children’s Health Exposure Analysis Resource Laboratory. DNA was extracted from fecal samples using the Qiagen DNeasy PowerSoil Kit (Qiagen; 12888) with the TissueLyser at maximum speed (30 beats per second) for 10 minutes. Composition of the gut microbiome was determined through sequencing and amplification of the V3-V4 hypervariable region of the 16S rRNA gene (PCR primer pair: 341F 5’-GTGCCAGCMGCCGCGGTAA-3’; 805R 5’-GACTACHVGGGTWTCTAAT-3’). Libraries were made from 12.5 ng of DNA following a standard 16S Metagenomic Library Preparation Workflow from Illumina, Inc. Libraries were pooled in equal amounts based on fluorescence quantification, resulting in a 630 bp amplicon. Final library pools were quantitated via qPCR (Kapa Biosystems; KK4824). The pooled library was sequenced on an Illumina MiSeq using MiSeq v3 600 cycle chemistry (Illumina MS-102-3003) at a loading density of 6-8 pM with 20% PhiX, generating roughly 20M, 300 bp paired-end reads. In addition to the 279 experimental samples, traditional negative, no template control (NTC) negative, positive, and ZymoBIOMICS mock microbial community controls were included in the assays (Pollock et al., 2018). The Emory Integrated Genomics Core performed the assays.

After sequencing, demultiplexed raw amplicon sequences were processed using QIIME2 (Quantitative Insights Into Microbial Ecology) version 2019.4 (Bolyen et al., 2019). Denoising and dereplication, including chimera removal and trimming of reads based on Phred quality scores, were performed using the Divisive Amplicon Denoising Algorithm 2 (DADA2) module (Callahan et al., 2016a). The majority of Phred scores were greater than 25. To improve sensitivity of the algorithm, we trimmed the first 30 base pairs and truncated the reads at position 290, where the mean read quality distribution dropped for the overall dataset. Amplicon sequence variants (ASVs) were inferred using DADA2 to increase resolution and allow for intrinsic biologic meaning (Callahan et al., 2017). Taxonomy was assigned using a naive Bayes classifier on the SILVA database (SILVA 132 release) (Quast et al., 2013). We excluded samples from girls that did not have breast assessments at 2PM (n=48). The end product yielded a total of 6,270 unique ASVs from 218 samples. Standard preprocessing (filtering, subsetting, agglomeration) with the R package ‘phyloseq’ was used to exclude undefined (NA) or ambiguous (uncharacterized) phyla, phyla below a prevalence of 1, and ASVs below a prevalence threshold of 2% of total samples (McMurdie and Holmes, 2013; Callahan et al., 2016b). The final ASV feature table comprised 18,628,903 total reads (mean per sample = 87,872, range 28,118 to 278,065) and a total of 1,600 unique ASVs across 218 samples.

Demographic, anthropometric, and nutritional data were collected by trained dietitians during the annual study visit at the health clinic. Age- and sex- adjusted body mass index (BMI; kg/m²) Z-scores were calculated using the World Health Organization growth reference data and categorized into Norma/Underweight (Z-score ≤ 1), Overweight (1 < Z-score ≤ 2), and Obese (2 < Z-score). Body fat percentage was estimated using Tanita-BC-418 MA bioelectrical impedance measurements (Tanita-Corporation, Tokyo, Japan) and categorized into Healthy, Overweight, and Obese based on Tanita children’s age- and sex-specific body fat reference curves (McCarthy et al., 2006; Cediel et al., 2016). Age at menarche was determined via phone interviews completed by study dietitians every three months during puberty and dichotomized (≤ 12 years, > 12 years). Mothers of the participants were present at the clinic visits and were asked to complete short questionnaires with information on their highest level of education (secondary education or less, post-secondary education), birth mode of the participant (cesarean, vaginal), and months of exclusive breast feeding (<3 months, 3-6 months, >6 months). Dietary 24-hour recalls were collected longitudinally beginning in April 2014 by trained dietitians using the USDA multiple-pass method. Food and nutrient information were obtained using a harmonization process that mapped Chilean foods to the USDA Food and Nutrient Database for Dietary Studies. Daily intake of five major food groups (vegetables [g], fruit [g], red or processed meat [g], yogurt [g], and whole grains [g]) and total energy intake (kCal) were averaged over the data collection period and up to the date of the 2PM breast composition assessment to reduce random-measurement error and to obtain a more accurate assessment of long-term diet (Willett, 2012). We also collected data on ethnicity (Mapuche [native Chilean indigenous] and non-Mapuche, according to last name), average daily hours of television during the week (≤1 hour, 1-3 hours, >3 hours) as a proxy of physical activity, and antibiotic use in the 6 months prior to fecal sample (yes, no). Missing covariate data were imputed using last observation carried forward for body fat percentage and dietary variables. We assumed the remaining covariates with missing values, including maternal education, ethnicity, breast feeding, birth mode, antibiotic use, and television hours were missing completely at random. We used mean imputation (breast feeding) or median imputation (maternal education, ethnicity, birth mode, antibiotic use, television hours) in the R package ‘mice’ (van Buuren and Groothuis-Oudshoorn, 2011). The variable with the highest proportion of missing values was television hours (n=22, 10.1%); all other covariates were missing at a proportion less than 3%.

Relative abundance at the phylum level was plotted for all samples and across %FGV and aFGV terciles. We estimated alpha diversity as the observed richness (the number of species per sample) and the Shannon index (a measure of richness and evenness) on unfiltered data (Shannon, 1948; Chao, 1987). Rarefaction without replacement was used to standardize library sizes and account for uneven sampling depth prior to estimating alpha and beta diversity metrics and has been suggested to be suitable for microbiome data (McKnight et al., 2019; Estaki et al., 2020). Overall differences in observed richness and the Shannon index by %FGV and aFGV terciles were tested using the Kruskal-Wallis (KW) followed by post-hoc pairwise testing with the Wilcoxon Rank Sum in cases where data were compatible with evidence to reject the KW null hypothesis. To evaluate whether our results were sensitive to body composition, we performed analysis of covariance (ANCOVA) and adjusted for body fat percentage. Beta diversity was visualized using principal coordinate analysis (PCoA) plots with Bray-Curtis dissimilarity. Single and multivariate permutational analysis of variance (PERMANOVA) were used to test for overall differences in microbial composition between %FGV and aFGV terciles and produced marginal effects using the adonis2 function in the ‘vegan’ package (McArdle and Anderson, 2001; Oksanen et al., 2007). Homogeneity of variance was used to test whether differences in community structure were due to dissimilar dispersions. Supplemental PCoA visualization and PERMANOVA tests were done using weighted UniFrac, which differs from Bray-Curtis in that it incorporates phylogenetic relationships (Lozupone et al., 2007). Dietary intake of whole grains, fruit, vegetables, and yogurt are important predictors of gut microbial composition and may influence breast density. Therefore, we performed supplementary PERMANOVA analyses with added terms for four food-specific intake groups, including average whole grain intake (g/day), average fruit intake (g/day), average vegetable intake (g/day), and average yogurt intake (g/day). Microbial diversity metrics were estimated using the R packages ‘phyloseq’ and ‘vegan’ (Oksanen et al., 2007; McMurdie and Holmes, 2013).

ASVs from the filtered final feature table were agglomerated to the genus level prior to associating microbial community features with breast density outcomes. A total of 223 genera remained after agglomeration of the 1,600 ASVs. We used the R package ‘MaAsLin2’ (Microbiome Multivariable Associations with Linear Models), which relies on a modified generalized linear model for compositional data, to identify differentially abundant microbe genera across %FGV and aFGV terciles (Mallick et al., 2021). All models specified minimum abundance of 0.01% in 10% of samples, total sum scaling normalization, and arcsine square root-transformation (Mallick et al., 2021). The models were adjusted for a set of potential confounders including age, body fat percentage, antibiotic use, maternal education, total calories, hours of TV watching, ethnicity, birth mode, and breast feeding. Where appropriate, the Benjamini and Hochberg (BH) correction method was used to control the false discovery rate (FDR) and produced q-values (Benjamini and Hochberg, 1995). Supplementary MaAsLin2 analyses were also performed with adjustment for the four food-specific dietary intake groups.

Functional metabolic pathways of microbial communities were predicted from 16S rRNA marker sequencing data using PICRUSt2 implemented as a QIIME2 plugin (Douglas et al., 2020, 2). Default parameters were used in the pipeline, including SATe-enabled phylogenetic placement, hidden-state prediction, and a distance cut-off of 2. PCoA ordination with Bray-Curtis dissimilarity was used to visualize predicted abundance of metabolic pathways present in the MetaCyc Metabolic Pathway Database. Differences in MetaCyc pathway abundance between %FGV and aFGV terciles were testing using multivariable PERMANOVA. Multivariable MaAsLin2 was used to evaluate differential abundance of 458 predicted MetaCyc pathways with default parameters and BH correction for multiple comparisons. As in prior analyses, we additionally adjusted for the four food-specific dietary intake groups in supplementary analyses.

The association between breast density and gut microbial composition and function was assessed among 218 GOCS participants. The mean age of the participants was 14 years at the time of breast assessment and 12 years at reported menarche ( Table 1 ). The participants were primarily non-indigenous (Non-Mapuche, 83.0%) and 78.4% of the mothers had a secondary school educational attainment or less. The mean %FGV among the sample was 49.5% (SD=14.5); for tercile 1 (T1), mean %FGV was 34.0%; tercile 2 (T2), 48.4%; and tercile 3 (T3), 65.3%. The mean aFGV among the sample was 217.0 cm3; for T1, mean aFGV was 140.2 cm³; T2, 210.8 cm³; and T3, 292.7 cm³. More than half the cohort was considered overweight or obese. Overall, participants had a mean BMI Z-score of 1.0 and a mean body fat percentage of 32.6%. Additional study population characteristics are reported in Table 1 .

The relative composition of the gut microbiome for each %FGV and aFGV tercile is displayed in Figure 1 . Overall, Firmicutes was the most highly represented bacterial phylum and comprised 66% of the sample abundance. Microbes from the Bacteroidetes and Actinobacteria phyla represented 18.1% and 10.4% of the abundance, respectively. We found minor differences in relative abundance across %FGV terciles for Bacteroidetes (T1: 20.4%, T2: 15.4%, T3: 18.5%), Actinobacteria (T1: 8.9%, T2: 10.9%, T3: 11.5%) and Euryarchaeota (T1: 1.4%, T2: 0.9%, T3: 0.9%) ( Supplementary Figure 1 ). There were no differences in relative abundance at the phylum level across aFGV terciles ( Supplementary Figure 2 ).

We calculated alpha diversity for each sample using observed richness and the Shannon diversity index. Overall, the mean number of observed species present in the fecal samples was 214 and the mean Shannon index was 4.0 ( Figure 2 ). We observed differences in median observed species richness across %FGV terciles (T1: 226; T2: 200, T3: 223). There were no differences in the observed number of species across aFGV terciles nor in the Shannon index for both %FGV and aFGV terciles. Results from ANCOVA did not suggest a difference in alpha diversity across %FGV nor aFGV terciles when adjusting for body fat percentage (data not shown).

Beta diversity analysis for microbial ASV abundance was performed using Bray-Curtis dissimilarity and visualized by principal coordinate analysis (PCoA) ( Figure 3 ). PCoA axes 1 and 2 represented 9.2% and 5.8% of the total variance, respectively. Results from multivariable PERMANOVA analyses suggested that a small but statistically significant portion of the variability in microbial composition was associated with %FGV (R² = 0.012, p=0.02) ( Table 2 ). No significant contribution to variability was associated with aFGV (R² = 0.09, p=0.53). We also did not note any visual shifts in global predicted microbial pathway abundances ( Figure 4 ) nor any significant predictors of Bray-Curtis dissimilarity in PERMANOVA analyses ( Table 3 ). In supplementary analyses of beta diversity using weighted UniFrac, rather than Bray-Curtis, we found no significant association with %FGV and aFGV ( Supplementary Table 1 ). Supplementary PERMANOVA analyses with additional adjustment for dietary predictors did not substantially modify the variation in microbial nor pathway diversity ( Supplementary Tables 2 , 3 ).

In multivariable linear modeling with MaAsLin2, no genera were associated with %FGV nor with aFGV after FDR correction ( Supplementary Table 4 ). These null associations remained after further adjustment for additional dietary predictors, including whole grain, fruit, vegetable, and yogurt intake ( Supplementary Table 5 ). Additionally, there were no major differences in abundance of predicted MetaCyc pathways by %FGV nor aFGV when adjusting for potential confounders ( Supplementary Table 6 ) or important dietary predictors ( Supplementary Table 7 ).