**Patterns of Diversity in the Symbiotic Mite Assemblage of the Mountain Pine Beetle,** *Dendroctonus Ponderosae* **Hopkins**

**Sneha Vissa 1,\* , Javier E. Mercado <sup>2</sup> , Danielle Malesky <sup>3</sup> , Derek A. Uhey <sup>1</sup> , Boyd A. Mori <sup>4</sup> , Wayne Knee <sup>5</sup> , Maya L. Evenden <sup>6</sup> and Richard W. Hofstetter <sup>1</sup>**


Received: 1 October 2020; Accepted: 14 October 2020; Published: 17 October 2020

**Abstract:** The mountain pine beetle, *Dendroctonus ponderosae* (Coleoptera: Scolytinae)*,* is an economically important bark beetle species with a wide geographic range spanning from the southwestern United States into northern Canada. This beetle causes extensive tree mortality to 13 pine species. Mites (Acari) are common and abundant symbionts of mountain beetles that may influence their fitness through positive and negative interactions. We present a unique assessment of the mite associates of mountain pine beetles using measures of alpha and beta diversity. We sampled phoretic mites from five beetle populations: Arizona, Colorado, South Dakota, Utah (USA), and Alberta (Canada) that varied in host tree species, local climate, and beetle population level. We collected 4848 mites from 8 genera and 12 species. Fifty to seventy percent of beetles carried mites in flight with the highest mite loads occurring in middle and southern populations; decreasing in northern populations. Mite assemblages (i.e., both richness and composition) varied along a south to north latitudinal gradient and were driven by species turnover (i.e., species replacement). Differences in mite composition increased with distance between populations. We discuss climatic variation, environmental filtering, and host tree differences as factors that could affect differences in mite composition between beetle populations and discuss implications for functional shifts. Our results could represent a model for estimating diversity patterns of mite symbionts associated with other major insect pests in coniferous forest systems.

**Keywords:** biodiversity; bark beetles; symbionts; species assemblage; beta diversity; forest ecosystems

### **1. Introduction**

Globally, forests face increased pressure from insect pests, such as bark beetles (Coleoptera: Scolytinae) [1–3]. Bark beetle impacts can be severe as they actively damage tree phloem tissue and vector phytopathogens, resulting in large events of tree mortality [4,5]. The mountain pine

beetle, *Dendroctonus ponderosae* Hopkins, has been responsible for tree mortality across more than 30 million hectares of forest in western North America in the last decade [4,5]. Previously constrained by temperature, a warming climate has allowed the mountain pine beetle to expand its habitable range to higher elevations and its geographic distribution northward [6–8]. The mountain pine beetle colonizes 13 pine species of which it most commonly attacks lodgepole pine (*Pinus contorta* Doug.), white pine (*Pinus strobiformis* Englm.), limber pine (*Pinus flexilis* E. James), ponderosa pine (*Pinus ponderosa* Doug. ex C. Lawson) [9], and more recently, jack pine (*Pinus banksiana* Lamb.) [7,10].

Most of the mountain pine beetle's life is spent within the phloem layer of its host tree where beetles tunnel, deposit blue-staining fungi, lay eggs, and complete development from larval to adult [11–13]. The microhabitats created by beetle activity are used by numerous organisms hereby referred to as 'symbiota' or 'symbionts' [14]. The beetle's symbiota can also influence its fitness via ecological interactions (e.g., alimentary mutualists and fitness antagonists) [15,16]. For example, fungal symbiota associated with the mountain pine beetle can positively influence beetle reproductive success by providing nutritional resources [13,17,18] or negatively via parasitism and competition [19,20]. Mites (Acari) are important associates of bark beetles [16,21] and perform various direct and indirect functional roles in bark-beetle galleries such as vectoring fungi, predating each-other and beetle young, etc. [14,22,23]. They may specialize on a single beetle host, act as generalists, and form complexes of cryptic species with varying degrees of host specificity [24,25]. As in other closely related bark beetles, mite functions influence a beetle's success both directly and indirectly [26–29]. Mites can indirectly affect tree death by facilitating the transmission of plant pathogens into a new tree [29,30]. For example, mites associated with the bark beetle, *Scolytus multistriatus* (Marsham), contribute to the transmission of the destructive vascular wilt disease agent, *Ophiostoma novo-ulmi* Brasier in *Ulmus* spp.; interestingly, this fungus can significantly reduce beetle fitness [31].

