Research ArticleMICROBIOLOGY

Coadaptation between host genome and microbiome under long-term xenobiotic-induced selection

See allHide authors and affiliations

Science Advances  05 May 2021:
Vol. 7, no. 19, eabd4473
DOI: 10.1126/sciadv.abd4473

Abstract

One of the most difficult experimental challenges today is testing the evolutionary dynamics shaping complex host-microbiome interactions. We investigated host-microbiome codiversification in response to xenobiotic-induced selection using an experimental evolution approach. To this end, we exposed the parasitoid wasp Nasonia vitripennis to sublethal concentrations of the widely used herbicide atrazine for 85 generations. Our results reveal that atrazine exposure not only mediated adaptive changes within the microbiome, which conferred host resistance to atrazine toxicity, but also exerted selective pressure on the host genome and altered host gene expression and immune response. Furthermore, microbiome transplant experiments reveal a decreased survival of adults from the control population after exposure to the evolved microbiome of the atrazine-exposed population, while no such decrease occurred in the reciprocal transplant. These results indicate that xenobiotic-induced selection mediated host-microbiome coadaptation, ultimately leading to a new host genome–microbiome equilibrium.

INTRODUCTION

All healthy animals and plants are colonized by many microorganisms, notably bacteria (1, 2). This symbiotic microbiome plays a notable role in host physiology and evolution. Specifically, the evidence that the microbiome modulates host development (25), health (6), nutrition (79), immunity (1012), and even speciation (1317) illustrates how the complex interactions between a host and its microbiome determine individual phenotypes. Furthermore, phylosymbiotic signatures (i.e., phylogenetic congruence between closely related host species and their associated microbiome) in mammals, insects, and corals (1721) support the existence of intimate and long-lasting host-microbiome relationships in diverse organisms. As a consequence, hosts and their associated microbiome are increasingly regarded as interconnected organisms, or “holobionts” (a host and its associated microbial community), forming an evolutionary unit (22, 23). This interconnectedness between the different members of a holobiont makes it difficult to investigate how these complex symbiotic systems evolve. Specifically, it remains challenging to disentangle the relationships between the multiple partners within the symbiotic system, as each microbe may interact not only with the host but also with multiple other microorganisms, resulting in intricate interaction patterns (24, 25). Consequently, gaining a better understanding of the evolutionary dynamics shaping host-microbiome interactions and the resulting phenotypes is one of the most crucial challenges in studying evolutionary biology today.

The microbiome represents a source of phenotypic and evolutionary novelty for metazoan hosts (26). For instance, the microbiome may be able to respond more rapidly to new selective pressure, e.g., changing environmental conditions, thereby facilitating the adaptation of the holobiont by providing a fitness benefit. There are numerous reported examples of rapid adaptive changes within the microbiome leading to a fitness benefit within different model systems. The bean bug Riptortus pedestris became resistant to an insecticide via the acquisition of a new symbiotic bacterium able to degrade the insecticide (27). Similarly, the increase of Actinobacteria, the community composition of the plant root microbiome, is implicated in the drought resistance of several plant species (28). These observations illustrate the potential contribution of the microbiome to host fitness via its capacity to adapt to selective pressure. However, drastic microbiome changes could also be problematic for the host, e.g., resulting in dysbiosis. For instance, a mismatch between host genotype and microbiome composition plays a role in interspecies hybrid lethality in Nasonia wasps (17). To avoid adverse fitness effects, the host needs to tolerate the rapid changes in the microbiome, e.g., through changes in regulatory mechanisms such as immunity. Hence, an external selective pressure might act on the host in two ways: (i) directly, since the host is exposed to the selective pressure, and (ii) indirectly via the changes in its microbiome.

Using the parasitoid wasp Nasonia vitripennis (Hymenoptera: Chalcidoidea), a tractable insect model to experimentally investigate host-microbiome interactions (3, 17, 29, 30), we recently demonstrated that continuous exposure to the toxic herbicide atrazine resulted in the enrichment of several atrazine-metabolizing bacteria within the wasp microbiome after only a few insect generations (31). The altered microbiome provided host resistance to atrazine (31). In the present work, we demonstrate that both the host and its microbiome respond to xenobiotic-induced selection, ultimately leading to a new host-microbiome equilibrium resulting from the host genome–microbiome coevolution. We experimentally evolved the microbiome of N. vitripennis by subtoxic atrazine exposure over 85 insect generations and monitored the impact of long-term atrazine exposure on microbiome composition and function as well as on host physiology and genetic differentiation. Our results confirmed that long-term atrazine exposure mediates adaptive changes within the microbiome, as the evolved microbiome metabolized atrazine more efficiently, which conferred resistance to atrazine toxicity to the host. In addition, long-term atrazine exposure also induced selective pressure on the host genome and altered host gene expression and immune response. Last, microbiome transplant experiments revealed a negative fitness effect when exposing adults from the control population to the evolved microbiome of the atrazine-exposed population. These results indicate that both the host and its microbiome adapted to atrazine exposure, either independently or in synergy.

RESULTS

Atrazine exposure mediates adaptative microbiome evolution in Nasonia

Two experimental populations were established from the laboratory line AsymCx of N. vitripennis: (i) a control population reared on sucrose solution (CCc) and (ii) an atrazine population (AAa) reared on sucrose solution supplemented with a subtoxic concentration of atrazine, 30 parts per billion (ppb) (Fig. 1A). We previously reported that atrazine exposure provokes changes in the microbiome that conferred resistance to atrazine to the wasp host (31). Here, 16S ribosomal RNA (rRNA) gene amplicon sequencing of the microbiome of both populations over 85 generations confirmed that the AAa microbiome gradually diverged from the parental CCc1 microbiome based on weighted UniFrac index, while the CCc microbiome remained relatively stable compared with the parental CCc1 microbiome over time (Fig. 1B). Using quantitative polymerase chain reaction (qPCR), we observed that bacterial density was higher in the AAa population compared to the CCc population at generation F35 but dropped back to control levels in later generations (fig. S1A) (31). Moreover, two bacterial strains (Serratia marcescens NVIT01 and Pseudomonas protegens NVIT02), which conferred host resistance to atrazine (31), were found to be enriched in the AAa population (fig. S1, B and C). Determining the median lethal atrazine concentration (LC50) at 19 time points throughout the 85 generations indicated that the AAa population developed resistance to atrazine, as it survived better to higher atrazine concentrations compared to the control population from the eighth generation onward (Fig. 1C). We confirmed that atrazine resistance was conferred by the microbiome, since feeding wasps from the AAa population an antibiotic sucrose diet (see Materials and Methods) at different generations eliminated atrazine resistance (Fig. 1D) and altered microbiome composition (fig. S1D), without other fitness costs on wasp longevity (fig. S1E). Likewise, axenic wasps from the AAa population obtained after germ-free (GF) rearing lost their atrazine resistance (Fig. 1D). These results support the conclusion that the microbiome of N. vitripennis plays a crucial role in atrazine detoxification and that this function is enhanced in the altered microbiome of the AAa population.

Fig. 1 Adaptive microbiome evolution mediates atrazine resistance in N. vitripennis.

