A population-based birth cohort was established in the city of Sabadell (Catalonia, Spain) as part of the INMA Project (Guxens et al., 2012). Between July 2004 and July 2006, 657 pregnant women who visited the primary health center of Sabadell for an ultrasound in the 1st trimester were recruited. Inclusion criteria were: age at least 16 years, intention to give birth in the reference hospital, no problems in communication, singleton pregnancy, and no assisted conception. Informed consent was signed, and the study was approved by the ethics committee of the Institut Municipal d’Investigacio Medica (IMIM), Barcelona, Spain.

Among 657 mother–infant pairs initially enrolled, 622 participated in the follow-up and conducted at the time of delivery. The mother–child pairs were later followed in the third trimester of pregnancy, at delivery, and at ages of 6 months and 1, 4, and 7 years (participation rate in the last follow-up, 76%). As part of the 7-year follow-up, 154 out of 473 children (32.6%) provided stool samples and were included in the present analysis. A total of 152 children yielded good quality microbiome sequencing data, but only fecal microbiome data from 151 children with the required determinant information were further analyzed. A graphic representation of the workflow is presented in Figure 1.

We considered all the variables collected as part of the INMA study from pregnancy to childhood and filtered them based on previous associations in the literature, prevalence, repeat measures (for diet and pollutants), missingness, and collinearity, resulting in 69 variables to be included in the association analyses (Figure 2), pertaining to six categories: pollutants (14 exposure variables), diet, demographics, quality of living/social environment, clinical records (49 non-genetic determinants), and genetics (6 genetic determinants). Full description of the variables is presented in Supplementary File 1; Table S1.

Some of them were related to birth and early life conditions, such as cesarean section delivery, predominant breastfeeding, siblings at birth, having pets at 14 months, attending day care center at age of 2 years, and maternal BMI during pre-pregnancy and pollutants. Others were related to factors at age of 4 or 7 years (number of servings/day of animal protein, dairy products, fruits and vegetables, fiber foods or sweet products, physical activity, asthma, vaccination, viral infections, and pollutants— mercury and tobacco smoke exposure sections) or specifically at 7 years of age [grade of adherence to Mediterranean diet measured by the Mediterranean Diet Quality Index (Serra-Majem et al., 2004)—KIDMED index, children BMI, and siblings]. General conditions such as sex, ethnic origin, and maternal education level were also collected.

Data for each non-genetic determinant and pollutant exposures were available for 151 children. In the case of missing values, the value was imputed, as explained in the statistical analysis section. Original genetic data were available for a subset of 107 children. The very specific dietary data obtained from the mothers (fish intake at 32 weeks of pregnancy), available for 148 children, was used only for adjustment in one specific mercury analysis, avoiding multiple imputation.

All statistical analyses were performed using R software. Descriptive statistics of the original metadata were calculated as median (first quartile, third quartile) for continuous variables and as frequency (percentage) for categorical variables.

To avoid calculations of β-diversity with missing values, continuous variables were transformed and standardized to logarithm (log), interquartile range (IQR), or log of IQR, as appropriate. Then, missing values for exposures and determinants (excluding genetic, FUT2, and fish intake at 32 weeks of pregnancy data) were imputed using multiple imputations by chained equations to generate 20 imputed data sets (mice R package). Original data were used for genetic, FUT2, and fish intake at 32 weeks of pregnancy.

β-Diversity association analysis: In the first place, a multi-determinant analysis using the envfit function (vegan R package, 999 permutations) was employed to identify determinants, which were significantly correlated with the gut microbiome β-diversity. We included the following dissimilarity matrices at the genus level: wUnifrac, Bray–Curtis, and Aitchison (Gloor et al., 2017), which were calculated as the Euclidean distance of centered log ratio-transformed data after Bayesian Multiplicative replacement of count zeros. The envfit function performs multivariate analysis of variance (MANOVA) and linear correlations for categorical and continuous variables, respectively. It calculates the effect size (R²) and significance (p-value) of each variable, comparing the difference in the centroids of each group relative to the total variation. Mean and standard deviation of the R² values obtained along the 20 imputations for 14 tobacco and mercury exposure variables and 49 non-genetic determinants were calculated. The mean of p-values was also calculated. An association was considered statistically significant if a p-value < 0.05 was obtained in 17 out of the 20 imputations. In addition, for a sensitivity analysis, genetic data were included in an additional envfit analysis considering the Aitchison distance in a subset of 107 children in the same conditions.