Mites associate with bark beetles and other forest insects via phoresy [32], i.e., by using the beetles as a means of passive transport between trees and under the bark within a tree. Very high loads of phoretic mites can impair host mobility and negatively affect the host or 'carrier' beetle [27,33]. In carrion beetles (Coleoptera: Silphidae) mites at abnormally large densities negatively affected their host beetle fitness [34]. Several factors are known to influence phoretic mite loads across beetle populations, including forest tree composition [16,35], temperature [36–38], geographic distance and dispersal (reviewed in [39]), and beetle abundance [29,40]. Since the mountain pine beetle shows varying preferences for certain pine species [12], with the abundance of preferred host trees, and local climate varying among forests [4–6], it is expected that mite assemblages differ across beetle populations, but this is rarely tested.

Disregarding mites in the study of bark beetles not only underestimates the total biodiversity associated with bark beetles, but also reduces our ability to interpret variation in beetle population dynamics [16,21,41]. Additionally, previous studies of mountain pine beetle mites have largely been restricted to single populations, lacking comparisons of different populations [5,21,42]; thereby limiting our understanding of how specific mite assemblages relate to their beetle host success or population stage.

Here we examine, describe, and analyze the phoretic mite assemblages associated with mountain pine beetles in five different beetle populations spread across a latitudinal gradient and varying in host tree species, local climate, and beetle population stage. We use measures of alpha and beta diversity to describe these compositional differences and discuss mechanisms that might be driving these differences; specifically host tree variability, climatic variability, beetle abundance, and mite functional group. We predict that increased latitudinal separation correlates with differences in mite assemblages, and that differences may be caused by climate and environmental factors. Our aim is to set the stage for future work that would examine broad scale patterns of mite assemblages associated with mountain pine beetles and other economically and ecologically important forest insect species.

### **2. Materials and Methods 2. Materials and Methods**

### *2.1. Mite Collection 2.1. Mite Collection*

Between 2009 to 2017, 2225 mountain pine beetle specimens were collected in Arizona, Colorado, Utah, South Dakota (United States), and Alberta (Canada) during the beetle's main emergence event, which usually occurs in the late summer months between July–August (Figure 1; Table 1). Beetles were captured using a combination of baited flight traps (Lindgren funnel traps; ref. [43]) and infested logs placed inside emergence chambers [44] (Table 1). Sampling size of beetles varied due to differences in funding and resources, and the sampling strategies of each author/collaborator. Host tree species, sampling years, and beetle population phase also varied with location (Table 1). This broad dataset, therefore, represents nearly all stages and climate types of mountain pine beetle populations across its North American distribution. Differences in sample size of beetles and associated mite loads were accounted for through the use of multi-variate non-parametric statistical analyses. Between 161 and 910 beetles from each population were examined for mites. Mites were removed from each beetle, identified by clearing and mounting on glass slides under a dissecting microscope, and stored using mite collection protocols and identification resources as described in Vissa, et al. [35] for Arizona and Utah populations, Mori, et al. [42] for Alberta, Reboletti [45] for South Dakota, and Mercado, et al. [21] for Colorado. These vary only subtly in mite storage (particularly location of specimen), and otherwise converge on the same collection, clearing and mounting protocols used for processing bark-beetle mites. Between 2009 to 2017, 2225 mountain pine beetle specimens were collected in Arizona, Colorado, Utah, South Dakota (United States), and Alberta (Canada) during the beetle's main emergence event, which usually occurs in the late summer months between July–August (Figure 1; Table 1). Beetles were captured using a combination of baited flight traps (Lindgren funnel traps; ref. [43]) and infested logs placed inside emergence chambers [44] (Table 1). Sampling size of beetles varied due to differences in funding and resources, and the sampling strategies of each author/collaborator. Host tree species, sampling years, and beetle population phase also varied with location (Table 1). This broad dataset, therefore, represents nearly all stages and climate types of mountain pine beetle populations across its North American distribution. Differences in sample size of beetles and associated mite loads were accounted for through the use of multi-variate non-parametric statistical analyses. Between 161 and 910 beetles from each population were examined for mites. Mites were removed from each beetle, identified by clearing and mounting on glass slides under a dissecting microscope, and stored using mite collection protocols and identification resources as described in Vissa, et al. [35] for Arizona and Utah populations, Mori, et al. [42] for Alberta, Reboletti [45] for South Dakota, and Mercado, et al. [21] for Colorado. These vary only subtly in mite storage (particularly location of specimen), and otherwise converge on the same collection, clearing and mounting protocols used for processing bark-beetle mites.