(A) Experimental design. Populations were fed control or atrazine diet for 85 generations (F85). At F25, populations were split to establish subpopulations. (B) Linear regression analysis of microbiome dissimilarity. The microbiome of the atrazine population became increasingly divergent compared with the parental CCc1 microbiome. n = 11 (CCc1), 10 (AAa8), 9 (CCc8), 7 (AAa12), 6 (CCc12), 8 (AAa20), 8 (CCc20), 16 (AAa30), 14 (CCc30), 8 (AAa31), 7 (CCc31), 7 (AAa32), 8 (CCc32), 8 (AAa35), 8 (CCc35), 8 (AAa36), 8 (CCc36), 8 (AAa45), 8 (CCc45), 8 (AAa51), 8 (CCc51), 10 (AAa58), 8 (CCc58), 6 (AAa69), 7 (CCc69), 8 (AAa75), and 7 (CCc75). (C and D) The mean LC50 value increased across generations in the atrazine-exposed population while remaining constant in the control population (C), and atrazine resistance was lost in wasps fed an antibiotic and in GF-rearing wasps (D). Different letters indicate significant differences based on the resistance ratio 95% CIs, which does not include 1.0. (E) Microbiome dissimilarity based on weighted UniFrac (WU) distance. n = 8 (AAa35), 8 (CCc35), 6 (AAc35), 8 (CCc36), 8 (AAa45), and 8 (CCc21). Scatterplot shows means with SD. ****P < 0.0001; ns, not significant; Student’s t test.

After demonstrating that atrazine exposure can alter the wasp microbiome and that the altered microbiome conferred a fitness advantage to the host, we next investigated how both microbiome structure and the induced host resistance evolved when atrazine exposure was removed. After 24 generations, both the CCc and AAa populations were split into two subpopulations, and the diet was switched for one of the subpopulations, such that one of the control subpopulations was fed the atrazine diet (CCa) and one of the atrazine subpopulations was fed the control diet (AAc) (Fig. 1A). As previously observed for the AAa population, once the CCc population was switched to the atrazine diet, it exhibited an increase in atrazine resistance (Fig. 1C). In addition, the two bacterial strains S. marcescens NVIT01 and P. protegens NVIT02 known to play a role in conferring this resistance were enriched in the microbiome of the CCa population 10 generations after the switch to the atrazine diet (fig. S1, F and G), and the susceptibility to atrazine could be restored by treatment with antibiotics (Fig. 1D). In the AAc population, atrazine resistance was observed for at least 17 generations after switching from the atrazine to the control diet (i.e., until F41), but resistance was completely lost 27 generations after the diet switch (i.e., at F52) (Fig. 1C). The observed resistance after atrazine removal was still microbiome dependent, as AAc wasps fed an antibiotic diet lost their atrazine resistance (Fig. 1D). Moreover, the AAc microbiome composition changed over time: 11 generations after the switch to the control diet (F35), the microbiome was more similar in structure to the AAa microbiome than to the CCc microbiome, but this was no longer the case 21 generations after the diet switch (F45) (Fig. 1E). Despite the gradual change in AAc microbiome composition, it is important to note that the microbiome did not revert back to the microbiome of the control population (Fig. 1E). Together, these results demonstrate that exposure to atrazine mediated a change in the microbiome of N. vitripennis, which, in turn, conferred an adaptive advantage to the host in the presence of atrazine. The acquired atrazine resistance was heritable for more than 11 generations after the removal of atrazine, thus providing a continued benefit to the host when exposed to atrazine during this period.

Changes in the AAa microbiome increase atrazine detoxification efficiency

To gain an in-depth understanding of how the microbiome mediated host resistance to atrazine exposure in the AAa population, we performed an untargeted metabolomic analysis by liquid chromatography–mass spectrometry (LC-MS) and compared metabolite profiles of the CCc, AAa, and their respective axenic GF populations at F85 after feeding on a 3-ppm (parts per million) atrazine diet. Sucrose-fed wasps of each population were included as controls. The mortality rate of the CCc population and the GF populations was significantly higher than that of the AAa population 12 and 24 hours (fig. S2, A and B, respectively) after exposure to atrazine. No difference in mortality rate between the populations was observed in sucrose-fed controls, indicating that the increased mortality was due to atrazine toxicity.

We detected 8573 metabolites in conventionally reared wasps and 1106 metabolites from the GF populations (table S1). Unexpectedly, only 0.94% (82 of 8753) of the metabolites were shared between the GF and conventionally reared populations, indicating the importance of the gut microbiome for Nasonia’s metabolism. Furthermore, the CCc and AAa populations, which harbored distinct microbial communities, also had distinct metabolic profiles (Fig. 2A), while their GF counterparts’ metabolomes were indistinguishable (Fig. 2B). This suggests that the difference between the CCc and AAa metabolomes was mainly mediated by their different microbiome.

Fig. 2 Microbiome changes affect the metabolome of N. vitripennis and increase atrazine detoxification efficiency.

(A) Principal components analysis (PCA) plot of total metabolomics between control and atrazine-exposed populations with 95% confidence ellipses. n = 54 for conventional wasps. (B) PCA plot of total metabolomics between control and atrazine-exposed GF populations. n = 58 for GF wasps. (C) Desethylatrazine metabolite concentration in conventionally reared populations after 12 and 24 hours of exposure to control or atrazine diet. Four to six biological replicates for each condition. (D) Desisopropylatrazine metabolite concentration in the conventionally reared populations after 12 and 24 hours of exposure to control or atrazine diet. Five to six biological replicates for each condition. *P < 0.05, **P < 0.01, and ****P < 0.0001, Student’s t test.

Next, we compared the metabolites detected in the control versus atrazine-fed groups of each experimental population. Two metabolites associated with metabolic processing of atrazine were of particular interest: Desethylatrazine, a secondary metabolite of atrazine (fig. S2C) (32), was elevated after 12 hours of exposure to 3-ppm atrazine diet in both the AAa (978-fold increase compared to sucrose diet) and CCc populations (333-fold increase compared to sucrose diet) (Fig. 2C) but was absent in GF wasps since the GF wasps cannot metabolize atrazine themselves (Fig. 2C). Desisopropylatrazine, another secondary atrazine metabolite (fig. S2C) (32), also increased significantly in the atrazine-exposed AAa and CCc populations (2.2-fold and 1.6-fold increase compared to sucrose-fed controls, respectively) (Fig. 2D) but was not detected in GF wasps since the GF wasps cannot metabolize atrazine themselves. After 24 hours of atrazine exposure, both desethylatrazine and desisopropylatrazine were still elevated in the CCc population (1292-fold and 5-fold increase compared to sucrose-fed controls, respectively) (Fig. 2, C and D), indicating that atrazine was still present in these wasps. In contrast, both compounds were only marginally detectable in the AAa population (Fig. 2, C and D), suggesting that the pesticide had been completely detoxified in the AAa population. These results indicate that AAa wasps harboring the evolved microbiome were more efficient in the metabolic detoxification of atrazine compared to CCc wasps harboring the ancestral microbiome.

Signatures of atrazine-induced selection in the host genome