In the second place, additional analyses were performed for each pollutant exposure variable without adjustment (crude model) for potentially important determinants. Exposures to tobacco were adjusted by maternal education, maternal BMI, and the envfit-selected variables. Exposure to mercury during pregnancy was adjusted by fish intake at 32 weeks of pregnancy, maternal education, and the envfit-selected variables. Exposure to mercury at 4 years was adjusted by fish intake at 4 years, maternal education, and the envfit-selected variables. A permutational multivariate analysis of variance (PERMANOVA) test by distance matrices with Benjamini & Hochberg (BH) correction for multiple testing using the adonis function (vegan R package, 999 permutations) was used. Mean and standard deviation of R² values obtained for the imputed datasets considering Aitchison distance were calculated. An association was considered significant when a BH-corrected q-value <0.05 was obtained for more than 17 out of the 20 imputed datasets. The mean of adjusted q-values was also calculated. The principal component analysis (PCA) plot with the clr-transformed features was obtained with the microViz R package for one imputation. Ellipses corresponding to 95% confidence intervals (CI) between categories were calculated.

Genus association study: At the genus level, taxa were normalized to relative abundance and regressed against individual tobacco and mercury exposures in a crude model, adjusting for the potentially important determinants selected for the PERMANOVA analysis, including the most important covariates selected by the envfit analysis, as described above (MaAsLin2 R package, by default conditions). Then, we selected individual exposures and their corresponding taxa from the raw MaAsLin2 output table and calculated the BH-corrected q-value. For the screening of significant taxa, a q-value <0.05 was considered statistically significant, as described previously (Shen et al., 2022).

After performing 16S rRNA sequencing and quality control, we acquired 34,444 to 285,423 sequencing reads per sample (mean ± SD: 83,143 ± 35,329 reads). Sequence data that support the findings of this study have been deposited in the EBI EMBL database (PRJEB69271). A total of 4,904 ASVs were detected, corresponding to 2 kingdoms, 13 phyla, 20 classes, 52 orders, 80 families, 228 genera, and 186 species. The core microbiota in Sabadell-INMA children at the phylum level, defined as phyla found in more than 85% of the samples, comprised 5 phyla with the following mean relative abundance: Bacteroidota 44.65%, Firmicutes 39.28%, Verrucomicrobiota 11.02%, Actinobacteria 2.74%, Proteobacteria 2.29%. The core microbiota at the genus level, defined as genera found in 95% of the samples, comprised 15 genera. The most abundant genera were Bacteroides (30.9%), Akkermansia (11.2%), Alistipes (11.1%), and Faecalibacterium (10.8%). We considered 207 ASVs for the association study corresponding to 2 kingdoms, 6 phyla, 9 classes, 20 orders, 30 families, 61 genera, and 58 species that were present in 20% of samples.

The description of the characteristics of the studied Sabadell-INMA birth cohort is shown in Supplementary File 2: Table S2. Exposures to pollutants (tobacco and mercury) were of particular interest. In the studied population, 46.6% of children were exposed to tobacco smoke in-utero, 30.3% during childhood based on the 4 year-old children data, and 17.9% along the full period (sustained exposure) based on questionnaires. The mean corrected cotinine levels were 224.1 μg/g for mothers and 17.5 μg/g for children. The median cotinine levels were 6.0 μg/g for mothers and 5.5 μg/g for children.

In this study, 65% of children had mercury levels in cord blood samples above 5.8 μg/L. At the age of 4 years, 38.4% of the children had hair mercury levels above 1 μg/g.

In the multi-determinant model, active smoking during pregnancy and SHS at 4 years based on detectable urinary cotinine levels (≥4 ng/mL) were associated with the β-diversity, explaining 2–4% of the variability (Aitchison distance, p-values: 0.0045 and 0.0322, respectively, Table 1). Having siblings at birth was also associated with the β-diversity (Bray–Curtis and Aitchison distances, p-values: 0.0214 and 0.0065, respectively). In the sensitivity analysis performed in a subset of 107 children considering also genetic data and Aitchison distance, active smoking during pregnancy and having siblings at birth remained associated with the β-diversity (p-values: 0.0016 and 0.0080, respectively). Genetic variables explained less than 5% of the β-diversity (Aitchison).

Only the statistically significant results are shown in Table 1. All results for the analysis of childhood and prenatal determinants for the three distances and the sensitivity analysis are presented in Supplementary Table S3. We found no significant associations of β-diversity with mercury or other tobacco exposure variables.

The corresponding mean of R² for 20 imputations obtained with the envfit analysis for all determinants and pollutant exposures are represented in Figure 3. All R² and p-values are presented in Supplementary File 3: Table S3 and Supplementary File 4: Table S4.

