Lifelines is a multi-disciplinary prospective population-based cohort study examining, in a unique three-generation design, the health and health-related behaviors of 167,729 people living in the North of the Netherlands. It employs a broad range of investigative procedures in assessing the biomedical, socio-demographic, behavioral, physical and psychological factors that contribute to the health and disease of the general population, with a special focus on multi-morbidity and complex genetics. For the purposes of this study, we used data available for a subpopulation of the Lifelines cohort, the Dutch Microbiome Project (DMP) (Gacesa et al., 2022). The DMP cohort is focused on establishing gut microbiome patterns in health and disease, and it includes additional questionnaires and biological samples collected between 2015 and 2016 alongside the Lifelines baseline data collection.

Within the DMP, we focused our analysis on 7,198 individuals for whom we have data on anthropometrics, antibiotic use, employment details, home urbanicity levels and occupational codes according to the International Standard Classification of Occupation (ISCO). In addition, all selected DMP participants also provided fecal samples, collected between 2015 and 2016.

Gut microbiome data was obtained from fecal samples through whole-genome shotgun sequencing. Samples were collected by the participants at home and frozen immediately. Upon arrival at the Lifelines biorepository, samples were stored at-80oC until DNA extraction. Information on the season of sample collection was also recorded and used in our analysis. The gut microbiome composition was characterized as previously described in Gacesa et al. (2022). In short, microbial DNA isolation was performed using the QIAamp Fast DNA Stool Mini Kit (Qiagen), according to the manufacturer’s instructions. Shotgun metagenomic sequencing was done at Novogene, using the Illumina HiSeq 2000 technology, and samples followed previously described preparation protocols (Gacesa et al., 2022). Metagenomes were profiled following the data analysis steps of the DMP (Gacesa et al., 2022), 1000IBD (Imhann et al., 2019) and Lifelines-DEEP (Tigchelaar et al., 2015) cohorts. In short, metagenomic profiles were filtered from human reads using the KneadData (v.0.10.0) tool. Taxonomic composition was obtained using MetaPhlAn3 (v.3.0.1) (Beghini et al., 2021). Lastly, samples with total read depth < 10 million and/or percentages of eukaryotic or viral abundance ≥25% were excluded.

Microbial alpha diversity was calculated at species-level using Shannon and Simpson indexes from R package “vegan.” To account for each participant’s unique gut microbiome composition, we focused our taxonomic association analyses on bacteria present in at least 20 individuals from the group of participants with exposure estimates to any class of pesticides, which resulted in selection of 168 bacterial species.

To derive estimates of cumulative and categorical occupational exposure to different classes of pesticides, we utilized data from work-related questionnaires, job codes and job-exposure matrices (JEMs). The job codes were retrieved based on the ISCO, using a Computer-Assisted Structured Coding Tool (CASCOT), following the process of Statistics Netherlands (CBS, 2008). These job codes provide a classification system for employment based on the tasks and duties undertaken during work. These codes were combined with the ALOHA+ JEM (Matheson et al., 2005; Skorge et al., 2009) to estimate the levels of exposure to different pesticide classes per job classification. The ALOHA+ estimates are based on the intensity and probability of exposure and yield a semi-quantitative ordinal scale that corresponds to three different levels: no exposure, low exposure and high exposure, numerically represented by a 0–1-2 scale. We obtained exposure estimates for herbicides, insecticides and fungicides and a more general class of pesticides we refer to as “general pesticides.”

For the majority of our population (66%), ISCO codes were retrieved based on employment information collected concomitantly with fecal collection. For participants for whom concurrent information was not available (34%), ISCO codes were retrieved from the baseline questionnaire data. Additional employment information, such as change in unemployment status and duration of employment/unemployment were used to check whether participants had any data discrepancies between both timepoints. If so, additional inspection of the data was made to confirm whether there were no inconsistencies between pesticide exposure estimates for each participant. Lastly, the participant’s ISCO codes, with coding from the 2008 data release, and the ALOHA+ JEM, with ISCO job coding from 1988, were linked via the publicly available correspondence table between distinct ISCO data releases (ISCO-08 and ISCO-88) (ISCO, 2008). Linkage between job codes with distinct pesticide exposure estimates, i.e., ISCO-88 codes that were mapped to multiple ISCO-08 codes that contained distinctive exposure estimates (or vice versa), were analyzed separately and final codes were selected based on similarity of job tasks (Koeman et al., 2013) (Supplementary Table S1).