Next, we investigated the impact of atrazine on the host genome. Using microsatellite analysis, we detected an increased rate of genetic divergence in the atrazine-exposed population (31). Here, we used whole-genome resequencing to investigate genetic divergence between the CCc and the AAa populations. To this end, we sequenced the genomes of 12 to 24 individuals from five populations: the ancestral population CCc1 and both populations at generation 35 (CCc35 and AAa35) and generation 70 (CCc70 and AAa70) (table S2). The fixation index (Fst) was slightly higher between the ancestral population CCc1 and the atrazine population AAa70 (Fst value of 0.0042) than between CCc1 and the derived control population CCc70 (Fst value of 0.0031) (table S2), which indicates a gradual genetic divergence in both populations. The low genome-wide divergence is not unexpected given the short time scale and high genetic homogeneity of the laboratory population (33). A maximum likelihood estimation of individual ancestries using ADMIXTURE (34) provided support for two clusters (fig. S3A), and we showed clusters from two to five (fig. S3B).

We identified candidate regions under atrazine-induced selection throughout the genome by identifying regions with maximum Fst values and minimum genetic diversity index (Pi) values. This approach identified 343 genes, representing a total of 1.7 Mb of genomic sequence, with maximum Fst values (top 10%) and minimum Pi1/Pi2 values (top 10%) in the AAa70 population compared to the ancestral CCc1 population (Fig. 3A). These candidate selective regions were detected when comparing both the AAa35 and AAa70 populations to the CCc1 population (table S2) but were not present in the same genomic regions when comparing the CCc35 or CCc70 population to the CCc1 population (table S2). This suggests that the candidate selective regions identified in the AAa population may not have been caused by random genetic drift but by selection either from atrazine or from the altered microbiome after atrazine exposure. The functions of most genes in the conserved selective regions are unknown (table S2).

Fig. 3 Genomic regions with selective sweep signals in the AAa population.

(A) SNPs in the genomes of the AAa70 population identified based on the Fst values and θπ ratios compared to the CCc1 ancestral population. n = 24 and 12 for CCc1 and AAa70, respectively. (B) Example of LOC100119293 with strong selective sweep signals in the genomes of the AAa35 and AAa70 populations. Fst values and θπ ratios are plotted using a 10-kb sliding window. The red line and blue line represent the 10% outlier threshold of Fst and Pi values, respectively. The gray bar indicates the selective region.

Atrazine exposure alters host immunocompetence and immune gene expression

Given the crucial role of the host’s innate immune system for health and in mediating host-microbiome interactions, we next investigated host immunity differences after long-term atrazine exposure using immune challenge experiments and RNA sequencing (RNA-seq). Immune challenge experiments were performed on female pupae of the CCc, AAa, and AAc populations at generation F69. By this time, the AAc microbiome no longer conferred the same level of atrazine resistance to the host (Fig. 1C). Yellow female pupae were inoculated either with sterile phosphate-buffered saline (PBS) (control), heat-inactivated Gram-positive bacteria (Staphylococcus lentus strain CICCHLJQ29 isolated from N. vitripennis), or heat-inactivated Gram-negative bacteria (S. marcescens strain NVIT01 isolated from N. vitripennis) using sterile insect needles dipped in the respective inoculum. Heat-inactivated bacteria are frequently used to elicit an immune response through the recognition of bacterial cell wall or surface components such as peptidoglycans or lipopolysaccharides, allowing to investigate the effects of immune activation while avoiding pathological effects due to bacterial proliferation (35, 36). The proportion of eclosed adults 120 hours after challenge was compared between the populations as an indicator of immune activation, reasoning that an induced immune response would channel energy resources away from metamorphosis. While no difference in adult eclosion was observed after inoculation with PBS, the proportion of eclosed adults was significantly different between the CCc and AAa populations after bacterial challenge, depending on the type of bacteria (Fig. 4A). While the proportion of eclosed adults was lower in the AAa population compared to the CCc population after challenge with heat-inactivated S. lentus, the opposite was observed after challenge with heat-inactivated S. marcescens NVIT01 (Fig. 4A). The latter was the same strain that can metabolize atrazine and was already enriched in the microbiome of the AAa population (31). Hence, it is conceivable that AAa wasps were better adapted to this strain than CCc wasps after 69 generations of atrazine-induced selection and thus did not mount an immune response against it. These results indicate differences in the immune response against two types of bacterial pathogens after long-term exposure to atrazine. The AAc population showed an intermediate eclosion proportion in both cases (Fig. 4A), indicating a long-lasting impact of atrazine exposure on host immunity, detectable even 45 generations after the last exposure to the herbicide.

Fig. 4 Atrazine exposure alters the host’s immune response and immune gene expression.

(A) Proportion of eclosed adults 120 hours after inoculation of female pupae with PBS (control) or heat-inactivated Gram-positive or Gram-negative bacteria. Different letters indicate significant difference (P < 0.05) based on ANOVA followed by Tukey’s multiple comparisons test. n = 4 for each treatment for every condition. (B) The global transcriptional response of the CCc and AAa populations 6 and 24 hours after immune challenge with Gram-positive or Gram-negative bacteria. Gene expression was normalized on the basis of PBS-inoculated wasps as controls. (C) Upset plot showing the overlap of DEGs between parental and F1 hybrid populations 6 hours after Gram-positive bacteria challenge. (D) Immune gene expression 6 and 24 hours after bacterial challenge in parental and F1 hybrid populations. (E) Upset plot showing the overlap of DEGs between parental and F1 hybrid populations 6 hours after Gram-negative bacteria challenge.

The atrazine selection–induced difference in immune response was further investigated using RNA-seq on pools of five female pupae (three to four biological replicates per condition) from the F69 generation. RNA-seq was performed at two time points: 6 and 24 hours after challenge with either heat-inactivated Gram-positive or Gram-negative bacteria, compared to inoculation with PBS as a control. Globally, challenge with the Gram-positive S. lentus induced a weaker transcriptional response than challenge with the Gram-negative S. marcescens in both the CCc and AAa populations (Fig. 4B). The transcriptional response to the same bacteria was very different between the two populations: Both populations had a similar number of differentially expressed genes (DEGs) 6 hours after challenge with heat-inactivated Gram-positive bacteria compared to PBS inoculation (91 DEGs in the CCc population and 114 in the AAa population) (Fig. 4B and table S3), but not a single DEG was shared between the two populations (Fig. 4C). Twenty-four hours after heat-inactivated Gram-positive bacteria challenge, the number of DEGs remained stable in the CCc population (96) but increased in the AAa population (367) (Fig. 4B and table S3). Only 48 DEGs (13.08% of the DEGs in the AAa population and 50.00% of the DEGs in the CCc population) were shared between the two populations (fig. S4A), indicating that the transcriptional response was still quite divergent, due to a high number of DEGs specifically regulated in AAa wasps (Fig. 4B).

Among the DEGs after heat-inactivated Gram-positive bacteria challenge, only six were related to immunity based on Gene Ontology (GO) analysis: two antimicrobial peptides (Defensin 1-1 and Defensin 1-2) and four peptidoglycan recognition proteins (PGRPs; LOC100119736, LOC100117452, LOC100121609, and LOC100121589) (Fig. 4D and table S3). Defensin 1-2 and the three PGRPs were significantly up-regulated in AAa wasps 6 hours after the heat-inactivated Gram-positive bacteria challenge, while none of them were differentially expressed in the CCc population at this time point (Fig. 4D and table S3). Twenty-four hours after the heat-inactivated Gram-positive bacteria challenge, the two antimicrobial peptides Defensin 1-1 and Defensin 1-2 as well as two PGRPs (LOC100119736 and LOC100117452) were up-regulated in both the CCc and AAa populations. In contrast, the other two PGRPs (LOC100121609 and LOC100121589) were up-regulated only in CCc wasps (Fig. 4D and table S3). These findings suggest a temporal shift in immune gene expression after Gram-positive bacteria challenge between the two populations since AAa wasps up-regulated immunity-related genes earlier than CCc wasps (Fig. 4D and table S3).