On the other hand, β-diversity measured by the Aitchison distance was not significantly associated with the 12 tobacco exposure variables or the 2 mercury exposure variables included in this study in the single pollutant model (PERMANOVA). An unadjusted model was considered and a model adjusting by the covariates above-listed (see Material and methods section), considering “having siblings at birth” as the covariate selected by the envfit analysis. The results obtained for each of the 20 imputations and both models are shown in Supplementary File 5: Table S5. Plots of PCA with the clr-transformed features are shown by the exposure to active smoking of their mothers at any time during pregnancy and having siblings at birth in Figures 4 and 5, respectively.

Tobacco smoking during pregnancy was associated with the microbiota at 7 years of age at the genus level. The significant results are shown in Table 2. Active smoking at any time during pregnancy was associated with a decrease in the relative abundance of Akkermansia both in the unadjusted and the adjusted models (coef = −3.76, q-value 0.004, coef = −3.75, q-value 0.005, respectively). Sustained maternal smoking from first to third trimester of pregnancy was associated with an increased relative abundance of Dorea both in the unadjusted and adjusted models (coef = 2.01, q-value 0.006, coef = 2.04, q-value 0.01, respectively). Figure 6 shows the changes in the abundance of Akkermansia and Dorea associated with tobacco smoking during pregnancy considering the adjusted model. In addition, maternal smoking categorized as non-sustained, low dose, or high dose was associated with an increased relative abundance of Dorea (coef = 1.79, q-value = 0.03), although only in the unadjusted model. Mercury or tobacco exposure from pregnancy to childhood measured by other variables considered in this study were not associated with significant changes in the abundance of genus in the children’s gut microbiome at 7 years of age. The best 10 results (genus) with the smallest q-values obtained by the MaAsLin2 analysis are shown in Supplementary material 6: Supplementary Table 6.

In the studied population, a high level exposure to tobacco smoke in utero and during childhood was found.

Median quantification of mercury from cord blood was equivalent to 4.5 μg/L in maternal blood, which was above the maximum mean value reported through the HELIX cohort (0.7–3.9 μg/L) (Aurrekoetxea et al., 2013). In this study, Mercury levels in cord blood samples were above 5.8 μg/L for a high proportion (65%) of children as found in communities with high fish intake such as Japan, Hong Kong, Taiwan, and Polynesia (geometric median 9.8, 8.8, 9.2, and 10.5 μg/L, respectively) (Ruggieri et al., 2017). Higher blood cord mercury levels have only been found in Canada or Greenland (geometric median 18.5 and 25.3 μg/L, respectively). Hair mercury levels above 1 μg/g were detected for a high proportion of children (almost 40%). These high values could be related to the place of residence of the studied population, where high fish and other aquatic food intake was the rule. Exposure to both tobacco smoke and upper quartile level of mercury was found in 12.1% of children at birth and 3.45% of 4-year-old children. For example, in comparison, only 8% of the children included in the DEMOCOPHES cohort implemented in 17 European countries exceeded this health-based value (Esteban et al., 2015).

The gut microbiota in healthy children of 7 years old showed a predominance of Bacteroides, Akkermansia, Alistipes, and Faecalibacterium, as described previously (Xiong et al., 2022). Based on the results of a multideterminant model, active smoking during pregnancy was the most associated determinant with the children’s gut microbiome β-diversity (Aitchison distance). Among all the other possible determinants included in the study, having siblings at birth was significantly associated with the β-diversity. Nevertheless, no significant associations between smoking variables with the β-diversity were observed in the univariate adjusted model.

We tried to determine the significant factors driving differences in the individual gut microbiota at the genus level. Maternal smoking during pregnancy was related to a decreased relative abundance of beneficial genus (Akkermansia) and an increased relative abundance of potentially pathogenic genus (Dorea), suggesting the importance of the environmental determinants during pregnancy, as a potential window of vulnerability. We did not find any other persistent effects at the genus level for tobacco exposure during childhood.