*Forests* **2020**, *11*, x FOR PEER REVIEW 3 of 17

**Figure 1.** Map of sampling locations for the five sampled beetle populations. **Figure 1.** Map of sampling locations for the five sampled beetle populations.

*Forests* **2020**, *11*, 1102



### *2.2. Assessing Climatic Conditions 2.2. Assessing Climatic Conditions*  To assess sampling year and population location factors, we compared sampling year average

To assess sampling year and population location factors, we compared sampling year average and 30-year annual averages of precipitation and temperature extracted from the PRISM (Parameter-elevation Regressions on Independent Slopes Model) Climate Group, Oregon State University database [46] using field coordinates listed in Table 1. The same data for Alberta, Canada was retrieved from the Canadian Climate Normals database (Environment and Climate Change Canada, 2019; ref. [47]). and 30-year annual averages of precipitation and temperature extracted from the PRISM (Parameterelevation Regressions on Independent Slopes Model) Climate Group, Oregon State University database [46] using field coordinates listed in Table 1. The same data for Alberta, Canada was retrieved from the Canadian Climate Normals database (Environment and Climate Change Canada, 2019; ref. [47]). Alberta, the northern-most population, had the lowest minimum and maximum temperatures

Alberta, the northern-most population, had the lowest minimum and maximum temperatures overall; and had particularly low average minimums in the sampling year (Figure 2). South Dakota had the highest average maximum and average minimum temperatures overall. Extracted climate data showed low precipitation in both regions. Colorado's climatic conditions were similar to that of South Dakota during the 2012 sampling year but experienced a marked increase in total precipitation during the 2013 sampling year. Utah experienced similar average maximums as those seen in Colorado, however, lower average minimum temperatures were observed. Arizona, the southern-most population, had higher maximum and minimum temperatures compared to Utah, although Arizona received less total precipitation compared to Utah (Figure 2). overall; and had particularly low average minimums in the sampling year (Figure 2). South Dakota had the highest average maximum and average minimum temperatures overall. Extracted climate data showed low precipitation in both regions. Colorado's climatic conditions were similar to that of South Dakota during the 2012 sampling year but experienced a marked increase in total precipitation during the 2013 sampling year. Utah experienced similar average maximums as those seen in Colorado, however, lower average minimum temperatures were observed. Arizona, the southernmost population, had higher maximum and minimum temperatures compared to Utah, although Arizona received less total precipitation compared to Utah (Figure 2).

**Figure 2.** Basic climate information for sampling locations and years compared to 30-year annual averages. Obtained from the PRISM Climate Group database [46], Oregon and the Canadian Climate **Figure 2.** Basic climate information for sampling locations and years compared to 30-year annual averages. Obtained from the PRISM Climate Group database [46], Oregon and the Canadian Climate Normals database [47].

### Normals database [47]. *2.3. Mite Data Analysis*