After a challenge with the heat-inactivated Gram-negative S. marcescens, the transcriptional response was generally stronger, especially in the AAa population (Fig. 4B). There were 437 and 229 DEGs in the CCc population 6 and 24 hours after heat-inactivated Gram-negative bacteria challenge, respectively (Fig. 4B and table S3), compared to 648 and 609 DEGs in the AAa population (Fig. 4B and table S3). Of these, 131 DEGs (20.22% of the DEGs in the AAa population and 29.98% of the DEGs in the CCc population) were shared between the two populations 6 hours after challenge (Fig. 4E). Similarly, 112 DEGs (18.39% of the DEGs in the AAa population and 48.91% of the DEGs in the CCc population) were shared between the two populations 24 hours after challenge (fig. S4A). Notably, the same two antimicrobial peptides Defensin 1-1 and Defensin 1-2 and the three PGRPs were up-regulated in both populations 6 and 24 hours after Gram-negative bacteria challenge (Fig. 4D and table S3). A third antimicrobial peptide (Defensin 2) and another PGRP (LOC100121589) were up-regulated only in CCc wasps 6 and 24 hours after Gram-negative bacteria challenge, respectively (Fig. 4D and table S3).

These results demonstrate the impact of long-term atrazine exposure on the global transcriptional response as well as on the expression of specific immunity-related genes within the first 24 hours after challenge with heat-inactivated bacterial pathogens. Independent of the type of bacteria injected, AAa female pupae had more DEGs (notably down-regulated genes) than CCc female pupae under most experimental conditions (Fig. 4B). Moreover, immune genes were up-regulated earlier in AAa female pupae compared to CCc female pupae after heat-inactivated Gram-positive bacteria challenge (Fig. 4D). It is noteworthy that AAa female pupae demonstrated a reduced survival to adulthood after challenge with heat-inactivated Gram-positive S. lentus (Fig. 4A). Still, the precise link between the observed transcriptional changes after bacterial challenge and the reduced survival 5 days later remains to be established. These results suggest profound changes in host physiology due to atrazine exposure, which may be due to one or more of these factors: (i) the presence of atrazine or its degradation products in the host’s body, (ii) selection of particular host genotypes in the atrazine-exposed population, or (iii) the altered microbiome influencing other aspects of host physiology as well, including the immune response.

Immune response and immune gene expression in F1 hybrids

Given that deregulation of the innate immune response is one of the notable factors mediating hybrid breakdown between different Nasonia species (17), we next investigated the immune response in F1 hybrids of the CCc and AAa populations. At generation F68, AAa females were mated with CCc males (population F1AC), and CCc females were mated with AAa males (population F1CA) (see below; Fig. 5A). No developmental defects were observed in the hybrids (fig. S5A), confirming the absence of any hybrid breakdown due to incompatibilities between the CCc and AAa host genotypes. However, the proportion of eclosed adults following the same immune challenge experiments in female pupae as above was significantly reduced in both F1 hybrid populations compared to the pure parental populations (Fig. 4A), indicating a stronger developmental impact of immune challenge in F1 hybrids.

Fig. 5 Host genome x microbiome interactions.

(A) The experimental crossing procedure to generate F1 hybrid populations. Black line, CCc genotype; red line, AAa genotype. (B) Microbiome dissimilarity based on weighted UniFrac distance showed that the microbiome is maternal-like. n = 10 (F1CC), 10 (F1CA), 11 (F1AA), and 11 (F1AC). Violin plot shows means with SD. Different letters indicate significant differences (P < 0.05) based on ANOVA followed by Tukey’s multiple comparisons test. (C) The LC50 value of parental and hybrid populations. Different letters indicate significant differences based on the resistance ratio 95% CI that does not include 1.0. (D and E) GO functional enrichment analysis of DEGs in female pupae 6 hours after PBS challenge. DEGs between F1AC and F1AA populations (D) and between F1CA and F1CC populations (E). The functional enrichment analysis was performed using the functional annotation tool DAVID. The innate immune response is in bold.

Gene expression analysis 6 hours after bacterial challenge further demonstrated differences in gene expression between the hybrids and the parental populations. First, the F1AC population had a higher number of DEGs than the F1CA population, independent of the type of bacteria injected (table S3): 182 versus 92 DEGs after heat-inactivated Gram-positive bacteria challenge and 197 versus 70 DEGs after heat-inactivated Gram-negative bacteria challenge (Fig. 4, C and E). This might indicate a maternal effect on gene expression, as AAa wasps generally had a higher number of DEGs than CCc wasps (Fig. 4B). Considering that both hybrid populations should have a similar genetic background (50% AAa and 50% CCc), the differences in gene expression could also be due to the microbiome, which is maternally inherited (see below). Second, only a low number of DEGs were shared between the hybrids and the parental populations (Fig. 4, C and E), suggesting divergent transcriptional responses in the hybrids. Notably, none of the seven immunity-related genes that were up-regulated in the parental populations after bacterial challenge were differentially expressed in the F1CA population (Fig. 4D and table S3). However, the expression level of these genes tended to be higher than that in the parental populations, even after PBS challenge (fig. S4B), indicating a higher constitutive expression level of these genes in the F1CA population. In contrast, the antimicrobial peptide Defensin 2 was up-regulated in the F1AC population after Gram-positive bacteria challenge, and all three Defensins were up-regulated after Gram-negative bacteria challenge (Fig. 4D and table S3). This indicates highly divergent immune gene expression patterns in both F1 hybrid populations compared to the parental populations. Considering that the host genetic background should be similar, this may be due to maternal effects or altered host genotype–microbiome interactions in the hybrids.

Host genome–microbiome interactions under long-term herbicide exposure

To investigate potential interactions between host genotype, gene expression, and the microbiome, we compared the female pupal microbiome of the two F1 hybrid populations F1AC and F1CA established at generation F69 (Fig. 5A; see above) with the microbiome of the parental populations at generation F69. The microbiome of both hybrid populations was maternal-like (Fig. 5B), due to maternal transmission of the microbiome (17), resulting in two hybrid populations with a similar genetic background (50% AAa and 50% CCc) but different microbiome. In line with this, the F1AC population was resistant to atrazine similar to the maternal AAa population, while the F1CA population was not, similar to the maternal CCc population (Fig. 5C).

Despite the absence of developmental defects in the hybrids (fig. S5A; generation 50 in fig. S5B), the observed differences in global transcriptional response, immune gene expression, and decreased survival to adulthood after bacterial challenge (Fig. 4, A and C to E) share similarities with the deregulated immune response resulting from host genotype–microbiome incompatibilities in interspecies hybrids between N. vitripennis and Nasonia giraulti (17). Functional enrichment analysis identified innate immune response as an enriched function among genes that were differentially expressed in F1AC female pupae (i.e., harboring an AAa-like microbiome) compared to AAa female pupae (Fig. 5D). The antimicrobial peptide Defensin 1-2 and three PGRPs (LOC100119736, LOC100117452, and LOC100121609) were constitutively up-regulated in F1AC compared to the AAa population, even in the absence of bacterial challenge (fig. S4B and table S3). In contrast, the innate immune response was not enriched among DEGs in F1CA female pupae (i.e., harboring a CCc-like microbiome) compared to CCc female pupae (Fig. 5E), with only one antimicrobial peptide Defensin 1-1 being constitutively up-regulated in F1CA (fig. S4B and table S3). This demonstrates asymmetric differences in constitutive immune gene expression between the hybrid and parental populations harboring the same microbiome. Notably, a higher number of immunity-related genes were constitutively up-regulated in the hybrid population harboring the evolved AAa microbiome, which may indicate genotype × microbiome interactions between the ancestral CCc genotype and the evolved AAa microbiome. On the other hand, the two hybrid populations did not differ in immune gene expression despite their different microbiome (fig. S4B and table S3). This may indicate that genotypic differences between the hybrid and parental populations partly mediate the observed differences in gene expression as well.