We calculated cumulative exposure estimates by linking self-reported length of employment with the ALOHA+ estimates (Matheson et al., 2005; van der Mark et al., 2015). Employment questionnaires had categorical answers on length of employment (“Less than 1 year,” “Between 1 and 5 years,” “Between 6 and 10 years” and “More than 10 years”). To quantify these variables numerically, we assigned each category its corresponding mean value (0.5, 3, 7.5 and 15 years). Finally, we multiplied the length of employment by the weighted numeric values of the initial exposure assessment, with 4 indicating high exposure and 1 indicating low exposure, in line with methods used in previous reports (Matheson et al., 2005). In cases where participants were exposed to at least one class of pesticide, at any level, and there was no information regarding their length of employment (24% of the exposed group), we assigned these participants the mean numeric length of employment among the exposed group (11.6 years).

We used covariate-adjusted models to investigate the association between occupational exposure to pesticides and microbial alpha and beta diversity as well as to the abundance of microbial species. The models included adjustment for age, sex, body mass index (BMI), season of fecal sampling, mean Bristol Stool score, DNA concentration of the sample and antibiotic use in the past 3 months (Table 1). Comparisons between the results from the adjusted and unadjusted regressions demonstrate that adjusting our model for the covariates increased the power for discovery in the statistical modeling and led to only minor changes in effect estimates. Associations of pesticide classes to alpha diversity and species abundance were calculated using linear regression. For regressions with categorical exposure estimates, an overall assessment of the statistical significance for all categories was produced using ANOVA. Associations to beta diversity (represented as Aitchison dissimilarity) were calculated using PERMANOVA (adonis2 function from R package “vegan”). Prior to association analysis, the abundances of microbial species were transformed using centered log-ratio transformation (adding pseudo count equal to ½ of the smallest non-zero abundance in the dataset, similar to the procedure implemented in the R package “microbiome”) to properly normalize skewed distributions and deal with the compositional nature of the data. For the cumulative exposure estimates, we used inverse-rank sum transformations. All statistical test results were adjusted for multiple testing using the Benjamini-Hochberg procedure.

We analyzed data from 7,198 participants of the DMP cohort who reported being employed in at least one timepoint of the DMP data collection. The mean age of the population was 50.7 years, with ages ranging from 20 to 84 years. The majority of the cohort (59%) is female (Figure 1A). Among the participants in our study, 4.1% (n = 299) were assigned a level of exposure to at least one type of pesticide (Table 2).

The exposed participants had a similar age distribution to the total population, with a mean age of 50.9 years. A higher proportion of males (69.2%) were estimated to be occupationally exposed to pesticides compared to females in the same category (30.8%) (Fischer’s exact test, p < 2.2×10⁻¹⁶, Odds Ratio = 3.42) (Figure 1B). Occupational exposure to general pesticides, herbicides, insecticides or fungicides was reported in a total of 19 distinct job classifications. Of our exposed group, 37% (n = 111) were dairy and livestock producers, 15% (n = 46) were freight handlers, 13% (n = 39) did general garden-based labor, 9% (n = 28) were mixed crop and animal producers, 5% (n = 16) were field crop and vegetable growers and 3% (n = 11) were mobile farm and forestry plant operators, among other less common occupations.

Overall, low pesticide exposure estimates were more prevalent compared to high exposure estimates, across all pesticide classes (Table 2). Apart from herbicides, which displayed a smaller percentage of exposed individuals (only 1.5% exposed), all other pesticide types had similar exposure rates across the population (3.2% for fungicides, 3.9% for insecticides and 4.1% for general pesticides), and these exposures were often concurrent (3.0% exposed to all three types), resulting in high correlation between these pesticide classes (R_min = 0.86 and R_max = 0.98) (Figure 1C).

Additionally, data on length of employment was available for 79.8% of the total population. The majority of the participants (53.6%) reported working for the same employer for over 6 years, while only 26.2% indicated being employed for a shorter period. A similar trend was observed when looking at the employment profile across the exposed individuals, where 58.9% declared employment for more than 6 years. The average length of employment within the general cohort was 9.7 years, whereas this was 11.1 years for the exposed group. These observations suggest that the individuals who are occupationally exposed to pesticides had long and potentially repetitive exposure profiles.

Previous reports have demonstrated that other exposome factors, such as geographical location, socio-economic status and urban and farm-related environments are associated with changes in the human gut microbiome (Segata, 2015; Tasnim et al., 2017; Hildebrand et al., 2021; Gacesa et al., 2022; Valles-Colomer et al., 2023). Because occupational pesticide exposures occur most often in farm, agricultural and garden-related jobs, which are often characteristic of a less urban environment, we sought to clarify whether the occupational pesticide exposure estimates were correlated with the urban density of participant’s residences.