*2.3. Mite Data Analysis*  Analyses were carried out in R. ver. 3.6.2 (R Core Team, 2019) using raw mite abundance data. Species accumulation curves for each population and the overall total sample pool were estimated using the *specaccum* function in the *vegan* package [48]. To test for differences in mite loads across populations, we used generalized linear models (GLM) using the *glm* function in the *MASS* package [49]. To test for pairwise differences in the assemblage of species across sampling locations based on the presence and absence of species and their abundances, a negative binomial GLM with pairwise comparisons using the *manyglm* function in the *mvabund* package [50] was used. The use of the negative binomial distribution was determined by model selection using the Akaike information Analyses were carried out in R. ver. 3.6.2 (R Core Team, 2019) using raw mite abundance data. Species accumulation curves for each population and the overall total sample pool were estimated using the *specaccum* function in the *vegan* package [48]. To test for differences in mite loads across populations, we used generalized linear models (GLM) using the *glm* function in the *MASS* package [49]. To test for pairwise differences in the assemblage of species across sampling locations based on the presence and absence of species and their abundances, a negative binomial GLM with pairwise comparisons using the *manyglm* function in the *mvabund* package [50] was used. The use of the negative binomial distribution was determined by model selection using the Akaike information criterion (AIC) [51] where Poisson, zero-inflated binomial and negative binomial distribution models were compared.

criterion (AIC) [51] where Poisson, zero-inflated binomial and negative binomial distribution models were compared. The multivariate species data were fit and analyzed with the *mvabund* package. The principle model fitting function, *manyglm* fits a GLM via resampling for each sampled species using species abundances. It allows for multiple species testing and uses a likelihood ratio test (LRT) and resampled

*p*-values to detect significance where the null-hypothesis (H0) considers beetle population to have no effect on mite species composition. To account for correlation in testing 1000 resampling iterations were used in our analysis via 'pit.trap' resampling in testing [50].

Beta diversity metrics for each population were calculated with the *betapart* package [52]. This package uses multiple-site dissimilarity measures to assess the spatial patterns of beta diversity. It also accounts for compositional heterogeneity across sites, thereby accounting for multivariate structure of dissimilarity in species composition [53]. Beta diversity calculations were conducted using both the Sorenson and Jaccard dissimilarity indices with no differences in results observed between the two. The Sorenson dissimilarity is used in all visualizations where: βsor is the overall Sorenson dissimilarity index, βsne is dissimilarity explained by nestedness, and βsim is dissimilarity explained by turnover. βsor, βsne, and βsim are measured on a 0–1 scale where βsor = 0 indicates all species are shared and βsor = 1 indicates no species are shared. βsor may also be expressed as (βsne + βsim). Degrees of nestedness were visually represented in a cluster dendrogram made within the *betapart* package. Finally, indicator species associated with each sampling location were assessed with the *indicspecies* package [54]. The indicator species analysis [55] highlights the occurrence of species that determine the community assemblage or its diversity within the given sampling area.

### **3. Results**

### *3.1. Mite Abundance, Taxa Found, and Species Richness*

A total of 4848 mites were collected across all sampled mountain pine beetle populations (Table 1). Mites represented 8 genera and 12 species (Table 2). No species were endemic to any population, while two species (*Tarsonemus endophloeus* Lindquist and *Proctolaelaps subcorticalis* Lindquist) were found in all populations. Mites from over 100 beetles were sampled from each population with species accumulation curves beginning to approach or fully approaching asymptote in all populations (Figure 3). The number of mites per beetle varied significantly (GLM analysis; *p* < 0.01) in all but three comparisons across populations (Table 3). Mite loads were highest in middle to southern populations; and decreased in northern populations. Alberta, the northern-most population, had significantly lower mite loads than all other populations. Percentage of phoretic mites per beetle was lowest in beetle populations in Alberta, Canada. Arizona and Utah had the greatest percentage of beetles carrying mites, followed by South Dakota where half the beetles carried mites (Figure 4). *Forests* **2020**, *11*, x FOR PEER REVIEW 7 of 17

**Figure 3.** Mite species accumulation curves for our five sampled mountain pine beetle populations. **Figure 3.** Mite species accumulation curves for our five sampled mountain pine beetle populations.

**Figure 4.** Percent beetles carrying mites in each sampled population.

**Table 2.** Average abundance per beetle (±standard error), functional and taxonomic information for mite species captured. Indicator species (*p* < 0.05) for each population are indicated with an asterisk (\*). Species with average abundance <0.1 are considered rare species. Significant difference between abundances listed below across populations were confirmed (see generalized linear models (GLM)

**Taxa Information Average Abundance Per Beetle (±S.E.)** 

*Tarsonemus ips* Lindq*.* Fungivore 0.34 (0.03) 3.26 (0.14) 3.82 (0.46) 1.19 (0.15) 0.79 (0.1) \*

**Dakota Utah Colorado Arizona** 

**Group Alberta South** 

analysis; Table 3).