We also detected six DEGs in the hybrid populations compared to the parental populations that were located in the selective regions identified in the genomes of the AAa population (table S3). Of these, five were differentially expressed between the F1CA and CCc populations (three up-regulated and two down-regulated in the hybrids), and one was down-regulated in F1AC hybrids compared to the AAa population. In addition, another gene was differentially expressed between the two hybrid populations (fig. S6A and table S3). Only three of the seven genes were differentially expressed between the two parental generations (fig. S6A), indicating that these changes in gene expression were not mainly due to atrazine-induced selection but occurred specifically in one of the hybrid populations. This suggests that host genotype x microbiome interactions may influence the regulation of these genes in the hybrid genetic background.

Microbiome transplants validate host genome–microbiome interactions

The interaction between host genotype and the altered microbiome after atrazine exposure was experimentally tested using intra- and interpopulational microbiome transplant experiments. To this end, axenic adults from the CCc and AAa populations (F69) were obtained using GF rearing. One-day-old adults were exposed either to sucrose solution as a control or to sucrose solution supplemented with microbiome from either the CCc or AAa population (F69) at two different concentrations. We observed that the survival of CCc adults decreased significantly after an exposure to the higher dosage of the AAa microbiome, while no such decrease occurred after exposure of AAa adults to the ancestral CCc microbiome (Fig. 6). These results demonstrate incompatibility between the evolved AAa microbiome and the CCc genotype when the two are put in contact without a period of coadaptation. The fact that the reciprocal transplant (CCc microbiome in AAa genetic background) did not produce negative fitness effects suggests that AAa wasps recognize the CCc microbiome as their species-specific microbiome. This asymmetric host genome–microbiome incompatibility mirrors the asymmetric differences in constitutive immune gene expression between the F1 hybrids and the parental population harboring the same microbiome (Fig. 5, D and E).

Fig. 6 Survivorship after intrapopulational and interpopulational microbiome transplants.

There was a significant increase in mortality in the CCc population after exposure to a high concentration of the AAa microbiome (log-rank test, P < 0.001).

DISCUSSION

This work used an experimental evolution approach to investigate the impact of an environmental selective pressure (the herbicide atrazine) on an insect host and its associated microbiome. This experimental design allowed us to identify selection and adaptation signatures at the level of the holobiont and the hologenome (22, 23). Using the parasitoid wasp N. vitripennis as a model system, we demonstrated that both the host and its microbiome responded to xenobiotic-induced selection. Regarding the microbiome, our results confirm that long-term atrazine exposure mediates adaptive changes within the microbiome after continuous exposure to subtoxic atrazine levels. The altered (or evolved) microbiome provided a functional benefit to the host as it metabolizes atrazine more efficiently, thereby conferring resistance to atrazine toxicity to the host. This is a tractable example of how microbiome adaptation can contribute to host (or holobiont) fitness by responding rapidly to selective pressure. Our results also showed that the altered microbiome was heritable over multiple generations, even after the removal of atrazine. The initial adaptation in response to atrazine happened quite fast (within eight generations) and lasted for more than 17 generations after the removal of atrazine. While the mechanisms underlying the subsequent loss of atrazine resistance are unknown, we hypothesize that it might be due to modified host regulatory mechanisms (e.g., gene expression and immunity) maintaining the altered microbiome or microbial interactions within the microbiome itself. Specifically, the enriched taxa might slowly decrease in abundance after the removal of atrazine, opening up new niches that can be recolonized by other members of the microbiome.

Our experimental evolution approach further demonstrated the impact of long-term atrazine exposure on the host at the physiological, transcriptional, and even genomic level. For instance, we identified candidate genomic regions with signatures of selection in the atrazine-exposed population. This may be due to direct atrazine-induced selection pressure on the host or reflect genomic adaptations to the altered microbiome after atrazine exposure. However, the data presented here do not allow us to disentangle the two scenarios.

In addition, we detected profound changes in the global transcriptional response and immunocompetence after challenge with bacterial pathogens. Atrazine did not systematically act as an immune suppressor, as previously demonstrated in vertebrate models (37, 38). Specifically, long-term atrazine exposure decreased host immunocompetence against Gram-positive bacteria but increased immunocompetence against Gram-negative bacteria. Two aspects are important to note in this context: First, N. vitripennis responded more strongly to Gram-negative bacteria even in the control population, indicating innate differences in recognition of different bacterial pathogens. After Gram-negative bacterial challenge, the global transcriptomic response was much stronger in the atrazine-exposed population even beyond immunity-related functions, indicating a global impact of atrazine on host physiology and metabolism. Second, the Gram-negative bacterial strain used for the immune challenge was one of the enriched atrazine-metabolizing taxa within the Nasonia microbiome (S. marscescens NVIT01). This may have played a role in the improved survival after immune challenge with this bacterium in the atrazine-exposed population, either through an increased host resistance to this bacterial species or through immune priming before bacterial challenge, mediated by the increased abundance of this bacterium in the microbiome. Furthermore, we observed an asymmetric host genome–microbiome incompatibility, i.e., reduced fitness after transplanting the evolved microbiome back into the control population, indicating that a period of adaptation is necessary for the host to tolerate the evolved microbiome. This is similar to the incompatibility between host genotype and microbiome involved in interspecies hybrid lethality in Nasonia wasps (17).

In summary, this work demonstrates adaptation to selective pressure at the hologenomic level, in that both the host and the microbiome responded to selection. The next challenge consisted of determining whether both partners responded independently to atrazine-induced selection or whether the adaptation to atrazine in one partner, in turn, mediated selection on the other partner. In other words, the rapid adaptive changes in the microbiome may have exerted selective pressure on the host to adapt to and tolerate the new microbiome, resulting in the observed genomic and transcriptional changes. Hybridization and microbiome transplant experiments between the atrazine-exposed and the control populations revealed an asymmetric incompatibility between the genetic background of the control population and the evolved microbiome of the atrazine-exposed population. This indicates that coadaptation was necessary for the host to tolerate the altered microbiome and that this coadaptation might have mediated (at least partly) the differential expression and fixation of certain host genes involved in microbiome regulation. Our work provides insights into the evolutionary dynamics of holobionts, demonstrating that selection may act at multiple levels: on the host, the microbiome, and both partners in synergy, ultimately leading to a new host-microbiome equilibrium.

MATERIALS AND METHODS

Chemicals and media

