ABSTRACT
Tyrophagus putrescentiae is inhabited by bacteria that differ among mite populations (strains) and diets. Here, we investigated how the microbiome and fitness of T. putrescentiae are altered by dietary perturbations and mite populations. Four T. putrescentiae populations, referred to as dog, Koppert, laboratory, and Phillips, underwent a perturbation, i.e., a dietary switch from a rearing diet to two experimental diets. The microbiome was investigated by sequencing the V1-V3 portion of the 16S rRNA gene, and selected bacterial taxa were quantified by quantitative PCR (qPCR) using group/taxon-specific primers. The parameters observed were the changes in mite population growth and nutritional status, i.e., the total glycogen, lipid, saccharide, and protein contents in mites. The effect of diet perturbation on the variability of the microbiome composition and population growth was lower than the effect induced by mite population. In contrast, the diet perturbation showed a greater effect on nutritional status of mites than the mite population. The endosymbionts exhibited high variations among T. putrescentiae populations, including Cardinium in the laboratory population, Blattabacterium-like bacteria in the dog population, and Wolbachia in the dog and Phillips populations. Solitalea-like and Bartonella-like bacteria were present in the dog, Koppert, and Phillips populations in different proportions. The T. putrescentiae microbiome is dynamic and varies based on both the mite population and perturbation; however, the mites remain characterized by robust bacterial communities. Bacterial endosymbionts were found in all populations but represented a dominant portion of the microbiome in only some populations.
IMPORTANCE We addressed the question of whether population origin or perturbation exerts a more significant influence on the bacterial community of the stored product mite Tyrophagus putrescentiae. The microbiomes of four populations of T. putrescentiae insects subjected to diet perturbation were compared. Based on our results, the bacterial community was more affected by the mite population than by diet perturbation. This result can be interpreted as indicating high stability of the putative intracellular symbionts in response to dietary perturbation. The changes in the absolute and relative numbers of Wolbachia, Blattabacterium-like, Solitalea-like, and Cardinium bacteria in the T. putrescentiae populations can also be caused by neutral processes other than perturbation. When nutritional status is considered, the effect of population appeared less important than the perturbation. We hypothesize that differences in the proportions of the endosymbiotic bacteria result in changes in mite population growth.
INTRODUCTION
Symbiotic intracellular bacteria have been traditionally classified as primary, including ancient and obligate bacteria, and secondary, including recent and facultative bacteria (1). The primary symbionts are those adapted to the host and vice versa. The associations between secondary symbionts and hosts are facultative although secondary symbionts are also adapted to their hosts. Moreover, secondary symbionts may be absent in individuals or populations (2). Variations in symbionts affect host structure and dynamics at the population level (3). The coexistence of infected and uninfected populations or the presence of populations with incompatible strains of symbionts indicates that symbionts form a metacommunity among the populations of insects (4). The endosymbiotic bacteria Wolbachia, Rickettsia, and Cardinium are present in 52, 24, and 13%, respectively, of all arthropod species (5). It was found that Wolbachia and Cardinium maintain prevalence within many arthropod taxa by implementing a variety of reproductive manipulations (6). Although Wolbachia and Cardinium are transmitted maternally, both taxa or their genes have been found to be transmitted horizontally among spider mites (7). More recently, horizontal Wolbachia transfer between Drosophila melanogaster stocks facilitated by Tyrophagus putrescentiae has been reported; moreover, this transmission occurred through the ingestion of Wolbachia in infested corpses (8). However, Rickettsia has not yet been identified in domestic mites (9, 10).
Approximately 3 decades ago, Bowman (11) demonstrated intraspecific differences in various populations of Acarus siro, Lepidoglyphus (Glycyphagus) destructor, Tyrophagus longior, and Rhizoglyphus sp. The intrapopulation differences in the physiological status of these mite populations in terms of their digestive enzymes were higher than the differences among species (11–13). Recently, differences in population growth were observed among four tested T. putrescentiae populations: three of the four populations exhibited density-dependent growth (14). One possible explanation for this discrepancy is the presence of different symbionts or differences in their proportions in mite populations. Analyses of the Sanger sequences of 16S rRNA genes from six T. putrescentiae populations support potential intraspecific physiological differences among populations of T. putrescentiae as a result of associated bacteria (15). The differences were caused by the endosymbionts Wolbachia, Cardinium, and Blattabacterium-like, and the putative symbiont Bartonella-like bacteria. Solitalea-like bacteria were present in all tested populations of T. putrescentiae (15) and are suggested as putative symbionts in A. siro (16). The proportion of bacteriocytes localized in the T. putrescentiae fat body (17) was low in the observed individuals and increased after perturbation (diet change) (15). This effect was attributed to the perturbation associated with different nutrient requirements and subsequent acquisition of bacteria (Brevundimonas vesicularis, Stenotrophomonas maltophilia, Serratia liquefaciens, and Serratia marcescens) from the environment as well as bacteriocyte formation under certain conditions (18, 19). These endosymbionts are not primary symbionts in T. putrescentiae mites because they are not present in all individuals (18, 19).
In this study, we examined whether the T. putrescentiae bacterial community is influenced more by perturbation (diet switch) or by mite population and whether there are correlations between fitness indicators and microbiome composition in T. putrescentiae populations. To address these questions, we analyzed four populations of T. putrescentiae mites. Population growth and nutritional status were analyzed as fitness parameters in samples of adult/juvenile mites. The bacterial community was characterized by sequencing the V1-V3 fragment of the 16S rRNA gene and performing quantitative PCR (qPCR) with group/taxon-specific primers.
RESULTS
The bacterial communities of the T. putrescentiae populations differed.The bacterial community (see Table S3 in the supplemental material) was more influenced by population than perturbation (Tables 1 and 2), as indicated by both the Illumina and qPCR data sets. According to permutational multivariate analysis of variance (PERMANOVA) and redundancy analyses (RDAs), regardless of whether the Bray-Curtis or Jaccard distances were utilized (P ≤ 0.001 for all cases), the F values clearly document this finding. The two-way PERMANOVA and distance-based RDA (db-RDA) results using the Jaccard index differed from the analyses of the Bray-Curtis matrix (Tables 1 and 2). For instance, the role of interactions between populations and perturbation was slightly lower in the Jaccard space but remained significant in terms of P values (Table 1). The presence/absence of operational taxonomic units (OTUs) was more influenced by population, but certain OTUs were influenced by perturbation, as indicated by the Venn diagrams (Fig. S3 and S4).
Comparison of the effect of populations and the tested variables and their interactions on T. putrescentiae microbiome compositiona
The effect of perturbation on the microbiome composition of T. putrescentiae populations tested separately by one-way PERMANOVA and db-RDAa
The similarities of the samples and contribution of OTUs were visualized by triplots (Fig. 1A and B) of the Illumina and qPCR data (see below for environmental variables and model selection). In both data sets, the most important OTUs distinguished among the population belonged to endosymbiotic bacteria (Cardinium, Wolbachia, Solitalea-like, and Blattabacterium-like bacteria) and to Bartonella-like bacteria. For the Illumina data set, OTUs of members of Bacillaceae and Enterobacteriaceae were important, i.e., Staphylococcus cohnii (OTU 4), Bacillus cereus (OTU 9), and Xenorhabdus (OTU 14). The random forest algorithm (errors = 0) analyses indicated that the following OTUs were responsible for differences: (i) OTUS of Micrococcaceae with low similarity (90%) to Acaricomes (OTU 15), Enterobacteriaceae with low similarity (92%) to Xenorhabdus (OTU 14), Cardinium (OTU 3), and Kocuria koreensis (OTU 19) were characteristic of the laboratory population; (ii) OTUs of Staphylococcus lentus (OTU 16) and B. cereus (OTU 9) were characterized as specific to the Koppert population; (iii) OTUs of Wolbachia (OTU 1), Bartonella-like (OTU 2), and Solitalea-like (OTU 6) bacteria were characteristic of both the Phillips and dog populations; and (iv) an OTU of Blattabacterium-like bacteria (OTU 8) was associated with the dog population only (Table S5). The associations were visible in the Krona projections (Fig. S2).
Triplot visualization of principal coordinates in redundancy analysis of the T. putrescentiae bacterial community. The analysis includes population, diet, and growth as environmental factors. A correlation triplot containing sample scores given by the weighted sums of OTUs was constructed. (A) Standardized Illumina data set transformed into Jaccard matrix. (B) qPCR data set transformed into a Euclidean matrix. The perturbation was coded as follows: 0, rearing diet, i.e., SPMd for all populations with the exception of the dog population (DDFd); 1, SPMd/DDFd or, for the dog population, DDFd/SPMd; 2, SPMd/OSRd or, for the dog population, DDFd/OSRd. Lab, laboratory population; Phil, Phillips population; Kop, Koppert population; P1, perturbation 1; P2, perturbation 2.
Perturbation of bacterial communities.A heat map (Fig. 2) was created to visualize the differences in the numbers of OTUs or taxa in the Illumina data set. There were differences in the log2-transformed data, indicating some contribution from low-abundance OTUs (Fig. S5 indicates cluster significance). When perturbation and population were incorporated into the random forest analyses, the error rate (0.222) was higher than for population only (Table S4). The most important OTUs were those of Micrococcaceae with low similarity to Acaricomes (OTU 15), Blattabacterium-like (OTU 8), S. cohnii (OTU 4), and Cardinium (OTU 3).
Heat map representing the distribution of samples and OTUs in the T. putrescentiae bacterial community in the standardized Illumina subsample data set. Data were log2 transformed. The samples are characterized by different colors, as follows: red, Koppert; blue, laboratory; green, Phillips; violet, dog population. The legend for the sample codes is provided in Table S2, and the legend for the OTUs is provided in Table S3.
The perturbation had different effects on the microbiome depending on the mite population. The changes were random rather than systematic in the qPCR data set (Fig. 3). The number of Wolbachia bacteria was increased by perturbation in the Phillips, Koppert, and dog populations but decreased in the laboratory population. In contrast, Cardinium and Solitalea-like bacteria were not influenced by perturbation in the above populations, apart from Solitalea-like bacteria in the Koppert population (Fig. 3).
Quantification of the bacterial community of T. putrescentiae based on qPCR with universal and taxon-specific primers (primers are given in Table S1). The numbers of copies were recalculated per single mite and log10 transformed for the universal and taxon-specific primers, as indicated. The perturbation was coded as follows: 0, rearing diet, i.e., SPMd for all populations with the exception of the dog population (DDFd); 1, SPMd/DDFd or, for the dog population, DDFd/SPMd; 2, SPMd/OSRd or, for the dog population, DDFd/ORSd.
Based on the Illumina data set, the laboratory population contained 13 core OTUs, and perturbation resulted in the addition of three OTUs: an OTU of Lactobacillus amylovorus (OTU 37) for perturbation 1 and two OTUs of Virgibacillus (OTUs 21 and 34) for perturbation 2. The dog population included nine core OTUs, and perturbation 1 resulted in a loss of three noncore OTUs (from L. amylovorus [OTU 37], Lactobacillus sp. [OTU 71], and Pantoea vagans [OTU 102]) and the acquisition of four OTUs (from Wolbachia sp. strain 2 [OTU 24], Prevotella paludivivens [OTU 45], Lactobacillus paracollinoides [OTU 106], and Pseudomonas putida [OTU 203]). The Phillips population included nine core species, and perturbation resulted in the disappearance of three OTUs that originally occurred with the control stored product mite diet (SPMd): P. paludivivens (OTU 45), L. paracollinoides (OTU 106), and Selenomonas lacticifex (OTU 202). Perturbation 1 provided three new OTUs: from L. amylovorus (OTU 37), Staphylococcus saprophyticus (OTU 38), and S. lentus (OTU 16). The Koppert population included 11 core OTUs, and perturbation resulted in the disappearance of OTUs 27, 30, and 31 (classified as Bacillaceae_2), whereas L. amylovorus (OTU 37), Lactobacillus sp. (OTU 71), and P. vagans (OTU 102) appeared after perturbation 1 (Fig. S4).
The observed F and P values resulting from the RDAs (Tables 3 and 4) suggest that the inverse Simpson diversity index was more influenced by perturbation than by population. The effect of the perturbation was apparent in all populations, but there were exceptions in the dog population (perturbation 1, dry dog food kernel diet switched to oilseed rape diet [DDFd/OSRd]) and Phillips population (perturbation 1, SPMd/DDFd). Alpha diversity appeared to be influenced more by perturbation, whereas the population showed a greater influence on beta diversity (Fig. 4).
Comparison of the effects of population and perturbation and their interaction on T. putrescentiae population growth and nutritional status
Comparison of the effects of population on T. putrescentiae growth and nutritional statusa
Comparison of the observed parameters of the T. putrescentiae populations under the control condition and after diet switch. The points are means; the bars are standard deviations. (A) Inverse Simpson diversity index (1/D) of bacterial community characterized by the Illumina data set. (B) Population growth, expressed as the logarithms of numbers of individuals after 21 days of growth from 10 individuals. (C to F) The nutrient contents, as indicated, were recalculated per milligram of nutrient to gram of fresh weight of mite. The perturbation was coded as follows: 0, rearing diet, i.e., SPMd for all populations with the exception of the dog population (DDFd); 1, SPMd/DDFd or, for the dog population, DDFd/SPMd; 2, SPMd/OSRd or, for the dog population, DDFd/ORSd.
Fitness parameters.Based on the results of RDA applied to normalized growth, population clearly had a greater effect than diet (P < 0.001 versus P = 0.015). The Koppert population differed somewhat from the others (Fig. 4). Perturbation 2 exhibited almost no effect in the Koppert population although its influence was clearly apparent in the other populations (Tables 3 and 4). The interactions between population and perturbation were also significant (P < 0.001). RDA also suggested a substantial influence of both population and perturbation on the overall distribution of the studied nutrients (Table 3). However, perturbation contributed more heavily to the mutual effects of these two factors. Indeed, the respective F values for interaction, population, and perturbation were 4.92, 3.5, and 11.51. All P values were less than 0.005. When the effects on the status of each nutritional component (i.e., glycogen, lipids, proteins, and saccharides) were studied separately (Table 3), all influences (including their interactions) were significant at least at the 0.05 level, with the exceptions of population in the case of glycogen status (F = 2.35, P = 0.07) and the interaction between population and perturbation in the case of protein contents (F = 1.38, P = 0.233) (Fig. 4).
Bacteria associated with nutrient contents.The combination of the forward selection of explanatory variables and variance inflation factor (VIF) values greater than 10 indicated that it is not advisable to add nutritional status and growth to the basic db-RDA model including population and diet together. However, when we focused on perturbation instead of the diets themselves, the VIFs decreased remarkably. In the case of the Illumina data set, lipid status alone could be added to the basic model obtained using the Jaccard distances. In the case of the qPCR data set, glycogen content showed similar results. Therefore, triplots were constructed with these variables (population, perturbation, and lipids for the Illumina data set and population, perturbation, and glycogen status for the qPCR data set) (Fig. 1A and B). According to the best RDA model, population was responsible for the highest variability in the Illumina data set (adjusted R2, cumulative, 0.8), followed by lipid status (adjusted R2, cumulative, 0.814) and perturbation (adjusted R2, cumulative, 0.829). Additionally, the analysis of variance (ANOVA)-like permutation test applied to the terms in the RDA model confirmed the significant effects of population, lipid status, and perturbation and the interaction between population and perturbation (all with P ≤ 0.01) on the OTU distribution. All populations differed from each other in terms of bacterial composition, as shown in the triplot visualization of principal coordinates in RDA (Fig. 1A). Lipids were closely related to perturbation 2. RDA models were constructed to describe the interactions between environmental variables and absolute numbers of 16S rRNA gene copies (qPCR data). The combination of the forward selection of the environmental variables with the VIF computation and the ANOVA-like permutation tests applied to the individual terms confirmed the inclusion of all populations, glycogen content, perturbation 2, and the interaction between population and perturbation in general (for all cases, P < 0.001). Based on the triplot (Fig. 1B) constructed for the model including all these mentioned variables, predominantly negative correlations were apparent for the qPCR data set. However, certain positive associations were also identified. Specifically, Cardinium was closely associated with the laboratory population and its interactions, along with perturbation 1. The vector of glycogen status formed very small angles with the vectors of Actinobacteria and Blattabacterium-like bacteria (and to some extent also with the vector of Wolbachia). Perturbation 2 explained a certain degree of variance among the other constituents of the microbiome.
DISCUSSION
Differences in bacterial communities among the mite populations.Previous analyses of insect species possessing various eating habits, i.e., omnivores, herbivores, and predators, have revealed significant influences of both diet and taxonomic origin on the spectrum of associated bacteria (20–23). This study examined four different populations (strains) of one astigmatid mite species. Our experiments were based on the assumption that the associated bacteria are influenced by mite origin (population) or different diets (perturbation). The observation that population is associated with microbiome composition indicates that population-specific factors, including genotype and metacommunity composition, influence the microbiome. According to the results of the RDAs, the different T. putrescentiae populations exerted a greater influence than perturbation on microbiome variability. We hypothesize that this result is caused by high stability of the symbionts in response to dietary perturbation in short-term adaptation. Once the microbiome is established, the symbionts appear relatively resilient to subsequent perturbations, such as diet changes. We cannot exclude different situations in the microbiome after long-term diet adaptation when the mite physiology is changed (24). Although this factor is a selective process, we cannot rule out neutral processes, including passive dispersion (25), during perturbation.
The endosymbionts Wolbachia, Cardinium, Blattabacterium-like, and Solitalea-like bacteria contributed to the explained variability in both the Illumina and qPCR data sets. Although Cardinium and Blattabacterium-like bacteria were detected in all populations by Illumina technology, they were under the detection limits for qPCR analysis (1 copy per mite). Therefore, we considered Cardinium and Blattabacterium-like absent, similar to the findings of a previous study in which we used 16S rRNA Sanger sequences (15). The next identified core taxa included Bartonella-like, S. cohnii, B. cereus, Brevibacterium siliguriense, and Enterobacteriaceae with low similarity to Xenorhabdus, which were all resistant to perturbation. Here, we also provide evidence that the mites acquire bacteria (Lactobacillus, Pantoea, Prevotella, and Virgibacillus) after perturbation, as previously observed in cultivation experiments (17, 18, 26). This is also supported by the beta diversity data according to the Jaccard index (Tables 1 and 2). The T. putrescentiae populations also differed in terms of the proportions of B. cereus, S. cohnii, Virgibacillus, Brevibacterium siliguriense, Micrococcaceae, and Enterobacteriaceae taxa. Analysis of a wide spectrum of insect species revealed the presence of Firmicutes as the main group of bacteria found in insect guts (23), which was confirmed in this study in the astigmatid mite T. putrescentiae. The presence of Enterobacteriaceae taxa in T. putrescentiae raises the question of whether these bacterial taxa are symbiotic in the mite. Previously, we observed sequences with 99% similarity to Xenorhabdus cabanillasii among clones from L. destructor and A. siro (10), whereas sequences similar to those of Erwinia bacteria comprised the majority of the microbiome of Blomia tropicalis (9). “Candidatus Erwinia dacicola” resides intracellularly in the gastric ceca of the larval midgut but extracellularly in the lumen of the foregut and the ovipositor diverticulum in adult flies (27). Moran et al. (28) described three symbiotic Enterobacteriaceae species, “Candidatus Serratia symbiotica,” “Candidatus Hamiltonella defensa,” and “Candidatus Regiella insecticola,” as sister groups to one another; moreover, they revealed a relationship to a species of Photorhabdus. These bacterial species occur in aphids and are localized to bacteriocytes, sheath cells, and the hemocoel (28). Based on these findings in the literature, the bacteria observed in this study may be located in both the gut and the parenchymal or reproductive tissues of the mite; future localization studies are necessary to confirm this assumption. The Enterobacteriaceae have been suggested to be beneficial based on their contribution to the longevity of the Mediterranean fruit fly (Ceratitis capitata) (29). The sequences obtained in the current analyses differed from those identified previously in the stored product mite Xenorhabdus (10). The Enterobacteriaceae sequences with low similarity (92%) to those of Xenorhabdus from T. putrescentiae should be identified in future studies using markers other than the 16S rRNA.
Differences in bacterial communities caused by perturbation and the effect on population growth.In insects, different nutritive factors, such as the protein-carbohydrate balance, influence the animal's life history, including reproduction (30, 31). Perturbation can lead to dysbiosis, i.e., a disruption of the proportions of gut core bacteria, and homeostasis, i.e., changes in the host physiology (32). In our experimental design, we used two initial (rearing) diets, SPMd for the majority of the T. putrescentiae populations and DDFd for the dog population, which might influence the results obtained. The perturbation factor overcame this issue, so the statistical analyses were not influenced by the difference in the initial diets.
In this study, we used population growth and increased nutrient contents in the mite samples as indirect indicators of fitness. Under suitable conditions, we assumed shorter development and a higher number of offspring than under unsuitable conditions, resulting in increased growth (33). All selected diets used in this study are natural for T. putrescentiae (34). The connection between diet perturbation and changes in the gut bacterial community is well documented (21, 23, 35). We observed that diet perturbation was associated with a rapid change in T. putrescentiae population growth, as previously reported (24), and we confirmed the strong effect of perturbation on both population growth and nutrient status. The effect of diet perturbation on the variability of the microbiome composition and population growth was lower than the effect induced by the mite population. In contrast, perturbation had a greater effect on nutritional status in mites than the mite population. Experiments with bacterial communities in Drosophila melanogaster have revealed that the bacterial community influences the relationship between energy storage and body mass and coordinates metabolism with body size (36). The RDA statistical model revealed the influence of lipid and glycogen status on the microbiome, which is correlated with Actinobacteria and Blattabacterium-like bacteria.
All populations contained Bartonella-like, Bacillus, Staphylococcus, Brenneria, or Kocuria bacteria. Perturbation resulted in changes in the numbers of these five taxa and the subsequent acquisition or disappearance of three or four taxa (Lactobacillus, Pantoea, and Virgibacillus). The Koppert population was associated with high numbers of Bacillaceae with similarity (90 to 99%) to Virgibacillus. The Sanger sequencing of the 16S rRNA previously identified sequences with similarity (94 to 99%) to Virgibacillus after a diet shift from SPMd to a fungal diet (37). These Bacillaceae include bacteria that produce enzymes facilitating the utilization of nutrients by the mites (18). Bacterial exoenzymes are active against difficult-to-digest structural polysaccharides such as cellulose and chitin (18) and are increased in the bacterial communities of T. putrescentiae on Fusarium diets (37). Bacillus cereus isolates originating from dog kernels infested by T. putrescentiae produce alkaline and neutral proteases (38). The participation of mite-associated bacteria in the digestion of skin, hair, nails, and feathers has been suggested (39) and could enhance the diet switch for T. putrescentiae.
A recent study by Serbus et al. (40) showed that the indirect effects of diet on Wolbachia titers in Drosophila ovaria include influencing cyclic AMP (cAMP)-regulated transcriptional coactivators in the TORC pathway (40). Experiments with the white fly Bemisia tabaci and the symbiont Hamiltonella showed that for insects reared on a standard-nitrogen diet, no cost or benefit was associated with Hamiltonella infection, whereas on low-nitrogen diets, Hamiltonella-infected whiteflies often grew better than uninfected whiteflies (41). In this study, diet perturbation influenced endosymbiotic bacteria, i.e., Blattabacterium-like and Wolbachia bacteria, whereas Solitalea-like and Cardinium bacteria were influenced marginally. However, the effect of perturbation on the absolute numbers of endosymbionts in the T. putrescentiae populations appeared to be random rather than systematic.
Multi-infection by Cardinium and Wolbachia is common at the individual level in spider mites (Acari: Tetranychidae) (42). At the population level, we observed different OTUs for Wolbachia in one population, in which low-abundance OTUs appeared after perturbation. In addition, the T. putrescentiae populations in our study were multi-infected with more taxa of endosymbiotic bacteria. The infection frequencies in the mite populations ranged from completely uninfected to a polymorphic state (43). Tetranychus truncates mites doubly infected with Cardinium and Wolbachia exhibit greater fecundity than uninfected individuals (44). Based on our data, we are unable to conclude whether perturbation causes dysbiosis in terms of the changes in the proportions of specific bacterial taxa in T. putrescentiae populations. Both the microbiome and population growth were more influenced by population than by diet perturbation, indicating the effect of endosymbionts or other bacteria on population growth. We hypothesize that the different proportions of Cardinium/Wolbachia or Blattabacterium-like/Wolbachia bacteria at the population level after perturbation could result in changes in population growth due to the various combinations of endosymbiotic bacteria in T. putrescentiae individuals. Similar changes have been observed in insects (2, 4) although the effects of these symbionts could also be the results of neutral processes (25). However, additional experiments are necessary to compare the effects of endosymbiotic bacteria on T. putrescentiae population growth.
MATERIALS AND METHODS
Mites and rearing diets.The bacterial communities of four populations of Tyrophagus putrescentiae (Schrank, 1781) were compared: (i) the laboratory population, which was sampled from a grain store in Bustehrad, Czechia, by E. Zdarkova in 1996; (ii) the dog population, which was sampled from dog food in the United States by J. Hubert in 2007; (iii) the Koppert population, which was obtained from a grain-derived diet by E. van Baal in the Netherlands in 2012; and (iv) the Phillips population, which was obtained from dog food by T. Phillips in Manhattan, KS, USA, in 2014. The mites were mass reared in Iwaki flasks and collected as previously described (45).
Diets utilized for perturbation.Except for the mites of the dog population, T. putrescentiae was reared on a wheat-derived diet designed for stored product mites (SPMd), which included a mixture of oat flakes, wheat germ, and Pangamin-dried yeast extract (Rapeto, Bezdruzice, Czechia) (10:10:1 wt/wt). The mites were transferred to the diets up to 6 months after collection and maintained until the experiments.
The dog population was maintained on a ground dry dog food kernel diet (DDFd) (Friskies Junior Life Plus Nutrition; Nestle Purina, Buk, Hungary) (see references 24 and 14). Both diets were heated at 70°C for 30 min and rehydrated. SPMd, DDFd, and the oilseed rape (Brassica napus L.) diet (OSRd) were compared in the dietary switch experiments. SPMd and DDFd are described in the previous section. The oilseeds were derived from a mixture of cultivars and were prepared according to a previously described protocol (46).
Treatments and samples.The experimental design included a rearing diet, SPMd, and two experimental diets, OSRd and DDFd. However, for the dog population, DDFd was the rearing diet, and OSRd and SPMd were the experimental diets. For the statistical analyses, the factor variable “perturbation” was coded in the following manner: 0, rearing diet, i.e., SPMd for all populations except the dog population (DDFd); 1, SPMd/DDFd for all populations or DDFd/SPMd for the dog population; 2, SPMd/OSRd for all populations or DDFd/OSRd for the dog population.
Unsexed mites (ca. 1,000 individuals) from approximately 1-month-old cultures were used in all experiments. Flasks containing mites and food were placed in Secador desiccator cabinets (Bel-Art Products, Wayne, NJ, USA) at 85% relative humidity and 25°C. The design included at least 10 chambers per diet and mite population. The mites that remained on the plug or the surfaces of the chambers were collected using a fine-tipped artist's paint brush. Then, the samples of mite populations were transferred into Eppendorf tubes and weighed on a microbalance (MS Mettler-Toledo, Greifensee, Switzerland); each sample had an accurate weight ranging from 0.015 to 0.025 g of wet weight. Multiple mites from each population were used for DNA extraction and the evaluation of protein, saccharide, lipid, and glycogen contents. DNA was extracted from the surface-cleaned and homogenized mites according to a previously described protocol (15). The extracted DNA was stored in a freezer at −28°C until analysis.
Growth test.The growth test began by switching 10 adult mites onto the two experimental diets and one rearing diet in Iwaki flasks. We performed 10 replicates per diet and mite population. The mites were maintained under the conditions described above. The experiment was terminated as described previously (38), and the mites were counted under a dissection microscope. To standardize the diverse levels of individual numbers associated with the diets, the vectors of the experimental diets were projected onto the vector of the rearing diet (SPMd or, in the case of the dog population, DDFd). The well-known procedure of vector projection is detailed, for instance, by Gentle (47). Then, a redundancy analysis (RDA) (48, 49) was applied to examine whether the perturbation exerted any influence or if the changes in growth were determined by populations. Using the vegan package (50) of the R statistical software (51), models with diets and perturbation were built. In addition, the interactions between population and perturbation were investigated. In particular, we aimed at changes possibly caused by perturbation within separate populations, where a special RDA was employed as well. Permutation tests (1,000 replications) were employed to determine the effect of the explanatory variables.
Nutritional status: protein, saccharide, lipid, and glycogen contents in mites.Using a method modified from Kaufmann (52), the mites were homogenized by beating in a FastPrep-24 homogenizer (MP Biomedicals, Santa Ana, CA, USA) in 200 μl of 2% sodium sulfate. Twenty-five microliters of each homogenate was used to estimate protein contents using the Bradford reagent (catalog number B6916; Sigma-Aldrich, St. Louis, MO, USA). The optical density (650 nm) was measured on a Secomam spectrophotometer (Secomam, Ales, France) according to the manufacturer's instructions. The reagent was first calibrated to a protein standard (catalog number P0834; Sigma-Aldrich). The remainder of the homogenates was precipitated in 800 μl of chloroform-methanol mixed 1:1 (vol/vol) and centrifuged (3,000 × g for 60 s). The pellets were retained for glycogen analysis. The supernatants were transferred to new Eppendorf tubes, and 600 μl of distilled water was added. After centrifugation (3,000 × g for 60 s), the top fraction (water-methanol) was collected for sugar analysis, and the bottom portion (chloroform) was used for lipid analysis. For lipid analysis, chloroform was evaporated on a heat block (110°C), and 200 μl of sulfuric acid (95 to 98%) was added for 10 min at 110°C. Next, 1 ml of vanillin-phosphoric acid reagent was added. The reagent contained 0.6 g of vanillin in 100 ml of hot distilled water and 400 ml of 85% phosphoric acid. The developed color was measured on a Secomam spectrophotometer at 450 nm. The reaction was calibrated to a lipid standard. For the sugar and glycogen assays, the tubes were placed on a heat block (110°C) to evaporate the solvent to a total volume of 100 μl, and 1 ml of anthrone reagent was added. The anthrone reagent contained 150 ml of distilled water, to which 385 ml of sulfuric acid (95 to 98%) and 0.75 g of anthrone were added. The samples containing the anthrone solution were mixed and heated for 17 min, and the optical density was then measured at 625 nm. The reaction was calibrated using a glucose standard (catalog number G6918; Sigma-Aldrich).
The data were recalculated to milligrams/gram of fresh weight. Two types of analyses were performed using R statistical software (51) and the vegan package (50). First, RDA models for individual nutrients were built similarly to those for growth as a response variable. Again, the interactions between the explanatory variables and changes within separate populations were studied as well. Second, an RDA model based on the correlation matrix computed from all nutrients together was utilized. Finally, permutation tests (1,000 replicates) were employed to determine whether population and perturbation exhibited important effects on the variability in nutritional status.
Amplification and sequencing.The design included four populations on three diets and six biological replicates (i.e., 72 samples). DNA samples were sent to MR DNA (Shallowater, TX, USA) for sequencing of the V1-V3 portions of the 16S rRNA genes using the universal primers 27Fmod and 519Rmod on an Illumina MiSeq platform, employing methods based on the bTEFAP (bacterial tag-encoded FLX amplicon pyrosequencing) process (53) in MR DNA.
Bioinformatic analyses.The multiplex sequences were processed using the mothur, version 1.36.1, software (54), according to the MiSeq standard operation procedure (SOP) (55) and UPARSE pipeline in USEARCH software (56). The mothur commands, which were used with default settings, are available online (http://www.mothur.org/wiki/MiSeq_SOP). Read lengths, both forward and reverse, were 300 bp. The sequences were demultiplexed using mothur. The list of samples and barcodes is provided in Table S2 in the supplemental material, and primers are given in Table S1. The sequences were assembled in UPARSE, and 6,927,385 sequences were obtained. The barcodes and primers were trimmed in mothur and fastq, and renamed fasta files were prepared. The renamed fasta file contained the sample name in the beginning of the first line. Then the both fastq and renamed fasta files were processed in UPARSE. Individual OTUs were constructed by binning sequences into clusters with 97% similarity, discarding putative chimeric OTUs (2.6%) and singletons (93.4% from unique sequences) in the process. The OTUs were classified according to Ribosomal Database Project training set 15 (57). Representative sequences for OTUs (from UPARSE) were processed using the blastn program on the NCBI platform (58) to identify closely related sequences in GenBank. We identified Cardinium, Solitalea-like, and Bartonella-like OTUs by comparison of the representative sequences to previously obtained, nearly full-length 16S rRNA cloned sequences from T. putrescentiae (15) (Table S3). The taxonomical features of the samples were visualized by performing KRONA projections (59).
The shared file was constructed and processed using mothur and later in PAST, version 3.06, software (60). The obtained sequences differed in the number of reads, which ranged from 25,071 to 152,889 (Fig. S1). The data standardization was based on the subsample created in mothur on 25,071 sequences. The subsampled data set (Table S4) was used for all further statistical analyses with Illumina data. The inverse Simpson index was calculated in mothur and then subjected to another RDA very similar to the previous one (see “Growth test” and the paragraph above on nutritional status). For Venn diagrams, OTUs were selected based on the criterion of presence in all replicates per treatment. Based on the selection, the Venn diagrams were constructed per population or treatment in mothur. The effects of population and diet on OTU distribution were tested by performing two-way and one-way permutational multivariate analysis of variance (PERMANOVA) (61) in the vegan package of the R platform, function adonis (50), using the Bray-Curtis and Jaccard distance matrices with 1,000 permutations. In the two-way analysis, we explored the interactions not only between population and diets (or perturbation) but also between population and growth and individual nutrients. In the one-way analysis, we studied how the diet or its switch influenced OTUs within every single population. All computations were performed using the vegan package in R (50). Furthermore, using the same package, several db-RDA models (occasionally referred to as models of the constrained analysis of principal coordinates) based on the Bray-Curtis or Jaccard dissimilarity matrices were constructed that included population and diet (eventually perturbation) as the explanatory variables. Fitness variables (growth and nutrients) were also added to explore whether they contribute to model improvement in terms of decreasing the P value. We simultaneously controlled the values of variance inflation factors (VIFs) (62) to avoid collinearity in the models. The R package packfor (63), which can indicate significant explanatory variables via the forward selection procedure, was used. To visualize the coordinates resulting from the best db-RDA models, triplots were created. Furthermore, a heat map depicting OTU abundance and clustering in dendrograms was produced using the gplots package in R (64). The OTU abundance data were logarithmically transformed as suggested by Anderson et al. (65). We prefer the log2 transformation, which permits better visualization of less abundant species. When constructing the dendrograms, we used the Bray-Curtis transformation. The similarity profile (SIMPROF) procedure (66) was applied to determine whether sample clusters were significant using the clustsig package in R (67).
A random forest algorithm with logarithmically transformed (log2) data and METASTATS (68) were applied to describe the effects of diets on OTUs for each population separately. The analyses were performed in mothur (1,500 trees and 1,000 replicates, respectively).
qPCR analyses.Amplifications were performed on a StepOnePlus real-time PCR system (Life Technologies, Grand Island, NY, USA) in 96-well plates using GoTaq qPCR master mix (Promega). SYBR green (Bio-Rad Laboratories, The Netherlands) was used as a double-stranded DNA (dsDNA) binding dye. The protocol and standard preparation were as described previously (16), and the primers and annealing temperatures are described in Table S1. The design employed six biological replicates, each analyzed with two technical repeats. To evaluate the effects of population and perturbation, the obtained numbers of copies were first logarithmically transformed (log10). Since the data were already transformed this way, we built on the Euclidean indices in subsequent analyses using the dissimilarity or distance matrices. Again, two-way and one-way PERMANOVA were performed together with db-RDAs. Whereas the two-way PERMANOVA technique provided specific insights regarding the relationships between interactions (population and perturbation), its one-way counterpart permitted the study of the influence of diets (or their switches) within each population separately. Using selected db-RDA models and indicating significance according to their P values obtained via the permutation tests (1,000 replicates), we mainly constructed triplots and visualized the situation. Experimentally, other explanatory variables such as growth and nutrients were added to the RDA models to determine whether they contribute to improvements in the P values.
Accession number(s).Sequences were deposited in the NCBI database under BioProject number PRJNA328005 and SRA number SRP078077.
ACKNOWLEDGMENTS
The authors are obliged to the referees for their highly useful comments on the previous version of the manuscript. We are grateful to Elmer van Baal and Tom Phillips for providing the T. putrescentiae populations, to Jan Kopecky for primer design and his advice and help with sequence processing, and to Martin Markovic for technical help.
This study was supported by project number GA15-09038S of the Czech Science Foundation.
FOOTNOTES
- Received 17 January 2017.
- Accepted 21 February 2017.
- Accepted manuscript posted online 24 February 2017.
Supplemental material for this article may be found at https://doi.org/10.1128/AEM.00128-17.
- Copyright © 2017 American Society for Microbiology.
