ABSTRACT
In a longitudinal agricultural community cohort sampling of 65 adult farmworkers and 52 adult nonfarmworkers, we investigated agricultural pesticide exposure-associated changes in the oral buccal microbiota. We found a seasonally persistent association between the detected blood concentration of the insecticide azinphos-methyl and the taxonomic composition of the buccal swab oral microbiome. Blood and buccal samples were collected concurrently from individual subjects in two seasons, spring/summer 2005 and winter 2006. Mass spectrometry quantified blood concentrations of the organophosphate insecticide azinphos-methyl. Buccal oral microbiome samples were 16S rRNA gene DNA sequenced, assigned to the bacterial taxonomy, and analyzed after “centered-log-ratio” transformation to handle the compositional nature of the proportional abundances of bacteria per sample. Nonparametric analysis of the transformed microbiome data for individuals with and without azinphos-methyl blood detection showed significant perturbations in seven common bacterial taxa (>0.5% of sample mean read depth), including significant reductions in members of the common oral bacterial genus Streptococcus. Diversity in centered-log-ratio composition between individuals' microbiomes was also investigated using principal-component analysis (PCA) to reveal two primary PCA clusters of microbiome types. The spring/summer “exposed” microbiome cluster with significantly less bacterial diversity was enriched for farmworkers and contained 27 of the 30 individuals who also had azinphos-methyl agricultural pesticide exposure detected in the blood.
IMPORTANCE In this study, we show in human subjects that organophosphate pesticide exposure is associated with large-scale significant alterations of the oral buccal microbiota composition, with extinctions of whole taxa suggested in some individuals. The persistence of this association from the spring/summer to the winter also suggests that long-lasting effects on the commensal microbiota have occurred. The important health-related outcomes of these agricultural community individuals' pesticide-associated microbiome perturbations are not understood at this time. Future investigations should index medical and dental records for common and chronic diseases that may be interactively caused by this association between pesticide exposure and microbiome alteration.
INTRODUCTION
Oral microbiome taxonomic diversity and community composition are variable both within and between individuals, but as a whole microbiomes are relatively stable in an individual adult developed host (1). Xenobiotic-microbiome interactions have been targeted as important to drug metabolism and biotransformation (2) which also makes the microbiome an important interaction term in the study of toxicology. Understanding of the microbiome, environmental exposure, and potential impacts on human health is currently growing: It has been elucidated that the microbiome can affect cognition (3–5), metabolic syndrome and obesity (6–8), and social development (9, 10). Many common environmental exposures (11–18) can change the microbiome composition; these can be as simple as dietary differences (13, 14) and can be induced by probiotics (16) and antibiotics (11). It has also been shown that community composition changes in the microbiome are persistent when they occur as the result of extinction events and are rescued only when both nutrients (diet as environment) and microbiome transplantation (reseeding or exposure) allow recovery (19). Additionally, changes in the microbiome can have particularly compelling effects on tertiary phenotypes, which can affect the well-being of individuals. For example a Danish cohort study of prenatal exposure to antibacterial agents found an increased risk of being obese and overweight later in life in individuals with prenatal exposure to systemic antibacterials (20). Another study, in Canada, showed that the alteration of microbial compositions affects the risk of childhood asthma (21). These reports point at the microbiome as a common determinant of other phenotypes with pleiotropic effects on human health. Understanding the environmental effects of exposure on the microbiome is thus particularly important for public health.
Organophosphate pesticides (OPs), including azinphos-methyl, are a class of toxicants that have had favor in the agricultural industry due to their relatively low persistence in the environment (22, 23). Acute OP toxicity in higher organisms operates primarily via a common mechanism of covalently binding a serine residue in the active site of the acetylcholinesterase enzyme (AChE); inactivating this site results in a systemic inhibition of enzyme activity (24). AChE cleaves the neurotransmitter acetylcholine, which is released in synapses when they fire. By inhibiting the enzyme cleavage of acetylcholine, synapses are not allowed to reset, which may result in death at high doses in adults and in neurodevelopmental abnormality in children at lower doses. Residential studies of exposure to the OP chlorpyrifos have shown reductions in working memory, reductions in intelligence, and brain anomalies in the babies of exposed mothers (25–29), as well as in laboratory animal experiments (30). Studies of acute exposure to these compounds have resulted in restrictions on the use of some OPs in both occupational and residential settings. As a result, the EPA has phased out any use of azinphos-methyl as of 30 September 2013 (http://www.epa.gov/opp00001/reregistration/Azinphos-methyl/phaseout_fs.htm).
The ∼200-household Yakima Valley community agricultural cohort studied here provides seasonal longitudinal measures of exposure to OPs using blood samples collected in the spring/summer of 2005 and winter 2006 for adults and children (31). Sequenced buccal samples from 117 adult individuals were used to determine proportional abundances of oral microbiome bacteria detected in buccal swabs. Blood concentrations of azinphos-methyl have been determined by mass spectrometry in this family-based community pesticide exposure study (31). Samples were collected in spring/summer 2005 to measure pesticide exposure during the spring-to-summer fruit-thinning season, when the workers had a high likelihood of OP exposure. The winter 2006 samples assessed exposure during the off-season, when the likelihood of exposure was low. Seventy individuals had buccal samples 16S rRNA sequenced from both longitudinal collection seasons. Bacterial detection was performed by DNA sequencing of the V5 and V6 regions of the 16S rRNA gene from 206 adult buccal swab samples. Counts of reads mapping to operational taxonomic units (OTUs) (32) with >97% sequence identity were converted to sample proportional abundances. The microbiota OTUs were assigned bacterial taxonomies by the Ribosomal Database Project (RDP) Classifier (33). Proportional taxonomic data fall into the category of compositional data analysis with simplex geometry, where the correlation structure of reductions in abundance of one bacterial taxon is autocorrelated with compensatory proportional increase in another taxon due to the remaining read counts being occupied by other compositional members (34). The “centered log ratio” of each sample's taxonomic proportional abundances applied here has been proposed as an appropriate transformation to correct for this compositional nature of the data (34–36). To avoid the complications due to limits of detection of rare member of the microbiota when read depths vary between samples, we have chosen to analyze associations in only common taxonomic genus members of the community composition (>0.5% mean proportion of sample read depth) as defined by the spring/summer 2005 samples. We adopted a slightly less stringent definition for common (>0.1% sample mean OTU proportion) when conducting exploratory principal-component analysis (PCA) at the OTU taxonomy resolution, in contrast to the genus-level taxonomy used in hypothesis testing of exposed groups. To model read depths and the probability that we would detect a taxonomic proportion, we used the negative binomial probability, Pr[missed taxon] = e−pn, where p is the proportion of reads identified to that taxon and n is the total number of reads for that sample (37, 38). This allowed us to establish criteria from the samples' taxon mean proportions (e.g., genera > 0.005, OTUs > 0.001) for which taxa to investigate at the read depth of sequencing sampled. Using only deeply sequenced samples (minimum threshold, 2,500 reads; median, 38,454; mean, 55,321), we fix a high probability at the median read depth, Pr[missed 0.001] < 1.99 × 10−17, that we will observe the identified common taxon in the Dirichlet process if the taxon exists in the sample above the proportion threshold. This common taxon analysis with deeply sequenced samples mitigates the false zero inflation of nondetects, which can occur when analyzing various depths of reads sampled and rarely detected taxa. We have also adopted this common genus detection limit threshold to limit the number of hypotheses tested to only common microbiota; to further reduce the possibility of bias, nonparametric Wilcoxon's rank sum tests were used to test for differences in microbial abundances. Using this approach, we identified robust false-discovery-adjusted associations between the identified microbiota abundance and azinphos-methyl exposure in this Hispanic Yakima Valley (Washington State) community agricultural cohort. We additionally extended this result to an exploratory PCA of the centered-log-ratio common OTU sample proportions present in buccal microbiomes of pesticide-exposed individuals and other members of the Yakima Valley Hispanic community.
RESULTS
Azinphos-methyl blood exposure measures.The blood concentrations of azinphos-methyl were measured by gas chromatography-linked high-resolution mass spectrometry (GC-HR-MS) (39). Of the 117 individuals presented here, azinphos-methyl was detected in 36 farmworkers and one nonfarmworker for whom we also had buccal microbiome data passing quality control. The detected concentration range in individuals was 0.021 to 6.192 ng azinphos-methyl/g blood serum. The limit of detection was 0.04 ng/g of plasma. One sample had a value less than the limit of detection that was nonzero (0.021); this sample was treated as a nondetect. The remainder of the assayed nondetects were zeros. In File S1 in the supplemental material we have formatted the range of the numbers such that values from the limit of quantification (0.04) to 1 ng/g (the median) are coded as 0.1 and values greater than 1 are coded as 1.1 ng/g plasma. This was to help ensure another level of privacy for the human subjects in the published form and does not affect the results, as we use the data categorically as detect or nondetect in this analysis. Only farmworkers had detected levels of azinphos-methyl in the spring/summer 2005 blood serum samples. All sampled individuals from the spring/summer collection had an azinphos-methyl blood assay result. Blood samples collected in the winter 2006 nonspray season included one nonfarmworker sample with a detected concentration of azinphos-methyl. This individual was a nondetect in the matched spring/summer 2005 blood sample. Thirteen of 95 individuals with winter microbiome data did not have blood sample azinphos-methyl assay results from the winter blood samples. They did have azinphos-methyl blood results from the spring/summer collection, and this value was used to categorize these individuals as to their blood detection category in the spring/summer season microbiome analysis. Since we did not have a winter blood azinphos-methyl assay value for these individuals, they were removed from the winter microbiome analyses. File S1 presents the azinphos-methyl blood detection for those individuals for whom buccal 16S rRNA sequencing also was performed.
Buccal bacterial taxonomic assignments.We observed DNA bands at ∼300 bp in size by agarose gel electrophoresis of the PCR amplicon of the 16S rRNA V5 and V6 regions, as expected. Sequencing of the PCR amplicons was performed with the Ion Torrent DNA sequencer. Reads were quality control discarded if the primers were not detected during trimming, if the read had more than three expected errors, or if chimeras were detected (32). The remaining mean read length was 247 bp (221 to 279 bp), putting the effective error rate at less than 0.012 error per base after quality control. The median sequencing depth of the 196 adult samples with more than 2,500 reads passing quality control was 38,454 total reads (mean, 55,321), making these samples deeply sequenced. Operational taxonomic unit (OTU) clustering by USEARCH (32) with 97% read DNA sequence identity resulted in 2,520 OTUs with measured differences between taxonomy members in excess of 7.41 nucleotide bases [(1 − 0.97) × 247 bp], with fewer than three of these being sequencing errors. These OTUs were classified to 286 bacterial reference taxonomies with the Ribosomal Database Project (RDP) Classifier java program (33). The RDP classifier reports to the lowest taxonomic level it can, down to the genus level. Among the common taxa we analyzed, suborders Micrococcineae and Actinomycineae and the family Coriobacterineae are identified as the lowest taxonomy level classified for these phylogenies. We refer to these taxa collectively as “genera” with the genera in this text. Sample bacterial DNA sequences were mapped onto the 2,520 de novo OTU reference taxonomies with USEARCH (32). Reads that did not match an OTU were assigned to the category OTU_NA and treated like the other OTUs. File S1 provides a summary of read sequence counts for each sample's OTUs as classified.
Nonparametric buccal bacterial analysis by azinphos-methyl blood detection group.The microbiome data from buccal samples in spring/summer 2005 (n = 101) and winter 2006 (n = 82) were categorized to factor groups by detection of azinphos-methyl in blood. Twenty-nine spring/summer individuals and 26 winter individuals for whom sequencing was performed had concurrent winter or previous spring/summer blood azinphos-methyl detection by mass spectrometry as evidence of pesticide exposure. Each buccal sample's genus proportion was transformed by the centered log ratio, and the common genera (n = 22 genera, by spring/summer sample mean genus proportion of >0.005) were tested by the Wilcoxon rank sum test (34, 35) between the sample groups in which azinphos-methyl was detected and not detected. Table 1 presents the false-discovery rate (FDR)-adjusted (40) significance, the estimated nonparametric location difference, and the unadjusted 95% confidence interval (CI) for each genus's centered log ratio abundance between these two sample groups for both seasons.
Wilcoxon's rank sum test of azinphos-methyl blood detection groups for compositional perturbation, showing significance results for common (>0.5%) generaa
Seven taxa (Streptococcus, Micrococcineae, Gemella, Haemophilus, Halomonas, Actinomycineae, and Granulicatella) were significantly reduced in the azinphos-methyl group at an FDR of <0.1 in the spring/summer 2005 data (Table 1, bold). Two of these (Streptococcus and Halomonas) maintained suggestive negative association significance (P values of ∼0.01, FDR of ∼0.14) in the following winter (Table 1, underlined). We observed winter values with the same direction but lower magnitude than the spring/summer collection values for Wilcoxon's location difference between groups showing reduced abundance in those with detected azinphos-methyl. This is not surprising, as of the 26 individuals in the winter azinphos-methyl detection group, only one (the nonfarmworker) had concurrent winter exposure measured in blood. This suggests a relaxation but persistence of the exposure effect into the winter. Mannheimia, Herminiimonas, Actinobacillus, Neisseria, and Prevotella appear also to be suggestively perturbed in the spring/summer azinphos-methyl group. Actinobacillus and Prevotella were suggestively (0.1 < FDR < 0.2) perturbed taxa that were decreased in the azinphos-methyl detection group, whereas the other three were increased in relative abundance. In Fig. 1 (spring/summer 2005) and Fig. 2 (winter 2006), each tested common taxon's mean proportional abundance in buccal samples by blood azinphos-methyl detection group is presented.
Spring/summer 2005 common bacterial genera (>0.5% of sample mean) from buccal swab sequencing, organized by spring/summer azinphos-methyl blood mass spectrometry detection into exposure groups with greater than and less than the limit of detection (0.04 ng/g of plasma). Genera are ranked by the mean proportion of reads. Box plots in the left panel are in log scale (with zero genus proportion samples set to 1e−6); black bars, boxes, and whisker bars show the median, the central quartiles, and the extreme values, respectively. The right panel and the central text report the group mean proportion to show the relative abundances of various genera in relation to each other. AZM, group with azinphos-methyl detected; ND, group with azinphos-methyl not detected. Listed among the genera, Micrococcineae and Actinomycineae are suborders and Coriobacterineae is a family (see the text for more information).
Winter 2006 common bacterial genera (>0.5% sample mean) from buccal swab sequencing, organized by azinphos-methyl blood mass spectrometry detection into exposure groups with greater than and less than the limit of detection (0.04 ng/g of plasma) in blood samples from either season. Common genera are defined and ranked by the mean proportion of reads from the spring/summer collection (Fig. 1) to facilitate comparisons between taxa in Fig. 1 and 2. Box plots in the left panel are in log scale (with zero genus proportion samples set to 1e−6); black bars, boxes, and whisker bars show the median, the central quartiles, and the extreme values, respectively. The right panel and the central text report the group mean proportion to show the relative abundances of various genera in relation to each other. AZM, group with azinphos-methyl detected; ND, group with azinphos-methyl not detected. Listed among the genera, Micrococcineae and Actinomycineae are suborders and Coriobacterineae is a family (see the text for more information).
Samples sequenced at very high depth have zero proportions of some genera. This can be seen in the whisker bars for many genera which extend to the 10−6 extreme value. The lowest common population genus proportions measured that were nonzero in individual samples were ∼10−5, two orders of magnitude below the mean proportion threshold and one order of magnitude above the zero order of the log scale. Two model assumptions were used in the application of the nonparametric Wilcoxon rank sum to hypothesis testing of the taxonomic proportion centered log ratios by azinphos-methyl detection groups: (i) we deeply sequenced samples with many reads, and (ii) only common population sample mean taxonomic proportion members were tested. These assumptions will mitigate the zero-inflation statistical effects that can occur in count-based proportional data at various depths sampled. As an example, Streptococcus was absent in three individuals (H161A, H153A, and H125A) with azinphos-methyl detected in the spring/summer collection with sequenced samples at the depths (2,695, 3,884, and 5,788 reads; Pr[missed at 0.005], 1.4 × 10−6, 3.7 × 10−9, and 2.7 × 10−13, respectively), whereas Streptococcus was detected in all individuals from the group with azinphos-methyl not detected. Two samples (H160A and H149A) from the winter collection with azinphos-methyl detected also had no Streptococcus detected in their buccal microbiota data. The depth of sequencing for these samples (84,158 and 30,611 reads, respectively) gave very low probabilities (∼10−183 and 10−67, respectively) at the mean proportion threshold (0.005) that we missed Streptococcus in these individuals, suggesting that the true proportion is indeed very low in many of these exposed individuals. This analysis suggests that taxon proportions that are very low in some individuals and proportionally different between the exposure groups are likely due to a true difference in the composition of the microbiota community structure. Streptococcus is a taxonomically large common oral microbiota member. Some members of Streptococcus have been shown to be beneficial to oral health (i.e., S. salivarius, S. oligofermentans and S. mitis) while others have been implicated in oral disease (i.e., S. mutans, S. pneumoniae, S. sanguinis, S. sobrinus, S. gordonii, S. agalactiae, S. tigurinus, and S. anginosus) (41–53). This split in the beneficence profile of Streptococcus causes a diversity of disease and positive health potentials for the observed effect when the compositionally abundant member Streptococcus is perturbed. Similarly, Haemophilus (54–62) and other genera identified here common to the oral context (63) can have a diversity of disease and health profiles within any taxonomic level and can be contextually associated with different coinfections, plasmids, species, and strains (64). This makes drawing specific conclusions about the perturbed community member contributions to health difficult without specific information on health outcome in this agricultural context and additional investigations.
PCA exploration of oral buccal microbiota proportional abundances.The buccal samples in the spring/summer 2005 and winter 2006 microbiome data revealed two primary clusters of oral microbiome types in each season. A two-dimensional PCA of the samples' centered log ratio of OTUs detected at >0.001 mean proportional abundance is plotted in Fig. 3. The left panels depict the spring/summer scores for PC1 and PC2. The right panels depict the winter OTUs in the spring/summer PC loading space as supplementary out-of-sample PC1 and PC2 projections. In the plots for both seasons, the point colors represent azinphos-methyl detection in blood as red and orange versus blue nondetects in the top panels, farmworker occupational status as red versus nonfarmworkers as blue in the middle panels, and the Gaussian model-based PC1 cluster calls as red for exposed (E) and blue for unexposed (U) clusters.
Oral buccal microbiome differences between individuals in bacteria projected by PC1 and PC2 of the centered-log-ratio-transformed proportions of common OTUs (>0.001 mean proportion of spring/summer sample reads). The top panels show mass spectrometry results for azinphos-methyl blood detection as colored points (red, >1 ng/g; orange, <1 ng/g; blue, none detected) overlaid at the PC coordinates of the buccal microbiome composition. The color coding in all panels represents azinphos-methyl detection in either season so that the spring/summer azinphos-methyl detection may be visualized with the winter microbiome PCA. The azinphos-methyl limit of detection is 0.04 ng/g of plasma. The middle panels show the occupational status of individuals. The bottom panels depict the PC1 Gaussian model-based cluster calls of the data projected in the PC coordinates. Two groups are observed in both seasons; the left cluster group is enriched for farmworkers and azinphos-methyl detection in blood samples. We interpret the left cluster as “exposed” based on the coclustering of individuals with blood detection of azinphos-methyl and the enrichment of farmworkers. The OTU percentage of variance explained (∼18%) by PC1 and PC2 is reported in the axis labels. The right winter panels are supplementary out of sample projections through the spring/summer PC1 and PC2 loadings.
We observed two clusters evident in both seasons. The less aggregated left cluster is visibly enriched for those individuals with detected azinphos-methyl in the blood (Fig. 3, top panels) and farmworker status (middle panels). A similar clustering pattern is also seen in the following winter 2006 buccal microbiomes (right panels), suggesting a persistence of the microbiome state induced by previous exposure. Azinphos-methyl was detected in only one nonfarmworker in the winter. If azinphos-methyl was detected in either season, the color in the plots for both seasons reflects this detection to track migration of pesticide-exposed subjects' microbiome types while those same individuals are azinphos-methyl nondetects in the winter. t tests of the PC1 scores for spring/summer and winter showed that the microbiome composition PC1 score was significant for azinphos-methyl detection in the left direction of PC1. The spring/summer mean azinphos-methyl group PC1 score was −2.8, versus 1.1 in the nondetect samples. In the winter projection, the azinphos-methyl group PC1 mean score was −1.6, while the nondetected samples had a PC1 mean score of 0.73. The spring/summer P value of 4 × 10−7 and winter P value of 1.7 × 10−4 show a persistent significant association between the buccal microbiome composition PC1 score and azinphos-methyl detection. The data points for the samples with lower blood concentration (<1 ng azinphos-methyl/g blood, the median) are colored orange to show that there is no simple evidence of a dose response within the exposed individuals. t test P values of >0.49 result from the comparison of PC1 values for high (red) and low (orange) azinphos-methyl concentrations in spring/summer and winter. The middle two panels of Fig. 3 show the farmworker-versus-nonfarmworker distribution in the OTU microbiota PCA. In the spring/summer, the t test of the PC1 scores for farmworker occupational status also showed that farmworkers had significantly lower mean microbiome PC1 scores (−1.01) than nonfarmworkers (1.16) (Welch's t test, P value of ∼0.01). This farmworker occupational PC1 microbiome exposure was not persistent in the winter (P value of ∼0.45). The bottom two panels present the mClust cluster call assignments using the spring/summer 2005 PC1 as input to Gaussian model-based cluster assignments (65). To train the model, the spring/summer 2005 PC1 data were used to determine a decision value of 0.72 based on the midpoint between the maximum PC1 value in the spring/summer left cluster and the minimum PC1 value in the right cluster. This midpoint (0.72) between clusters was then used to classify the winter 2006 buccal samples in the respective clusters based on the projection PC1 of winter OTUs through the spring/summer PC1 loadings. This allowed us to compare PC1 cluster memberships between seasons in the space defined by the spring/summer data. The left spring/summer and winter clusters are assumed as “exposed” based on the coclustering of samples in which azinphos-methyl was detected in blood. Categorical analysis of occupation status and cluster membership showed that the spring/summer 2005 buccal microbiomes' left exposed PC1 cluster contains 57 individuals (38 farmworkers and 19 nonfarmworkers) and the right, more tightly grouped cluster contained 44 individuals (16 farmworkers and 28 nonfarmworkers). The difference of cluster membership was significant by Fisher's exact test for occupational farmworkers status (P value of ∼0.003; odds ratio, ∼3.5, 95% CI, 1.4 to 8.7). When the 29 farmworkers and one nonfarmworker with both buccal microbiome data and detected azinphos-methyl concentrations in blood were compared with spring/summer buccal microbiome cluster membership, the left cluster, enriched with farmworkers, contained nearly all (27/30) of the individuals with spring/summer microbiome data with azinphos-methyl detected (including the one nonfarmworker with detection in the winter) and 30 individuals in whom none was detected. The right cluster contained almost exclusively individuals with nondetects of azinphos-methyl (41/44), except for three farmworkers with azinphos-methyl detected in blood. Those with detected levels of azinphos-methyl were significantly enriched in one of the two clusters from the spring/summer 2005 buccal collection (P value of ∼6.5 × 10−6; odds ratio, ∼12.0; 95% CI, 3.2 to 67). We refer to the left cluster enriched with farmworkers and containing the majority of individuals with azinphos-methyl detection as the “exposed” cluster and to the right cluster enriched with nonfarmworkers as the “unexposed” cluster. In winter 2006, 48 individuals (28 farmworkers and 20 nonfarmworkers) cluster in the left exposed PC cluster and 34 individuals (14 farmworkers and 20 nonfarmworkers) in the right unexposed cluster (Fisher's exact test P value of ∼0.18; odds ratio, ∼1.98; 95% CI, 0.74 to 5.4). Comparison of azinphos-methyl blood detection in any season to the winter 2006 microbiome PC clusters showed that azinphos-methyl detection in blood was still significantly enriched in the exposed cluster (Fisher's exact test P value of ∼0.008; odds ratio, ∼4.4; 95% CI, 1.4 to 17.2), with 21 individuals having azinphos-methyl detected and 27 with none detected, compared to the unexposed cluster containing five individuals with azinphos-methyl detected and 29 with none detected.
Many individuals with azinphos-methyl not detected in the blood (spring/summer, 30/57; winter, 27/48) cluster with individuals in whom it was detected. This would be expected, as azinphos-methyl is metabolized and excreted with a half-life clearance of ∼30 h (66), making the temporal window in which one must capture an exposure important. We analyzed the persistence of cluster membership using the 70 individuals with both spring/summer 2005 and winter 2006 microbiome sequence data. The majority (94%, 66/70) of individuals sampled remained in their original cluster the following winter. The one nonfarmworker with azinphos-methyl detected in the winter of 2006 remained as an exposed cluster member in both seasons. This overall pattern suggests a persistence to the exposed microbiome state as detected by model-based PC1 clustering state. The two individuals who migrated from the exposed to the unexposed microbiome cluster were a farmworker and a nonfarmworker. The farmworker who changed from the exposed to the unexposed cluster had azinphos-methyl detected in blood in the spring/summer of 2005 but none detected in the winter. Among just the exposed spring/summer cluster individuals, the composition effect of the spring/summer azinphos-methyl-exposed cluster as measured by PC cluster membership at the later winter collection date had 95% (39/41) remaining in the exposed cluster for both seasons. This cluster membership concordance in combination with the extremely low bacterial taxon abundances measured for some exposed individuals suggests that these microbiomes have undergone extinctions of the affected microbiota with minor regrowth. The consequence of bacterial extinction with little regrowth has been shown in mice when gut recolonization is not facilitated by microbiome seeding (19). Two nonfarmworkers with no azinphos-methyl detected switched from unexposed to exposed microbiome PC clusters, suggesting exposure events since the last spring/summer's pesticide applications which we did not ascertain. We speculate that more of the individuals from the unexposed cluster would have been found to migrate to the exposed cluster if sample collections had continued to include a following spring/summer 2006 agricultural cycle, when pesticide spraying begins again. The nonfarmworkers' microbiomes that cluster with individuals with detected azinphos-methyl blood concentrations suggest that other, unobserved exposures have occurred throughout the spring/summer and fall and into winter, as observed in the one nonfarmworker with detected azinphos-methyl in the winter 2006 collection. We suspect that stochastic differences in exposure in the home and community environment, as well as differential geographic location, azinphos-methyl dose, and clearance duration, may explain the heterogeneous clustering of nonfarmworkers with farmworkers in the respective clusters. File S1 contains all reported sample microbiome PC coordinates, cluster membership calls, occupation categories, and azinphos-methyl concentrations.
Abundances of commonly detected OTUs by PC loadings and clusters.The top-ranked 21 PC loadings of bacterial OTUs for spring/summer are presented in Fig. 4.
The ranked PC1 and PC2 vector loadings of the top 21 OTUs in spring/summer 2005 are colored based on the modeled cluster (red, exposed; blue, unexposed) and with number-labeled vector groups (1 to 4) based on which in graph direction the OTUs are the highest in abundance. File S1 in the supplemental material reports the vector loading values.
In the spring/summer collection there are two primary vector groups, labeled 1 and 3, among exposed individuals; representative taxa present in these individuals are depicted in Fig. 4. Vector group 2 includes taxa that are increased in individuals for whom we did not detect azinphos-methyl in the blood. The high-loading-rank OTUs in the spring/summer 2005 data are representatives of several genera also detected in the nonparametric analysis of the common genus composition in the azinphos-methyl blood detection group. These include the genera Streptococcus, Gemella, and Haemophilus and the suborder Micrococcineae. The Streptococcus OTUs are still present, but at lower abundance, in the exposed cluster and appear more in the top portion of the plot for vector group 1 (spring/summer). Streptococcus is the most increased in samples in the direction of vector group 2 (spring/summer) to the unexposed cluster. Neisseria is a common driver in the left exposed cluster direction of PC1. While not associated by the nonparametric analyses, it is possible that Neisseria is less susceptible to the exposure and, as the most common microbiota detected, is among the majority of the reads left when other taxa are not present in the compositional relative abundance. This leaves Neisseria as even more common in the exposed individuals. Neisseria is a large genus with many diverse members, many of which do not normally cause disease (67). The PC cluster discerning taxa is split as to which Fusobacterium OTUs are indicator taxa for the respective exposure clusters, reflective of the Wilcoxon test, where it is not significant (Table 1; Fig. 1 and 2). OTUs in loading vector groups 1 and 3 are pointed in opposite directions of PC2 and not at the central mass in the left direction of PC1. This heterogeneity in PC2 OTU genera may be explained in part by complex relationships in the composition. The perturbed taxa may have lower abundance but are still present in some exposed individual subjects, as seen for Streptococcus in the nonparametric test. This PC2 heterogeneity may be interpreted as a less tractable community ecology and homeostatic dysbiotic state that would occur when a normal apex community member such as Streptococcus is depleted. This community perturbation would have effects on many members of the community that are reliant on amino acids from syntrophic metabolic niches (68) generated by Streptococcus. The winter exposed cluster projected from these spring/summer OTU loadings (Fig. 3) also shows this less aggregated “dysbiotic” pattern, where individuals are more heterogeneous in PC2 if they have a low score in PC1.
Alpha diversity.Using both the azinphos-methyl detection group subject identities and the cluster identity calls from PC analysis, we investigated the alpha species richness diversity of the samples using the method of Chao and Shen (69). The Chao diversity metric estimates the total number of organisms in a sample (count of OTUs) given an adjustment based on the number of rare members observed. Figure 5 presents box plots of the seasons by azinphos-methyl detection group and exposure cluster definitions for Chao alpha diversity values.
Chao diversity metrics for exposure categories. The left panels depict the diversity of each season's azinphos-methyl blood detection (AZM) (red) and nondetection (ND) (blue) groups. The right panels depict each season's PCA cluster group diversity, labeled as “exposed” (red) and “unexposed” (blue). The black bar is the median, with the boxes as the central quartiles. Whisker bars are the extreme values. The Chao alpha diversity metric is an estimate of the total number of organisms given the number of low-frequency OTUs detected. Significance was determined by Welch's two-sample t test.
There is no significant difference in the OTU richness or diversity in the individuals with azinphos-methyl detected in blood in either season. However, within both the spring/summer and winter sample collections, the exposed PC cluster had significantly reduced counts of organisms compared to the unexposed cluster as seen in the Chao diversity (Welch's t test, P < 0.009).
DISCUSSION
We identify both a direct association and a principle-component clustering pattern of individuals indicating that the compositional diversity profile of the oral buccal microbiome is associated with environmental exposure consistent with an agricultural signal. This agricultural signal was measured in members of the Yakima Valley Hispanic community in Washington State by blood concentrations of the organophosphate insecticide azinphos-methyl and self-reported status as apple and pear orchard farmworkers. These exposed individuals have reduced abundances of several common oral indicator genera in this molecular epidemiology analysis of the buccal microbiota profile of agricultural pome fruit orchard workers and community members. The spring/summer 2005 sample collection timing targeted occupational farmworker and community exposures at times when seasonal agricultural applications of azinphos-methyl insecticide were known to occur. The microbiome PC cluster group membership of nonfarmworkers and farmworkers exposed to azinphos-methyl suggests that many individuals in the community may have unobserved exposures to pesticides which we did not capture in blood. We observed one such exposure in a nonfarmworker with blood azinphos-methyl detection from the winter 2006 collection season, when exposure is not likely based on the seasonal nature of pesticide applications. This underlines the need for data sets with fine temporal resolutions to discover modes of exposure for community pesticide exposures. Previously published household dust data from this same cohort suggest that azinphos-methyl exposures may occur via take-home pathways in cars and houses, where dry, stale air and low UV light conditions can allow environmental persistence (70).
Only three nonfarmworker households had self-reported resident smokers, making smoking unable to explain the scale of exposure difference associated here with the microbiome composition. Soil bacteria have been known to degrade organophosphates since the early 1970s (71–75) and have been well studied with the intent of developing a bioengineered solution to environmental remediation and industrial uses (76). The earliest report (2013) that OPs such as azinphos-methyl may cause changes in the composition of a host microbiome were conducted with chlorpyrifos in an intestinal in vitro simulation reactor model and in vivo in the rat (77). Azinphos-methyl, chlorpyrifos, and other OPs are closely related compounds with the same covalent mechanism of primary toxicity, cholinesterase inactivation. The same authors showed that chronic chlorpyrifos exposure changes the gut permeability and permits translocation of bacteria across the rat intestine (78), and they quickly followed this in vivo rat model study with a perinatal exposure study in which they also showed that dysbiotic microbiota and exposure affected the intestinal development of the rat pups (79). Two other very recent articles show that other OPs cause changes in the gut microbiome. One of the OPs, trichlorfon, induced intestinal dysbiosis and reductions in many taxa at low doses in the Japanese quail (Coturnix japonica) (80). In another study, low-level dosing of mice with the organophosphate diazinon induced gut microbiota and metabolomic changes in both sexes (81). These orthogonal studies provide evidence to support our finding that human exposure to OPs is associated with the composition of the human buccal microbiome. These studies also found significant reductions in the abundances of Lactobacillus, the name-defining genus member of the order Lactobacillales, of which Streptococcus is also a common phylogeny member. We also found another common genus member of the order Lactobacillales, Granulicatella, to be significantly decreased in abundance (Table 1, spring/summer) (−0.89; 95% CI difference, −1.64 to 0; P value of ∼0.026; FDR, ∼0.08) and an indicator OTU in cluster analysis (Fig. 4). This suggests that the order Lactobacillales may in general be susceptible to organophosphate pesticide exposure. The in vivo rat studies reported recently (77–79) also showed increased intestinal permeability, reductions in tight-junction gene expression, and increased translocation of bacteria to distal organs in chlorpyrifos-dosed animals. The in vitro digester microbiome gut model did not contain a membrane equivalent of the intestine or an immune system and showed differences in the bacterial composition of human feces used to inoculate the culture system upon organophosphate exposure to chlorpyrifos (77). These orthogonal findings suggest that at least two independent mechanisms of microbiome perturbation exist, host membrane permeability and direct microbial action. While our findings are limited to the oral-buccal body site, these other authors have found that similar OP-induced perturbations can affect the gut microbiome and gut-derived neurotransmitters that are systemically circulated. Gut metabolomic analysis of mice dosed with the OP diazinon indicates that there are changes in gut-derived neurotransmitters and amino acid metabolism products. These organophosphate microbiome-mediated neurotransmitters can affect both host behavioral phenotypes and the niches of other gut bacteria present (81). Future functional validations of antibiotic and other etiological properties by in vitro and in vivo methods could distinguish between the membrane permeability, immune and inflammatory responses, and potential for direct antimicrobial properties, as these other studies suggest. The systemic bacterial ecology and interaction with the human immune system and tight-junction tissue sites may contribute to the holistic compositional effect in the buccal context. In any of these cases, agricultural pesticide exposure is associated here with significant community composition changes in these Hispanic individuals' buccal microbiomes. The community ecology of organisms significantly dependent on Streptococcus will likely be affected, as it is the second most common oral member detected with significant reductions in abundance.
We specifically timed the sampling of the population in the spring/summer and winter seasons to capture the azinphos-methyl exposure well in the respective seasons to well inform a molecular epidemiology exposure effect study such as this. The farmworkers were actively working in the recently azinphos-methyl-sprayed orchards, thinning the fruit and pruning, in the April to July 2005 spring/summer sample collection time, as shown by direct blood detection as the measure of exposure. The health impacts of these agricultural exposure-associated bacterial community structures and composition in the oral context and at other body sites (i.e., the gut, nasal passages, etc.) have not been investigated in this cohort to date. Hispanics have a disproportionate disease burden for many common diseases related to vasculature (i.e., diabetes and cardiovascular disease) (82–90) that have also been shown to have microbiome-related disease variables in many other experimental contexts (91–104). Dental health, bacterial infection, and peripheral vascular disease have many independent associations (105), and some of the etiology could be affected by systemic loosening of tight junctions as induced by organophosphates (78). In many other contexts, environmental pollutants have been shown to affect the microbiome (106), and here we show in humans the association between the composition of the microbiome and agricultural pesticide exposure.
Conclusions.Given the large microbiome differences at the oral buccal site, the Hispanic cohort studied in this paper should continue to be followed for repeated measures of buccal samples, with additional sampling, dental history, and investigation of the pleiotropy to other health outcomes that may have resulted from these agricultural exposures. These future investigations should also include molecular functional studies on individual bacterial taxa and in community ecology system biology contexts for their susceptibility to organophosphates. Additionally, more evaluation of the oral microbiome composition as a simple-to-collect molecular biomarker of past exposures as a whole for use in epidemiological studies would be able to greatly inform public health and practice.
MATERIALS AND METHODS
Agricultural cohort.The individuals in the cohort are ∼50% self-identified Hispanic farmworkers in apple and pear orchards and ∼50% Hispanic nonfarmworker community members. The cohort was established in 2005 and consisted of 200 adults and 200 children from the same households (31); 117 of these individuals were studied for this buccal microbiome analysis report. The Center for Community Health promotion of the Fred Hutchinson Cancer Research Center collected the samples from the human subjects in the Yakima Valley agricultural region of Washington State. Samples are housed by the University of Washington's Center for Child Environmental Health Risks Research.
Seasonal sampling design.Human oral buccal swabs and intravenous blood samples were collected in spring/summer 2005 and winter 2006. The spring/summer 2005 sample collection was timed for spring through early summer (April to July 2005), when farmworkers thin the fruit trees and were potentially exposed to pesticides recently applied to control insects. The winter 2006 sample collection was timed for the off-season (December 2005 to March 2006), when pesticides were not in use. Each participating adult provided written informed consent. The Fred Hutchinson Cancer Research Center's Institutional Review Board reviewed and approved the collection procedures for this study (file IR 5946).
Measurements of azinphos-methyl blood concentrations.Isotope dilution gas chromatography-linked high-resolution mass spectrometry (GC-HR-MS) (39) was used to determine blood serum concentrations of azinphos-methyl at the Pesticide Laboratory in the National Center for Environmental Health at the Centers for Disease Control.
Collection of oral buccal mucosa samples.The Catch-All sample collection swab (Epicentre, Madison, WI) was used for oral buccal mucosa sample collection. Participants were asked to not smoke, eat, or drink (except water) for 30 min before sampling. The oral cavity was rinsed twice with water before sampling. The swab was moved circularly inside the mouth on each cheek 20 times and then stored in a 15-ml polypropylene tube with the tip down in 1 ml of RNAlater (Qiagen, Venlo, Netherlands). Samples initially were stored in the field office at −10°C and then were moved to a −80°C freezer at the University of Washington laboratory.
DNA preparation.Buccal swab samples were thawed and diluted in 9 ml of phosphate-buffered saline (PBS). The swabs were agitated in a circular motion continuously while the technician viewed the PBS for oral buccal debris turbidity, and then the swab was discarded. The PBS suspensions were then centrifuged at 400 × g at room temperature for 10 min, and the supernatants aspirated and discarded. The manufacturer's protocol for the Maxwell 16 buccal swab LEV DNA purification kit (Promega, Fitchburg, WI) for sample DNA extraction and elution was performed. Quantification of the sample extract's double-stranded DNA (dsDNA) concentrations was performed with the Quant-iT PicoGreen dsDNA assay kit (Thermo Fisher Scientific, Waltham, MA) according to the manufacturer's instructions.
16S rRNA PCR and bead DNA cleanup.Eurofins (Eurofins Scientific, Luxembourg) primers, 10 ng of sample DNA extract template, and EmeraldAmp GT PCR master mix (Clontech Laboratories, Inc., Mountain View, CA) were combined in PCR tubes with nuclease-free water to primer concentrations of 0.9 μM in a volume of 50 μl. Molecular-grade water negative-control blanks were also thermocycled with the same reagents. The 16S ribosomal subunit RNA gene primers were designed following the gene placement described by Cai et al. (107). Primer sequences targeted the 16S rRNA gene variable region 5 (V5) (ATTAGATACCCNGGTAG) and variable region 6 (V6) (CGACAGCCATGCANCACCT). Primer design was bidirectional with four oligonucleotides in total, the two above-mentioned sequences and their reverse complements. Primer synthesis included the Ion sequencing adaptors and DNA barcodes. All primers with barcodes and sequencing adaptors are listed in File S1 in the supplemental material. PCR mixtures were thermocycled using an MJ Research PTC 200 Peltier thermal cycler (MJ Research, St. Bruno, Quebec, Canada) with the temperature cycle protocol (step 1, 95°C for 10 min; step 2, 95°C for 30 s; step 3, 50°C for 30 s; step 4, 72°C for 1 min 15 s; step 5, go to step 2, 1×; step 6, 95°C for 30 s; step 7, 68°C for 30 s; step 8, 72°C for 10 min; step 9, go to step 6, 34×; step 10, 72°C for 10 min; step 11, 4°C) and is also in File S1. PCR amplicons were purified with the Agencourt AMPure XP PCR purification system (Beckman Coulter, Brea, CA).
Buccal DNA PCR amplicon pooling.Sample PCR amplicon fragment length and concentration for sequencing were quantified with the Agilent high-sensitivity DNA kit (Agilent Technologies Inc., Santa Clara, CA). Amplicons with peaks greater than the high-sensitivity DNA reference markers were diluted and reanalyzed. Amplicon templates were diluted to a concentration of ∼26 pM in low-Tris-EDTA (low-TE) buffer and pooled for sequencing.
Template preparation, enrichment, and sequencing.DNA template preparation of the PCR amplicons was performed using the Ion PGM Template OT2 400 kit (Life Technologies, Foster City, CA) on the Ion OneTouch 2 system (Life Technologies), and sequencing was on the Ion Torrent PGM System (Life Technologies) with the Ion PGM Sequencing 400 kit (Life Technologies) and Ion 318 Chip v2 (Life Technologies). The default Ion Torrent Browser template program for 16S metagenomics target sequencing was used.
Read sequence preprocessing.Sequences for each barcoded sample were provided in individual FASTQ files by the Ion Torrent native software suite. Quality control included removing reads without both 5′ and 3′ primer sequences, primer trimming, and removing reads with greater than three expected errors based on the PHRED base quality score (32). This reduced the error per read on an Ion Torrent system to less than 0.012 (3/247 bp mean read length).
De novo 16S rRNA OTU reference generation.A single FASTQ file was concatenated from all samples and quality controlled as stated above, resulting in a FASTA file. Exact duplicate-read dereplication was performed via the USEARCH software (32) with the -derep_fullength and -sizeout flags. The FASTA file was organized in descending order of number of observations, with singleton sequences being removed using the USEARCH -minsize 2 -sortbysize functions. Sequence reads were then clustered to operational taxonomic unit (OTUs) of 97% sequence identity, and the selected cluster consensus sequence representative for each OTU cluster was output to an OTU FASTA file with the USEARCH -cluster_otus function. Chimera read removal using the USEARCH -uchime_ref -strand plus flags and the -db gold.fa reference file removed PCR artifacts from the OTU FASTA file. This file was labeled with the USEARCH fasta_number.py script. The alignment target file for later use in individual sample read sequence alignments to OTUs was then made with the USEARCH -makeudb function. OTU reference sequences assignments to bacterial taxonomies down to the genus level were done with the default settings of the RDP Classifier v2.6 java program (33). At a read sequence alignment of 0.97 in OTU generation and with a mean read length of 247 bp, there will be in excess of ∼7.41 bases different between OTUs, with fewer than 3 of the base differences between OTUs being sequencing errors inflating the count of rare members of the OTU taxonomy. This margin between OTU sequence identity definition and known sequencing error profile gives a balanced taxonomic resolution in the face of known errors that all sequencers have and can be modeled efficiently with a PHRED score (108, 109).
Sample 16S rRNA read sequence alignments.Read sequences were quality controlled as stated above. Two-step chimera detection with the USEARCH -uchime_denovo and the USEARCH -uchime_ref -strand plus flags with the -db gold.fa reference file removed PCR artifacts. Each sample's reads were aligned by USEARCH with the command-line flags -usearch_global -strand both -id 0.97 to the de novo OTU references to result in a count of reads matching each OTU.
Data analysis.Custom Perl scripts were used to organize taxonomy, sample OTU count, and sample categorical descriptor tables from the 16S rRNA gene USEARCH alignment files and RDP Classifier into formats imported to R for association and categorical analysis. The R package phyloseq (110) was used with these tables for species richness Chao alpha diversity measures. Analyses in phyloseq were conducted at the OTU classification level. When conducing total diversity measures that depend on the numbers of singletons, doubletons, and so on in the tail of the distribution, no proportional frequency cutoffs were employed. Whenever analyses required statistical comparisons between orthogonal detection groups (i.e., azinphos-methyl blood detection), OTUs were pooled by genera identified with the RDP classifier; genera identified at >0.5% mean proportional abundance in the spring/summer 2005 sample were used in subsequent analyses. PCA was performed using the R base function prcomp() with centered-log-ratio-transformed proportions via the R package “compositions” (111). OTUs detected at >0.1% mean proportion of the spring/summer 2005 sample reads were included in the PCA. This PCA allowed an exploratory examination of the beta diversity differences between individuals within the cohort. Cluster analysis of the spring/summer PC1 scores was performed with the mclust() R package (65). PC loadings from the spring/summer data were used to project the winter OTU centered-log-ratio-transformed proportions into the spring/summer data space generating projected PC1 and PC2 scores and to assign cluster based on the spring/summer fit for the winter data. The base R implementations of statistical tests (Fisher's exact test, Wilcoxon's rank sum test, and Welch's t test) were used where noted to identify enrichment of categorical variables.
Accession number(s).Sample FASTQ files are publicly available at the NCBI Sequence Read Archive under study accession number SRP058027.
ACKNOWLEDGMENTS
We thank the Yakima Valley community and families who participated in this study and the staff members from the Fred Hutchinson Cancer Research Center who performed sample collections.
We disclose that we have filed for a record of invention pertaining to the associated effects of azinphos-methyl on microbiota.
I.B.S., J.C.W., and E.M.F. conceived of the project and experiments. I.B.S., F.H.G., C.S.W., and S.H. designed experiments. I.B.S., C.S.W., F.H.G., J.T., M.K., and S.H. performed experiments. E.M.V., W.C.G., S.H., B.T., and E.M.F. curated and collected samples. I.B.S. wrote computer code and implemented analysis. I.B.S., J.C.W., T.W., E.M.V., W.C.G., A.S., and E.M.F. managed and analyzed data. I.B.S., E.M.F., J.S.M., C.S.W., A.S., and J.C.W. wrote the paper. All authors reviewed and commented on the paper.
FOOTNOTES
- Received 20 July 2016.
- Accepted 24 October 2016.
- Accepted manuscript posted online 11 November 2016.
Supplemental material for this article may be found at https://doi.org/10.1128/AEM.02149-16.
- Copyright © 2016 American Society for Microbiology.
REFERENCES
- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵
- 84.↵
- 85.↵
- 86.↵
- 87.↵
- 88.↵
- 89.↵
- 90.↵
- 91.↵
- 92.↵
- 93.↵
- 94.↵
- 95.↵
- 96.↵
- 97.↵
- 98.↵
- 99.↵
- 100.↵
- 101.↵
- 102.↵
- 103.↵
- 104.↵
- 105.↵
- 106.↵
- 107.↵
- 108.↵
- 109.↵
- 110.↵
- 111.↵