Atrazine analytical standard (CAS no. 1912-24-9, ≥97% purity) was purchased from TCI America (Portland, USA). Lysogeny broth (LB), nutrition broth, ampicillin (CAS no. 69-52-3), gentamicin (CAS no. 1405-41-0), kanamycin (CAS no. 25389-94-0), and penicillin-streptomycin were purchased from Sigma-Aldrich (St. Louis, USA). Other chemicals used in this paper were analytical grade reagents, and those used for high-performance liquid chromatography analyses were of LC-MS grade. Preparation of LB solid medium included the addition of agar (16 g/liter). The pH of media was regulated to 7.2 before sterilization. Sterilization conditions were 115°C for 30 min.

Nasonia rearing

Some of the following methods are similar to those previously published (31). The N. vitripennis strain AsymCX(u) (Wolbachia free) (33) was used to establish experimental populations. The AAa population were fed 30 ppb of atrazine in a 10% (w/v) sucrose solution while the CCc population were fed a control diet of 10% sucrose. Wasps were fed for 48 hours, then the food was removed, and the adults were provided with 25 Sarcophaga bullata pupae (fly hosts) for 48 hours. Within each lineage, all pupae were collected, mixed, and randomly sorted into groups of 50 females and 15 males for the next generation. The effective population size (Ne) is 42 based on Ne = 9 NmNf/(4Nm + 2Nf) (39). All N. vitripennis were reared in 25°C incubators with constant light. The S. bullata fly pupae were from the Seth R. Bordenstein laboratory (Vanderbilt University, USA) and Carolina Biological Supply (Burlington, NC). All fly host color and firmness were checked before providing them to adult N. vitripennis. The generation time for N. vitripennis under these rearing conditions is approximately 2 weeks. After 25 generations (F25), each linage was split into additional sublineages, and the diet was switched for subsequent generations.

GF rearing

Some of the following methods are similar to those previously published (31). All GF rearing was conducted following the protocol set forth by Brucker and Bordenstein (17) with slight modifications; host S. bullata fly pupae were first strained through a 100-μm mesh before centrifugation, no cell culture media or antibiotics were used, and volume differences were made up using sterile water. GF media was stored at 4°C for 2 weeks maximum.

Microbiome preparation

Some of the following methods are similar to those previously published (31). Microbiome was purified from L4 larvae by homogenization of larvae in sterile 1× PBS. The larval homogenate was then centrifuged at 1000 rpm for 3 min to remove large cellular debris, and the resulting supernatant was filtered through a 5-μm filter. The filtrate was centrifuged at 6000 rpm for 5 min, and the supernatant was removed. The pellet was resuspended in sterile PBS. The optical density (OD) of the bacterial solution was measured on a Multiskan FC microplate photometer (catalog no. 51119000, Thermo Fisher Scientific) using 620-nm filters (OD620).

Bacteria culture

Bacterial strains were grown in LB medium at 30°C at 250 rpm in baffled shake flasks for 16 hours and later centrifuged at 6000 rpm for 5 min. The resulting cell pellet was resuspended in the same volume of PBS media and centrifuged at 6000 rpm for 3 min. The resulting cell pellet was resuspended in PBS for the immune challenge. The OD of the bacterial cultures was measured on a Multiskan FC microplate photometer (catalog no. 51119000, Thermo Fisher Scientific) using OD620. Bacteria were heat-inactivated at 95°C for 10 min.

Antibiotic feeding

Mixture of penicillin-streptomycin (1000 μg/ml), gentamicin (1500 μg/ml), and kanamycin (20,000 μg/ml) was added into the sugar solution or atrazine food. The feeding protocol is identical with that in Nasonia rearing above.

Fecundity assay

Fifty virgin females and 15 virgin males were provided with 20 fly hosts in one vial to allow for mating and female ovary development. After 24 hours, batches of four females were transferred to new vials and provided with four new fly hosts for oviposition. After 24 hours, the females were removed, and offspring was counted throughout development: eggs (within 12 hours after oviposition), L4 larvae (6 days after oviposition), yellow pupae (10 days after oviposition), and adults (17 days after oviposition). Vials that contained dead females after the 24-hour oviposition period were discarded to insure that each batch had four female foundresses.

Measurement of LC50

Some of the following methods are similar to those previously published (31). Twenty to 50 yellow/black pupae were collected in a vial. After 24 hours, all dead pupae and those that had not yet emerged were discarded. The remaining newly emerged adults were fed varying concentrations of atrazine (0, 3, 10, 30, 100, 200, 300, 600, 800, and 900 ppm) for conventionally reared wasps and (0 ppb , 300 ppb, 3 ppm, and 10 ppm) for GF-reared wasps. At least two replicates were established per condition. Mortality was recorded after 0, 24, 48, 72, 96, 120, 144, and 168 hours following atrazine exposure. A concentration-response curve was generated and the best fit determined by probit analysis using Polo Plus-PC (40). Resistance ratios were calculated by dividing the LC50 value of the AAa population by the LC50 value of the same generation’s CCc population using Polo Plus-PC software version 2.0. Confidence intervals (CIs) for resistance ratios were calculated using the method described by Robertson et al. (40). Resistance ratios were compared across conditions to test the significance of resistance ratios at a 95% level CI. A difference between compared values is considered significant if the 95% CI does not include 1.0 (40).

Immune challenge assay

Yellow pupae (11 days after hosting) were sterilized by rinsing three times in 10% bleach and three times in PBS (30 s each time). Pupae were then allowed to air-dry for 1 hour at room temperature. Twenty-four pupae were stuck on double-sided tape in a small petri dish. Insect pins were sterilized by rinsing once with 10% bleach, once with 70% ethanol (EtOH), and once with PBS. The insect pins were then used to make a small hole in the abdomen of the pupae, and PBS or heat-inactivated bacteria solution was then injected into the hole using a micropipette. The petri dish containing the injected pupae was placed in a pan with some water to maintain humidity, and the pan was kept in the Nasonia rearing incubator as above. Adults and dead pupae were counted 5 days after immune challenge. The eclosion proportion was calculated as follows: Eclosion proportion=Number of adultNumber of initial pupae *100.

Microbiome transplantation

GF wasps were reared as outlined above. One-day-old adults were fed either sugar solution or sugar solution supplemented with microbiome from CCc or AAa population. Mortality was recorded after 0, 24, 48, 72, 96, 120, 144, and 168 hours following microbiome exposure.

Sample collection and DNA extraction

Some of the following methods are similar to those previously published (31). The founders for each generation of N. vitripennis were collected and maintained at −80°C until DNA and RNA extraction. For DNA extraction, each individual was put in a 1.5-ml tube and rinsed once with 1 ml of 70% EtOH and 1 ml of 10% bleach and twice with 1 ml of sterile water. Samples were then frozen and homogenized in liquid nitrogen. DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. DNA was quantified using the Qubit double-stranded DNA (dsDNA) HS assay kit on the Qubit 2.0 fluorometer (Life Technologies).

PCR, library preparation, and sequencing for microbiome analysis