This study assessed cord blood and hair mercury levels as markers for mercury absorption and, consequently, mercury exposure. These levels elevate in correlation with the mother’s and child’s fish consumption, respectively, as well as the mercury content in the consumed fish species. Therefore, fish intake was regarded as a confounding factor in the analysis. It has been well established that intestinal microbiota can influence methylmercury intestinal absorption and the induction of intestinal dysbiosis (Pinto et al., 2020). On the other hand, methylmercury could have adverse effects on gut bacteria and accelerated accumulation of mercury in organs due to disruption of gut microbiota. Nevertheless, a sustained influence of early life mercury exposure could not be demonstrated in this study at the genus level. As it was suggested that mercury could decrease the growth of Lactobacillus in different proportions depending on the species (Seki et al., 2021), an additional study at the species level using a shotgun strategy could be relevant. In addition, as changes in the overall composition of the microbiome have been described for mercury exposure during early gestation, but not at late gestation (Rothenberg et al., 2016), it would be interesting to study other indicators of early exposure to mercury. Finally, the exposure to mercury at the most vulnerable stage of growth, during brain development, could lead to neurotoxicity in the unborn infant by a mechanism independent of permanent changes in microbiota. Higher childhood or perinatal blood mercury has been related to higher relative abundance of potential pathogens (e.g., Flavonifractor plautii), beneficial species (e.g., Bifidobacterium longum, Faecalibacterium prausnitzii), and both potentially pathogenic and beneficial species (e.g., Bacteriodes vulgatus, Eubacterium rectale), but these results were not found in our study (Shen et al., 2022).

Cigarette smoking has been consistently linked to alterations in the gut microbiome among adults (Cicchinelli et al., 2023). The impact of maternal smoking during gestation on the prenatal and perinatal gut environment is also suspected. These alterations may be influenced by shifts in the mother’s own microbiome or variations in the in utero chemical environment. Consequently, such changes could potentially reshape the trajectory of the microbiome development of offspring, rendering them more susceptible to inflammation. For example, changes in the abundances of several genera of the child microbiome of the finger, nose, mouth, and ear canal have been reported after exposure to third-hand smoke. These genera are also affected by active smoking and second-hand smoke (e.g., Corynebacterium, Staphylococcus, Streptococcus) (Kelley et al., 2021). Similarly, the abundance of Serratia, Moraxella, Haemophilus, and Staphylococcus aureus in tracheal aspirate changed after tobacco exposure (Leroue et al., 2023). The longitudinal effects of environmental tobacco smoke exposure on gut microbiota of young children have been previously reported. Tobacco was correlated with Megasphaera abundance. Breastfeeding was correlated with Lactobacillus abundance and negatively correlated with Clostridium_sensu_stricto_1 abundance (Xie et al., 2021).

Previous microbiome studies have also shown the inverse relationship between the abundance of Akkermansia muciniphila, a mucus-layer-degrading bacterium, and diseases such as inflammatory bowel disease, obesity, and diabetes. Its potential as immunomodulatory probiotic for autoimmune and chronic inflammatory diseases has been explored in experimental models. This species seems to slow down the development and progression of diabetes, obesity, and IBD in mice (Rodrigues et al., 2022). On the other side, Phascolarctobacterium has been previously associated with colorectal cancer (Bucher-Johannessen et al., 2023). This suggests a long term effect of tobacco early-life exposure on gut microbiome of children.

It was previously reported that at 1, 3, and 6 months of age, Staphylococcus, Ruminococcus, Akkermansia, Lachnospiraceae, and Bacteroides were more dominant in the gut microbiota of infants and neonates exposed to tobacco smoke (McLean et al., 2019). Further studies would be necessary to explore the effect of smoking on specific species or strains of these genera by using shotgun metagenomic analysis, allowing a more detailed view of the microbiome and its function.

Several strengths of our study can be highlighted. First, highly sensitive biomonitoring methods were used for tobacco and mercury assessment, combined with a longitudinal design (repeated measurements). In addition, this population has particularly high exposure to tobacco and mercury, which makes it more likely to detect potential effects, even for a relatively small size in the omics field. The INMA children are a well-characterized, general population, where the width of parameters was measured, which is relevant for gut microbial diversity description. It was particularly important to be able to compare the effect of pollutants of interest, a novelty in the field of microbiome determinants, with established factors, such as diet and clinical factors.

The main limitations of this study are that samples were stored at 4°C for more than 1.5 h, and that the information of antibiotic intake in the previous month was not recorded. Considering that strong clustering by individual is expected in fecal samples because inter-individual variability could be greater than the variability associated with storage time (Holzhausen et al., 2021) and based on the results of alpha diversity of this dataset (data not shown), we considered that our samples were adequate for an exploratory analysis. Besides, unknown exposures that could be influencing the microbial composition as other chemicals, endocrine disruptors, heavy metals, air pollution, pesticides, food additives, air pollution in the place where children were raised were not included.

The causal link between the determinants and gut microbiota composition could not be inferred from this association study, as associations might be confounded by other factors. Exposure data collected at the age of 4 years might have changed at microbiota assessment age.