**Mite Species Functional** 

**Table 2.** Average abundance per beetle (±standard error), functional and taxonomic information for mite species captured. Indicator species (*p* < 0.05) for each population are indicated with an asterisk (\*). Species with average abundance <0.1 are considered rare species. Significant difference between abundances listed below across populations were confirmed (see generalized linear models (GLM) analysis; Table 3).


**Table 3.** GLM test statistics for total mite abundance per beetle across five sampled populations. Asterisk (\*) shows significane at *p* ≤ 0.05.


analysis; Table 3).

**Figure 4.** Percent beetles carrying mites in each sampled population. **Figure 4.** Percent beetles carrying mites in each sampled population.

**Table 2.** Average abundance per beetle (±standard error), functional and taxonomic information for mite species captured. Indicator species (*p* < 0.05) for each population are indicated with an asterisk (\*). Species with average abundance <0.1 are considered rare species. Significant difference between abundances listed below across populations were confirmed (see generalized linear models (GLM) A pairwise comparison of the total mite species richness per beetle showed no difference between populations (*p* > 0.05; Table 4). Specific indicator species in each population varied, with some overlap occurring across Alberta, South Dakota, and Colorado. There was no overlap in indicator species between the southern-most and northern-most populations (Table 2).


**Table 4.** Pairwise population comparisons of total mite species richness per beetle.

### *3.2. Partitioning Compositional Di*ff*erences Using Measures of Beta Diversity*

The mite assemblages associated with the five bark beetle populations differed geographically based on the presence and absence of species and their relative abundances (Table 5; LRT = 7159; *p* = 0.001). *Tarsonemus endophloeus* was present in mid-latitude populations (Utah, Colorado, and South Dakota), with highest abundances in South Dakota and Colorado. *Tarsonemus endophloeus* was rare in Utah and absent from the southernmost (Arizona) and northernmost (Alberta) populations in our study (Figure 5; Table 2). *Dendrolaelaps quadrisetus* Berlese was found only in Arizona and Utah; however, these were captured in 2016 and not in 2017. *Trichouropoda utahensis* Hirschmann and Wisniewski, in particular, was unique to middle and southern mountain pine beetle populations from Utah, Arizona, and in Colorado (where *T. maeandralis* Hirschmann also was found). It is also indicative of Utah where it was most abundant (Figure 5; Table 2). Mites in the genus *Trichouropoda* are not reported in Alberta in our study location; however, they have been observed on mountain pine beetles from the southwest of the province near Banff National Park [56]. Some rare species, which include *Histiogaster arborsignis* Woodring, *Nanacarus* sp., and *Parawinterschmidtia* sp., were captured in very low numbers.

**Table 5.** (**A**) Multivariate test statistic showing likelihood ratio test (LRT) values for population effect on mite composition based on the presence/absence of species and abundances where D<sup>f</sup> = degrees of freedom; (**B**) Wald test statistics for population differences shown in (**A**); (**C**) LRT for differences in mite species across all populations. Asterisk (\*) shows significance at *p* ≤ 0.05.


captured in very low numbers.

indicative of Utah where it was most abundant (Figure 5; Table 2). Mites in the genus *Trichouropoda* are not reported in Alberta in our study location; however, they have been observed on mountain pine beetles from the southwest of the province near Banff National Park [56]. Some rare species,

**Figure 5.** Mite assemblages shown as percentage of each mite species (excluding extremely rare species) sampled across five populations. Asterisk (\*) indicates significantly different assemblages between populations (for full statistics, see Table 5). **Figure 5.** Mite assemblages shown as percentage of each mite species (excluding extremely rare species) sampled across five populations. Asterisk (\*) indicates significantly different assemblages between populations (for full statistics, see Table 5). (*T. maeandralis* Hirsch. + *T. utahensis) Proctolaelaps subcorticalis* Lindq. 65.354 <0.01 \* *Histiogaster arborsignis* Woodr. 27.042 <0.01 \*

*Histiogaster* sp. 24.011 0

**Table 5.** (A) Multivariate test statistic showing likelihood ratio test (LRT) values for population effect on mite composition based on the presence/absence of species and abundances where Df = degrees of freedom; (B) Wald test statistics for population differences shown in (A); (C) LRT for differences in mite species across all populations. Asterisk (\*) shows significance at *p* ≤ 0.05 **(A) Multivariate Test—Analysis of Deviance Table**  Df Df Difference LRT *p*-value Differences in mite composition between populations can be further explained by measures of nestedness and turnover (species replacement). On a scale of 0–1 where 1 = complete dissimilarity, the total dissimilarity in mite composition (βsor) across populations was 0.54. The majority of dissimilarity resulted from turnover (βsim) which explained 0.43 of the total dissimilarity, while nestedness (βsne) accounted for 0.114 of the total dissimilarity. *Schweibea* sp. 18.026 <0.01 \* Other/Rare spp. 9.689 0.04 Differences in mite composition between populations can be further explained by measures of nestedness and turnover (species replacement). On a scale of 0–1 where 1 = complete dissimilarity, the total dissimilarity in mite composition (βsor) across populations was 0.54. The majority of

(Intercept) 2030 4 7159 0.01 \* Population 2026 **(B) Wald Test Statistics by Population**  Population Wald Value *p*-Value Intercept (Arizona) 6.311 <0.01 \* Alberta 6.646 <0.01 \* South Dakota 4.937 <0.01 \* Colorado 3.922 <0.01 \* Utah 5.31 <0.01 \* The highest rates of turnover occurred between the southernmost beetle populations (Arizona and Utah) and the other populations, particularly Alberta (Figure 6a). High rates of turnover were also seen between southern-most population (Arizona) and each of the mid-latitude populations (Colorado and South Dakota) as well as between the northernmost population (Alberta) and the mid-latitude populations (Figure 6a). The mite assemblages of Arizona and Utah were completely nested, i.e., the mites found in Utah were a subset of those found in Arizona. Alberta mite assemblages were a subset of those found in South Dakota and other middle populations (Figure 6a,b). dissimilarity resulted from turnover (βsim) which explained 0.43 of the total dissimilarity, while nestedness (βsne) accounted for 0.114 of the total dissimilarity. The highest rates of turnover occurred between the southernmost beetle populations (Arizona and Utah) and the other populations, particularly Alberta (Figure 6a). High rates of turnover were also seen between southern-most population (Arizona) and each of the mid-latitude populations (Colorado and South Dakota) as well as between the northernmost population (Alberta) and the midlatitude populations (Figure 6a). The mite assemblages of Arizona and Utah were completely nested, i.e., the mites found in Utah were a subset of those found in Arizona. Alberta mite assemblages were a subset of those found in South Dakota and other middle populations (Figure 6a,b).

**Figure 6.** (**a**) Sorenson dissimilarity indices of Beta (β) diversity for mite communities partitioned into nestedness (βsne) and turnover (βsim) from southern-most (Arizona [AZ]) to northern-most (Alberta, Canada [CAN]) sampled populations; (**b**) cluster dendrogram showing degree of nestedness based on overall dissimilarity (βsor). **Figure 6.** (**a**) Sorenson dissimilarity indices of Beta (β) diversity for mite communities partitioned into nestedness (βsne) and turnover (βsim) from southern-most (Arizona [AZ]) to northern-most (Alberta, Canada [CAN]) sampled populations; (**b**) cluster dendrogram showing degree of nestedness based on overall dissimilarity (βsor).

evident between Alberta and all other populations. Turnover and nestedness patterns varied

although species richness did not show any significant variability across populations.

Geographically distant mountain pine beetle populations have different symbiotic mite

**4. Discussion** 

### **4. Discussion**

Geographically distant mountain pine beetle populations have different symbiotic mite assemblages despite some commonalities reflected in species overlap. Sorenson's dissimilarity values for beta diversity indicate that the mite assemblage of beetle populations geographically closer to one another are more similar than populations that are geographically more distant. This was most evident between Alberta and all other populations. Turnover and nestedness patterns varied although species richness did not show any significant variability across populations.

### *4.1. Di*ff*erences in Species Richness and Abundance*

No significant differences were detected in mite species richness, presumably because richness per beetle was generally very low for all populations. Species accumulation curves approached asymptote in all populations (although not at the same rates), indicating that our sampling efforts adequately captured the overall mite richness associated with these populations. However, the disparities in species accumulation curves between populations may also be due to the possible presence of cryptic species that could not be uncovered. The number of mites per beetle (mite load) did not vary on a latitudinal gradient (as predicted), suggesting that other factors may better explain these differences. We observed that post-outbreak and 'endemic' stage populations had higher mite abundances, and observably higher total richness (Table 2) suggesting that mite loads may be indicative of beetle population phase. Hofstetter, et al. [57] show that higher mite loads dominated by *Tarsonemus* mites correlated with decreased beetle progeny densities because of mite introduction of antagonistic blue-stain fungi. While our sampling efforts could not test particular hypotheses, our observed mite load trends may provide insight on population specific ecological dynamics in mountain pine beetles. High mite loads can affect beetle flight [27,33] or influence their reproductive success within trees [58]. The specific effects of different mite levels on mountain pine beetles has not been tested. Our analysis of mite community patterns provides a foundational step towards understanding these effects.

The finding that *Tarsonemus ips* and *Proctolaelaps subcorticalis* occurred in all populations suggests that these mites may not be affected by latitude, pine tree host, climate, or beetle population stage. *Trichourpoda utahensis* was only found in Arizona and Utah; and although also found in Colorado, occurred alongside a second *Trichouropoda* species, *T. maeandralis* (which was not the case for Arizona and Utah). *Trichouropoda* mites carry fungi that differ from the mutualistic fungi carried by the mountain pine beetles, and its relationship with collapsing beetle population has been suggested [16].

### *4.2. Beta Diversity Patterns*

As expected, populations relatively closer together shared more species between them than those geographically distant. High nestedness at these sites may indicate a low overall phoretic mite biodiversity [59]. Utah and Arizona (and likewise South Dakota and Alberta) showed nested patterns where mite species were likely removed rather than replaced. This could be due to the large latitudinal and subsequent climate shift between these populations resulting in loss of species with changes in local climate.

High rates of species turnover are explained by environmental differences causing some species to be replaced by others more suited to the new environment [60]. Although not specifically tested in this study, we hypothesize that environmental filtering, specifically temperature, serves as a likely explanation for differences in mite community assemblages across latitude. Beta diversity patterns driven by turnover are associated with temperature seasonality, which may be experienced across a latitudinal gradient [61]. Other complementary abiotic factors may contribute to high levels of turnover, particularly dispersal ability. Hill, et al. [62], showed that spatial patterns and environmental heterogeneity play a large role in explaining high levels of dissimilarity [in community assemblage] among sites in passively transported invertebrates and are better explained by increased turnover. However, dispersal limitations are often defined for arthropods by temperature [63], which may be the ultimate driver of the latitudinal patterns.

Differences in host tree species may also help explain mite composition patterns across bark beetle populations. Our data analysis could not accommodate testing for differences in mite composition by host tree because host tree and location are confounding factors except for South Dakota and Colorado which shared the same host tree species (Table 1), and beetles collected via flight traps may come from a variety of tree species. Fungal symbiota associated with the mountain pine beetle induce different secondary chemical responses in different host tree species which may contribute to the tree's attractiveness as a host environment for the beetle and its associates [64,65]. In our study, although Colorado and South Dakota populations shared the same host (*Pinus ponderosa*), the overall mite compositions associated with these populations differed significantly. However, it should be noted that sampling years varied. However, these two populations did differ in their beetle population stage—Colorado was a collapsing population from an epidemic while the South Dakota population was increasing to epidemic status. Vissa et al. [35] reported vastly different mite loads associated with beetles attacking two different tree species within a single region in Portugal, suggesting that dominant host tree type can also affect beetle-associated mite abundances. Here we found that the population from limber pine (Utah) had beetles with the greatest mite loads. Tree genotype has also been shown to drive different defense production, insect preference, and long-term population dynamics [66,67].

Temperature affects patterns in species diversity across latitudinal gradients [68]. Our colder population of Alberta had the lowest mite diversity while the others, similarly warmer populations, had contrasting higher mite diversities. Evans, et al. [69] showed that a *Trichouropoda* sp. associated with the southern pine beetle had the highest reproductive success in warmer temperatures. This was reflected in our results where populations with higher average minimum and maximum temperatures also had more mites per beetle. Species of *Trichouropoda* were also absent in Alberta but prevalent in Colorado, Utah, and Arizona, where both temperature and precipitation were higher during development years; however, this may be confounded with the population stages in these populations.

### *4.3. Implications for Shifts in Functional Groups of Mites*

The concept of environmental filtering in symbiotic mite communities associated with bark beetle systems is not well explored. Swenson et al. [70] hypothesized that abiotic (i.e., environmental) filtering not only alters species assemblages in plants, but also alters the functional diversity. Thus, in climatically 'extreme' regions, not only are there fewer species [70], but the functional groups can be altered as well [71]. While we were unable to identify specific shifts in the functional groups of mites, the observed compositional differences may contain vital information on the resulting functional shifts.

High taxonomic turnover does not directly translate to a significant functional turnover [72]. Fish faunas with high taxonomic beta diversity were associated with low functional beta diversity, suggesting that although species assemblages varied, they maintained shared functional attributes [72]. Our methods were unable to tease apart the functional diversity patterns associated with symbiotic mite populations or allow us to specifically identify which mite species are being replaced and by whom; however, we did observe some uniformity in the abundant functional groups across all populations. Fungivores (particularly *Tarsonemus* mites) were abundant in each epidemic population, followed by a generalist predator (*P. subcorticalis*), suggesting that these species likely affect beetle population dynamics (either by promoting the spread of fungi or by direct predation).

### **5. Conclusions**

To our knowledge, this is the first study to examine and partition the differences in the beta diversity of symbiotic mite communities associated with an economically important North American bark beetle (the mountain pine beetle) across a wide range of its distribution; and may serve as a model system for estimating diversity patterns of symbiota associated with other major insect pests in pine forest systems. As anticipated, the composition of mite biodiversity associated with mountain pine beetles varied geographically although no mite assemblage was entirely unique, and environmental differences driven by geography can best explain these differences. Our assessment of mite species composition patterns provides foundational information for further steps to identifying key species interactions shaping forest insect systems.

**Author Contributions:** Corresponding author S.V., R.W.H. and D.A.U. contributed significantly to the conceptualization, mite collection, sample processing and editorial process of this manuscript. D.A.U. additionally provided significant assistance in statistical analyses. Authors D.M., J.E.M., and B.A.M. conducted extensive sampling and mite data processing for South Dakota, Colorado, and Alberta respectively; and also actively contributed to the written materials in this manuscript. Author W.K. served as the primary identification resource for mite data in this manuscript and also contributed as a manuscript editor. Finally, author M.L.E. conceptualized experimental design of the Canadian component, provided some editorial help, and facilitated the Canadian research work by obtaining the necessary funding. All authors have reviewed this manuscript in its entirety and approve of its submission. All authors have read and agreed to the published version of the manuscript.

**Funding:** The Canadian component of this manuscript was largely funded by the National Science and Engineering Research Council (NSERC) Discovery Grant.

**Acknowledgments:** We would additionally like to acknowledge the contributions of the following organizations, researchers, technical support individuals: The United States Department of Agriculture (USDA)/Forest Service (FS) and Rocky Mountain Research Stations (RMRS), Northern Arizona University (NAU), Utah State University, Black Hills National Forest, Alberta Agriculture and Forestry, John Moser, Elizabeth Alden, David Soderberg, Devin Letourneau, Pam Melnick, Jeffrey Mallete, Madison Martz, Karis Miller and Roy St. Laurent (NAU Statistical Consulting Lab). We additionally thank Heather Proctor, University of Alberta for her assistance with mite identifications and for editing.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