Some of the following methods are similar to those previously published (31). A portion of the 16S rRNA gene was PCR-amplified from 4 μl of DNA using the 27F and V4R (5′-GGACTACHVGGGTWTCTAAT-3′) primers. Duplicate reactions per sample were performed using NEBNext High-Fidelity 2X PCR Master Mix in a total reaction volume of 20 μl, 50°C annealing temperature, and 25 cycles. The resulting PCR products were used for a nested PCR reaction using modified dual-indexed primers (18). PCR was performed using 2 μl of the previous PCR product in a 20-μl reaction using the NEBNext High-Fidelity 2X PCR Master Mix, with only 12 cycles. PCR products were purified using Agencourt AMPure XP beads, quantified using the Qubit dsDNA HS assay kit (Life Technologies, Waltham, USA) on the Qubit 2.0 fluorometer (Life Technologies, Waltham, USA), and equimolar amounts of each sample were pooled together. Each pooled library was run on the Illumina MiSeq using the MiSeq Reagent Kit V3 for paired-end reads. Sequencing was performed by the Dana-Farber Cancer Institute, Harvard Medical School.

Microbiome analysis

Some of the following methods are similar to those previously published (31). Raw Illumina sequence reads were processed using QIIME 2 v.2019.4 (41). Sequences were joined using the VSEARCH plugin (41) and filtered using Deblur (41). Open-reference operational taxonomic unit (OTU) picking was performed using VSEARCH and clustered at 99% identity using SILVA’s 16S rRNA QIIME database (v.132) (42). Low-abundance OTUs (frequency < 10) were filtered and removed across all samples. QIIME2’s Naïve Bayes classifier was trained using the SILVA 16S rRNA database (v.132, 99% consensus taxonomy, seven levels) using the q2-feature-classifier (41). The full-length 16S-trained classifier was used to generate taxonomic classification using the classify-sklearn feature-classifier plugin (41).

The taxa plugin was used to remove all unassigned, chloroplast, mitochondria, and Wolbachia reads from the dataset (taxa filter-table command). Any library containing less than 1000 reads was removed, and the remaining libraries were normalized by rarefying the sequencing depth to correct for uneven sequencing depth. Representative sequences for each OTU were used for alignments with MAFFT (alignment MAFFT plugin) (41), with unconserved and highly gapped regions masked (alignment plugin, mask command). Maximum likelihood phylogenetic trees were generated with FastTree 2 (phylogeny plugin) and midpoint rooted (phylogeny plugin, midpoint-root option). Following the generation of OTUs, taxonomy assignment, gene alignment, and diversity metrics were generated for each sample using QIIME 2’s diversity plugin (core-metrics command, phylogenetic option). Simpson’s index was calculated separately using the diversity plugin (command alpha with Simpson metric).

qPCR analysis

Some of the following methods are similar to those previously published (31). qPCR was performed using a StepOnePlus Real-Time PCR system (Thermo Fisher Scientific, Waltham, USA). All primers are listed in table S4. Each reaction contained 3 μl of template, 10 μl of GoTaq qPCR Master Mix (Promega, Madison, USA), 1 μl of primer mix (10 mM), and 6 μl of sterile water for a total volume of 20 μl. A no-template control was included in each run to check for reagent contamination. A melting curve analysis was performed after each run to confirm the amplification specificity. The thermal protocol was as follows: 40 cycles at 95°C for 10 s, 55° or 60°C for 30 s, and 72°C for 20 s. Each sample was run in duplicates.

To prepare standard solutions for the quantification of target genes, PCR amplicons of the target gene were purified using Agencourt AMPure XP and cloned into Invitrogen pCR2.1-TOPO vector with TOPO10 One Shot chemically competent cells. The plasmids were then prepared with the Qiagen Plasmid Mini Kit (Qiagen, Hilden, Germany) and quantified with the Qubit dsDNA HS assay kit on the Qubit 2.0 fluorometer. Tenfold dilution series from 107 to 103 copies were prepared for each standard and used to calculate the gene copy number.

Metabolomic sample preparation

Before feeding and 12 and 24 hours after feeding, wasps were flash-frozen in liquid nitrogen and then stored at −80°C until extraction. Five individuals were pooled for each biological replicate for conventionally reared wasp, and an individual wasp was used for each biological replicate for GF-rearing wasp. Spearman correlation was used to exclude the outlier samples (correlation r > 0.8 for conventionally reared wasps and correlation r > 0.7 for GF-reared wasps). Extraction was performed in bead beater vials containing 1 ml of methanol and three 2-mm steel beads. Samples were homogenized for 15 min at 50 Hz in a TissueLyser LT (Qiagen). The homogenates (including wasp fragments) were mixed with 1 ml of methanol and 4 ml of chloroform in 8-ml glass vials and incubated for 5 min in an ultrasound bath. Two milliliters of water was added to each sample, followed by another 5-min incubation in an ultrasound bath. The samples were then centrifuged at 2000g for 10 min, and the upper aqueous phases were transferred to new vials. The samples were dried under N2 flow before being resuspended in 100 μl of 70% acetonitrile and transferred into glass microinserts.

Metabolomic data acquisition and statistical analysis

Untargeted LC-MS metabolomics was performed at the Harvard Center for Mass Spectrometry. Samples were analyzed on a Thermo ID-X mass spectrometer, coupled to a Vanquish LC (both Thermo Fisher Scientific). Five microliters of samples was injected on a Zic-pHILIC column (150 mm by 2.1 mm, 5-μm particle, Sigma-Aldrich). The solvents used were as follows: A, 20 mM ammonium carbonate in water with 0.1% ammonium hydroxide, and B, acetonitrile 97% in water. The gradient was as follows: from 100% B to 40% B in 20 min, then to 0% B in 10 min, and followed by 5 min at 0%. The column was then re-equilibrated by bringing B back to 100% in 5 min and kept at 100% B for 10 min. The flow rate was 0.2 ml/min. A pool sample was created by combining a small fraction of each sample. The pool sample was run every 10 samples as quality control (QC). All samples were measured in polarity switching as MS1 scans, at 120,000 resolution. AquireX Deep Scan was used on the pool sample to acquire MSMS data. Two sets of deep scan (one for each polarity) were run, with five-sample deep coverage.

Compound discoverer (CD version 3.1, Thermo Fisher Scientific) was used to extract the data from the raw files. The data were analyzed in two batches (the five wasp samples and one wasp sample separately). The QC samples were used to normalize each individual compound and compensate for instrumental drift. Compound identification was obtained from comparing the MSMS scan with the MZCloud library as well as internal MZvault libraries. All compounds with library match were manually inspected for proper integration and good match with the library entries. The data were gap-filled to allow proper statistics using CD gap-filling node. A final median-centering normalization was applied to the data to compensate for small biomass differences. The significance was calculated by multiple t test, and the limit of significance was set at a fold change of >1.2 (up-regulated) or <0.83 (down-regulated) and P < 0.05. We used the R random forest package to perform the feature selection, and the MeanDecreaseGini value shows the importance of the feature.

Whole-genome resequencing DNA extraction and library preparation

For DNA extraction, each individual was put in a 1.5-ml tube and rinsed once with 1 ml of 70% EtOH and 1 ml of 10% bleach and twice with 1 ml of sterile water. Samples were then frozen in liquid nitrogen with additional mechanical homogenization. DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. Samples were quantified using the Qubit dsDNA HS assay kit on the Qubit 2.0 fluorometer (Life Technologies). The libraries were prepared at the Bauer Core Facility at Harvard University using KAPA HyperPrep Plus and enzymatic fragmentation. The DNA input amount was 20 ng. Illumina adapters containing dual 8–base pair (bp) indices (Integrated DNA Technologies Inc.) were attached to the DNA fragments using six PCR cycles.

Molecular marker discovery