For the DMP population, urbanization levels were determined based on the home address at the time of fecal collection, as described in Gacesa et al. (2022), where the surrounding urban density was represented by a categorical scale from 1 (very urban, ≥2,500 addresses registered per km²) to 5 (very rural, ≤500 addresses per km²). To investigate the relationship between urbanicity levels and pesticide exposure, we calculated Spearman correlation coefficients between the density scale and all exposure estimates for all pesticide classes for both categorical levels and cumulative estimates (Figure 1C). While exposed occupations in our cohort were primarily linked to farm-related jobs, the urbanicity levels of their residential environment are not strongly correlated to any of the occupational exposure estimates (max. R_spearman = 0.13). This suggests that, in the northern region of the Netherlands, residential urbanicity indices do not have a substantial effect on occupational exposures.

To investigate the impact of occupational pesticide exposure on the gut microbiome, we compared the alpha diversity of the gut microbiome, measured using the Shannon and Simpson indexes, across the occupational exposure estimates to different pesticide classes and levels. Our results show that Shannon diversity was significantly decreased upon exposure to most pesticide classes. Fungicide exposure was found to be strongly associated with decreased taxonomic diversity (P_adj = 1.11×10⁻⁶), followed by general pesticides (P_adj = 1.12×10⁻⁶) and insecticides (P_adj = 2.15×10⁻⁵) (Figures 2A–C). In contrast, exposure to herbicides did not show any significant associations with changes in microbial alpha diversity (P_adj = 0.23) (Figure 2D). Similar results were obtained using the Simpson index (Supplementary Figure S1) The lack of statistical significance for herbicide exposure is most likely due to the much lower number of exposed participants compared to other classes of chemicals, as we still observe a trend of decreased alpha diversity.

We also observed that low exposure estimates were consistently more statistically significant across our associations compared to high exposure estimates. Here, the statistical power may again play a role in the significance of results, as only a small number of participants had high estimates of occupational exposure to pesticides. When looking at cumulative exposure estimates, we observed similar results to the categorical analyses, with exposure to fungicides, general pesticides and insecticides significantly associated with decreased microbial diversity but again no statistically significant association to herbicide exposure (Table 2).

To further analyze the gut microbiome beta diversity across individuals with distinct exposure profiles, we performed PERMANOVA analyses. This allowed us to explore the dissimilarities in the microbial community structure given the exposures to different classes of pesticides. In line with the alpha diversity results, exposures to general pesticides, fungicides and insecticides were significantly associated with a distinct gut microbiome composition (Table 3), while herbicide exposure was not significantly associated with overall microbial composition.

Our findings provide evidence that occupational pesticide exposure is associated with changes in the diversity and composition of the gut microbiome. We demonstrate that the impact of such exposures may follow different patterns according to distinct classes of pesticide exposure.

To inspect the effect of these exposures on specific gut microbial features, we explored associations between the different types of pesticide exposures and the relative abundance of the 168 bacterial species that passed our selection criteria. Our analysis included both the categorical levels of exposure and cumulative exposure estimates, and all statistical tests were corrected for age, sex, BMI, season of sampling, mean Bristol Stool score, DNA concentration of the sample and antibiotic use in the past 3 months (Table 1). Rather than relying solely on statistical associations for each level of exposure, we compared the difference between all categories of exposure for the same pesticide class using an additional ANOVA test.

In total, 39 bacterial species were significantly associated with exposure to at least one class of pesticide (Supplementary Table S2). Many of the significant associations between specific taxa and pesticide exposure estimates were consistent across the different classes of pesticides, with 26 bacterial species displaying significant associations with general pesticides, fungicides and insecticides.

We observed that exposure to all types of pesticides was significantly associated with the increased presence of Allisonella histaminiformans, Bacteroides coprophilus, Mitsuokella multacida and Parabacteroides sp. CAG 409 (Figures 3A,B). Conversely, the abundance of Barnesiella intestinihominis was decreased for exposures of all types (Figure 3C). Most of the impact of exposures was related to changes in the presence or absence of low-abundance bacteria like Acidaminococcus fermentans and Megasphaera elsdenii, but less commonly it was related to changes in the abundance of more abundant bacteria, such as decrease of Bacteroides dorei, Akkermansia muciniphila and Alistipes finegoldii (Figure 3D).

Of note, in the associations with specific taxa, we could see effects of herbicide exposure. However, from our analyses, only two bacterial species, Phascolarctobacterium faecium and Roseburia faecis, were significantly associated with herbicide exposure alone, and these signals were only marginally significant (P_aov-adj = 0.04 and 0.03, respectively). Additionally, these and other associations, such as with Allisonella histaminiformans and Barnesiella intestinihominis, only remained statistically significant after multiple-testing correction in tests with cumulative exposure to herbicides.