Low-quality reads (i.e., less than 90% of all bases having a quality score of ≥30) were trimmed or discarded using fastp software (43). Trimmed reads were mapped to the N. vitripennis reference genome Nvit_2.1 (GCF_000002325.3) using BWA MEM (44) with default parameters. Molecular markers were called using the GATK Best Practices pipeline, which considers indel realignment and mark duplication, and calls variants across all samples simultaneously through the Haplotype Caller program in GATK 3.8.0 (45). Variants were filtered using standard hard filtering parameters according to GATK Best Practices pipeline. Specifically, single-nucleotide polymorphism (SNP) and insertions and deletions (Indels) were obtained by mapping quality of >37 and quality depth of >24. Last, variants with more than 70% call rate and sequence depth more than twofold were used for downstream analysis.

Selective region detective

A maximum likelihood tree was constructed using the raxml (version 8.0) (46) software on the basis of a distance matrix, using the whole-genome SNPs shared by all the accessions. The population structure was calculated using the Bayesian clustering program admixture (34). The average pairwise divergence within a population (θπ), the Watterson’s estimator (θw), and Tajima’s D value were estimated for all genomes of each population using a sliding window of 10 kb. In each window, these parameters were calculated using the vcf-tools (47), and Fst value was calculated to measure the differentiation between the two populations.

RNA extraction and RNA-seq library preparation

For RNA extraction, five individual wasps were pooled and crushed to fine powder under liquid nitrogen. Total RNA was extracted using the RNeasy Plus Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. RNA yield was quantified using the Qubit RNA HS assay kit (Life Technologies, Waltham, USA) on the Qubit 2.0 fluorometer (Life Technologies, Waltham, USA), and RNA quality was verified using a NanoDrop 2000/2000c spectrophotometer (Life Technologies, Waltham, USA). The RNA-seq libraries were prepared at the Bauer Core Facility at Harvard University. Starting with 500 ng as input, poly(A) selection and stranded library preparation were performed using the KAPA mRNA HyperPrep system. The library target insert size is 200 to 300 bp.

RNA-seq analysis

All samples were run multiplexed on an Illumina NovaSeq 6000 (paired-end, 150-bp reads) at the Bauer Core Facility at Harvard University (Cambridge, MA). Raw reads were trimmed using TrimGalore-0.6.0 (https://github.com/FelixKrueger/TrimGalore) and mapped to the N. vitripennis genome Nvit_2.1 (GCF_000002325.3) using STAR 2.7.2b (48). Table S3 provides a summary of reads per sample and the percentage of mapped reads. The gene raw read count and expression statistics are provided in table S3. Cufflinks v2.2.1 (49) was used to assemble the reads into transcripts. The fragments per kilobase of transcript per million mapped fragments were used to normalize the expression levels for different genes and transcripts. Raw read counts were imported into DESeq2 (50) to determine DEGs. A gene with a fold change of >1.2 (up-regulated) or <0.83 (down-regulated) and unadjusted P < 0.01 was considered to be a DEG. GO enrichment analysis and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis were conducted using the functional annotation tool DAVID 6.7.

Validation of RNA-seq gene expression via qPCR

The expression profiles of six transcripts that were differentially expressed in the AAa population based on RNA-seq data were verified using qPCR on the same samples that had been used for RNA-seq library preparation. Complementary DNA (cDNA) was synthesized using the Verso cDNA Synthesis Kit (Invitrogen) from 1 μg of total RNA as input template after treatment with deoxyribonuclease I (Invitrogen). Specific primers for the target transcripts were designed using Primer5 (http://premierbiosoft.com/primerdesign/). Gene expression was normalized using the delta CT method against the 60S acidic ribosomal protein P1 (60sp1) as a reference gene (table S4 for primer) (51). The qPCR protocol was identical to the standard protocol for qPCR analysis outlined above. Differences in gene expression between populations were tested using Mann-Whitney U test. The comparison of qPCR and RNA-seq is shown in (fig. S6B) and demonstrates the concordance between RNA-seq and qPCR results.

RNA-seq saturation analysis

To validate that the RNA-seq sequencing depth was sufficient for gene expression analysis, we performed a sequencing saturation analysis on the sample (88) for which we had obtained the lowest number of reads (table S3). Figure S6C shows that the sequencing depth was sufficient even for this sample, as the curve indicating gene discovery rate reaches a plateau, indicating sequencing saturation.

Statistical analysis

General. Some of the following methods are similar to those previously published (31). Data were considered statistically significant when P < 0.05, for log-rank test, Mann-Whitney U test, Student’s t test, and permutational multivariate analysis of variance (PERMANOVA) test, as indicated in the figure, figure legend, or experimental methods. Asterisks denote corresponding statistical significance: *P < 0.05; **P < 0.01; ***P < 0.001. Asterisk indicates the resistance ratio calculated by 95% CI that does not include 1.0 using Polo Plus-PC software for LC50. Data are presented as the means ± SD, means ± SE, and LC50 ± 95% CI where appropriate from at least three independent biological replicates, unless stated otherwise in the figures, figure labels, or experimental methods. Statistical analysis was performed using GraphPad Prism 7 software or R, PERMANOVA test in QIIME2, as indicated.

Quantification of gene expression. For RNA-seq analysis, a gene with a fold change of >1.2 (up-regulated) or <0.83 (down-regulated) and unadjusted P < 0.01 was considered to be a DEG.

Quantification of metabolomic expression. For metabolomic analysis, a metabolite with a fold change of >1.2 (up-regulated) or <0.83 (down-regulated) and P < 0.05 (multiple t test without correction) was considered to be differentially regulated.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/19/eabd4473/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank A. Shultz (Harvard University) for assistance in the early analysis of the whole-genome resequencing dataset. We thank A. T. Tate (Vanderbilt University), R. Harris (James Madison University), M. Levy (University of Pennsylvania), K. D. Kohl (University of Pittsburgh), N. Cira (Rowland Institute at Harvard University), B. M. Berdy (Broad Institute of MIT and Harvard), and E. O. Martinson (University of New Mexico) for reviewing early drafts of the manuscript and in the preparation of this manuscript. We thank C. Vidoudez from the Harvard Center for Mass Spectrometry (Harvard University) for providing the metabolomic raw data report. We thank C. Stokes (Rowland Institute at Harvard) for building a tool for the immune challenge assay and C. Reardon (Harvard University Bauer Core Facility) for assistance with the RNA-seq library preparation. We thank Y. Liu from Oebiotech for assistance with metabolomic analysis. We thank E. Ottum and M. Ali supported by Harvard summer program to assist experiments. Funding: This work was supported by the Rowland Institute at Harvard Fellowship awarded to R.M.B. Author contributions: R.M.B. and G.-H.W. conceived the study. R.M.B. supervised the study and wrote the manuscript. G.-H.W. performed experiments, analyzed the data, and wrote the manuscript. J.D. wrote the manuscript and assisted with analysis. B.D. and L.H. analyzed the data and assisted the manuscript writing. All authors contributed to the final version of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The RNA-seq raw reads, whole-genome resequencing, and the 16S rRNA data were deposited in the Sequence Read Archive (PRJNA623128) of the National Center for Biotechnology Information. Additional data associated with this paper have been deposited at Mendeley Data at doi: 10.17632/mtwxm75prd.1. The code supporting the current study is available from the corresponding author on request.

Stay Connected to Science Advances

Navigate This Article