Weak spatial-genetic structure in a native invasive, the southern pine beetle (Dendroctonus frontalis), across the eastern United States

General approaches to assess lack of structure

Sex-biased dispersal

Different degrees of site fidelity between males and females can affect detection of population differentiation when using biparentally inherited microsatellite markers. Impacts are most pronounced when sex-specific dispersal rates and distances are highly asymmetric, and when sampling includes many individuals that have dispersed but not yet reproduced (Prugnolle & De Meeûs, 2002). To evaluate this, we analyzed male-only (n = 185 individuals from 9 sites) and female-only (n = 70, 8 sites) partitions, in addition to combined datasets (n = 255, 9 sites; Fig. 1).

Development and validation of new microsatellite loci

In the present study, we extended upon the assessment of a suite of microsatellite loci identified by Havill et al. (2019). Specifically, nine additional loci were selected from among those identified using the software QDD (Meglecz et al., 2014) following low coverage whole genome sequencing of one male D. frontalis from Homochitto National Forest, Mississippi (NCBI BioProject Accession Number PRJNA493650). To assess the suitability of these new loci, multiplex polymerase chain reaction (PCR) amplifications were performed for sets of three loci, each with a different 5′ fluorescent label on the forward primer and mostly non-overlapping allele size ranges (Table 1). Reactions were conducted in 15 µL volumes that contained 7.5 µL Type-it^® Microsatellite PCR Kit master mix (Qiagen, Valencia CA), 1.5 µL dH₂O, 0.75 µL of each primer (10 mM), and 1.5 µL genomic DNA. The following “touchdown” thermal cycling conditions were used: 95 °C for 2 min (1 cycle), 95 °C for 30 s, 61 °C minus 2 °C per successive cycle for 30 s, 72 °C for 45 s (5 cycles), 95 °C for 30 s, 51 °C for 30 s, 72 °C for 45 s (30 cycles), and 60 °C for 30 min (1 cycle).

Amplified fragments were run on an ABI 3730 sequencer (Life Technologies) with a Liz 500 ladder (Gel Company, San Francisco CA) at Yale University’s DNA Analysis Facility on Science Hill. Genotype scoring was performed using GENEIOUS v10.0.5 (Kearse et al., 2012). Following Havill et al. (2019), we assessed Mendelian inheritance patterns for each of the new loci in two local populations (i.e., 27 individuals from Sicily Island, Louisiana, and 28 individuals from Tombigbee, Mississippi). For each locus and population, evidence for departures from Hardy-Weinberg Equilibrium (HWE) were examined via exact tests (Guo & Thompson, 1992), with significance assessed using 10,000 Markov Chain permutations. Observed (H_O) and expected (H_E) heterozygosity were also calculated. Independent segregation of alleles among loci was assessed using all D. frontalis samples (255 individuals) by testing for departures from Linkage Equilibrium (LE) for all possible pairs of 33 loci (i.e., nine reported here, plus 24 from Havill et al., 2019), again using 10,000 Markov Chain permutations. Given that nine different sampling sites were included, Fisher’s combined probability test was used to calculate P-values across populations for each locus pair. These analyses were performed in GENEPOP v4.5.1 (Rousset, 2008), with sequential Bonferroni correction (Holm, 1979) of P-values associated with LE tests.

Identification and omission of loci that may elevate stochastic noise

MICRO-CHECKER v2.2.3 (Van Oosterhout et al., 2004) was used to determine whether there was evidence for null alleles, and if so, to estimate their frequencies. For each locus in each local population, a null allele was inferred to exist if there was a significant excess of homozygotes that were evenly distributed across all homozygote classes. This was determined by Fisher’s combined probability test, with P-values determined using 1,000 permutations. Frequency of the null allele (r) was then was estimated using the method of Chakraborty et al. (1992), which assumes that missing data are not attributable to homozygous nulls, and that populations are panmictic. To identify inbred populations for removal, we relied on the tendency for inbreeding to generate homozygote excess at many loci (cf. locus-specific impacts of null alleles). Here we defined outbred populations as those in which fewer than seven loci (i.e., <20%) showed significantly (P < 0.05) positive F_IS values. For each locus and population, the null hypothesis of panmixia was compared to an alternative hypothesis of heterozygote deficiency using a U-test (Rousset & Raymond, 1995), performed in GENEPOP with 10,000 Markov Chain permutations. After omitting outbred populations, we rank-ordered the loci based on the level of noise contributed by putative null alleles. To ensure that both prevalence and severity of the null allele at a locus were considered, the r value from each population with a significant excess of homozygotes was summed (r_cumulative).

The number of contiguous repeat units may be a useful proxy for microsatellite mutation rate, given that DNA structure affects the opportunity for slipped-strand mispairing (Bhargava & Fuentes 2010, and references therein). Here we used the reported positive association between repeat number and mutation rate to identify which loci might be predisposed to having alleles identical in state but not by descent. Briefly, based on one sequence per locus generated during marker development (Schrey et al., 2007; Havill et al., 2019; this study), we first estimated the PCR amplicon length for each locus given the locations of forward and reverse primer annealing sites. Next, the observed number of contiguous repeats was scaled to reflect what should be present within the median-sized allele (i.e., a reference allele identified from the eastern United States gene pool as whole). Rank-ordering of loci based on the potential level of stochastic noise contributed by homoplasy followed the inferred number of repeats of each reference allele.

The omission of loci that may elevate overall noise was based on joint consideration of null alleles and homoplasy. Along each of these two axes, we iteratively omitted the lowest ranked locus one at a time, until a two-sample t-test (one-tailed) determined that a significant (P < 0.05) reduction in the mean value of the original metric used for rank-ordering (compared to that of the 33-locus dataset) had been achieved. To assess whether the low noise dataset was comprised of a subset of loci with inadvertently low information content, we calculated Paetkau et al.’s (1995) Probability of Identity (PI; the chance that two individuals randomly drawn from a panmictic population have the same genotype) using GENALEX v6.5 (Peakall & Smouse, 2006; Peakall & Smouse, 2012). We then compared mean PI between the reduced vs. full dataset, using a t-test as above.

Identification and retention of loci that may enhance overall signal

For species in which population divergence is shallow, “outlier” microsatellite loci (e.g., those displaying unusually strong differentiation, potentially due to effects of divergent selection on linked genomic regions) may be particularly useful for identification of population structure (Russello et al., 2012). This is because recent divergence (e.g., over post-glacial timescales or shorter) may be slow to register in most loci, particularly when it is coupled with large effective populations sizes. To explore the potential benefit of focusing on loci with the strongest signatures of spatial-genetic structure in D. frontalis, we calculated global F_ST (Weir & Cockerham, 1984) using the augmented dataset, for each of the 33 loci separately, and then plotted the distribution of these values in rank order. Natural breaks were identified qualitatively, based on the magnitude of change (i.e., ΔF_ST = F_ST locus n −F_{ST locus n+1}, where n is the rank order position). These breaks defined the thresholds for retention vs. exclusion of loci from the high signal dataset.

Genotypic clustering analyses

To infer the number of natural genetic groups (K) and their members, we analyzed microsatellite datasets using three approaches: STRUCTURE v2.3.4 (Pritchard, Stephens & Donnelly, 2000), BAPS v6.0 (Corander, Waldmann & Sillanpää, 2003), and Discriminant Analysis of Principal Components (DAPC; Jombart, Devillard & Balloux, 2010). The first two are Bayesian approaches and assume HWE and LE within “true” clusters. However, STRUCTURE is strictly individual-based and does not make use of geo-referenced samples (although priors on population of origin can be applied), but it has been shown to tolerate clinal variation relatively well (Chen et al., 2007). BAPS can be run using either individuals or groups of individuals as the basic units of analysis, and it explicitly considers spatial coordinates, which can help when differentiation is weak. DAPC can also be used to identify natural groups of individuals a posteriori, but this method does not attribute the underlying cause of such groups to population genetic processes, nor does it incorporate spatial information. In all clustering analyses, we assessed K values from one up to and including the total number of sampling sites (i.e., maximum K = 8 or 9, depending on the dataset under consideration).

STRUCTURE. All runs employed the correlated allele frequency and admixture ancestry models. Estimated -ln likelihood scores were obtained for each value of K, with 20 replicates each. A burn-in of 50,000 Markov chain Monte Carlo (MCMC) generations and run length of 200,000 MCMC generations were used, with other parameters set as default. The best-fit value of K was chosen using Evanno, Regnaut & Goudet’s (2005) ΔK method implemented in STRUCTURE HARVESTER v0.6.94 (Earl & vonHoldt, 2012). However, the ΔK method is unable to assess support for K = 1, and it may be subject to bias with respect to identifying only the top (e.g., K = 2) level of hierarchical structure (see Janes et al., 2017). Accordingly, we also used a visual assessment of Ln Pr(X—K) plots, from which the smallest K that captured the major structure in the data was accepted, following Pritchard & Wen (2003).

BAPS. To enable spatial clustering of individuals, we created non-redundant GPS coordinates via the addition of fine-scale jitter (i.e., a random number at the 4th and 5th decimal places, equivalent to ≤ 11.1 m at the equator). For spatial clustering of groups of individuals, each of the 8 or 9 sampling locations (depending on the dataset) defined an a priori group, and the original spatial coordinates were used. Accordingly, the complexity of Voronoi tessellations differed considerably between the two approaches. Spatial data were analyzed with associated microsatellite genotype data, with 10 replicates per run. The best-fit K was identified via log marginal likelihood scores (Cheng et al., 2013).

DAPC. Individual-based multilocus microsatellite genotypes were first transformed via principal component (PC) analysis to identify uncorrelated variables to be used for subsequent discriminant analysis. The best-fit number of clusters was determined via k-means clustering, with the lowest Bayesian Information Criterion (BIC) score used to choose K. Following Miller, Cullingham & Peery (2020), an initial preliminary DAPC was conducted to determine the optimal number of PCs to retain, followed by a second final DAPC using the chosen number of PCs. In the initial assessment, the optimal number of PCs was determined using the cross-validation procedure, considering 20 to 300 PCs (in increments of 20) with 30 replicates at each level of PC retention, and 90% of the data comprising the training set. Pre-defined groups based on sampling location information were used only for the purpose of quantifying the number of de novo clusters that were inferred to be present at a given sampling location. DAPC was implemented using ADEGENET v2.1.2 (Jombart, 2008) in R v3.6.1 (R Core Team, 2020).

Isolation-by-distance analyses

To assess evidence for subtle gradients in genetic variation across the landscape, we conducted IBD analyses of our microsatellite datasets using either sampling sites or individuals as the basic units. The rationale for these two approaches was that population allele frequencies typically can change rather gradually over time and therefore provide insights into historical processes, whereas diploid genotypes are much more labile owing to reshuffling of alleles every generation in sexually reproducing species, and therefore have the potential to closely track recent or on-going processes (Sunnucks, 2000; Garrick, Caccone & Sunnucks, 2010; Epps & Keyghobadi, 2015).

Population-based metrics. Three measures of genetic distance were used: Cavalli-Sforza & Edwards’s (1967) chord distance (D_c), which has desirable properties for detecting IBD in the presence of null alleles (Séré et al., 2017); Nei’s (1972) standard genetic distance (D_s), which increases linearly with geographic distance if mutation follows an infinite allele model and occurs at a constant rate; and Weir & Cockerham’s (1984) estimate of F_ST linearized as F_ST/(1- F_ST) following Slatkin (1995), to account for the high mutation rate of microsatellite loci.

Individual-based metrics. Three metrics for genetic similarity were used: the kinship coefficient (D_kf, Cavalli-Sforza & Bodmer, 1971); proportion of shared alleles (D_ps, Bowcock et al., 1994), or relatedness (rˆ, Lynch & Ritland, 1999). Briefly, for a given locus, D_kf is the probability that a randomly chosen allele from one individual is identical-by-descent to that of second individual (summed over all loci and alleles). Conversely, when comparing two mutlilocus genotypes, D_ps is calculated as number of shared alleles summed over loci / (2 × number of loci compared). Unlike the other two metrics, Lynch & Ritland’s (1999) methods-of-moments estimator of relatedness makes use of information on population allele frequencies (e.g., a shared rare allele is considered more likely to be identical-by-descent than a shared common allele). Given that rˆ is asymmetrical, following Ritland (2000), we used the average of both directions.

MICROSATELLITE ANALYZER v4.05 (Dieringer & Schlötterer, 2003) was used to calculate D_c, D_kf, and D_ps, whereas GENALEX was used to calculate D_s, F_ST and rˆ, as well as linear pairwise geographic distance matrices based on latitude and longitude coordinates, using a modified Haversine Formula. For all datasets and distance metrics, the significance of IBD was assessed via Mantel (1967) tests with 999 permutations.

Development and validation of new microsatellite loci

For the nine new loci developed here, amplification success across all eastern United States D. frontalis samples was reasonably high (mean of 4.2% missing data per locus). Three loci showed significant departures from HWE in the direction of homozygote excess. However, in all cases this was limited to just one of the two local collection sites (Table 1), suggesting that there were no intrinsic issues with allelic inheritance (i.e., the new microsatellite loci appear to be autosomal, diploid, and single copy). Levels of with-population polymorphism were modest (mean of ∼4.3 alleles per new locus for the two exemplar local populations). Across all eastern United States D. frontalis samples, there was a mean of 8.1 alleles per new locus (Table 1). Based on our assessment of LE for all 33 loci (528 locus pairs), 22 pairs were significant at the P <0.05-level, but none remained so after sequential Bonferroni correction to account for multiple tests. Accordingly, there was no strong evidence for non-independence among loci in the augmented dataset. Preliminary analyses using all loci indicated weak genetic structure among sampling sites (F_ST = 0.009, P = 0.001).

Identification and omission of loci that may elevate overall noise

Null alleles. Two local populations had seven or more loci with significantly positive F_IS values: Homochitto, MS, and Woolford, MD. These were considered inbred and therefore omitted from consideration when assessing null alleles. MICRO-CHECKER identified 13 loci that may have a null allele (Table 2). However, prevalence was low (i.e., restricted to a single population) for 10 of these. Of the remaining loci, two had null alleles in four populations, and one had nulls in two populations. Across the 20 cases of putative null alleles (out of 33 loci ×7 retained populations = 231 potential cases), mean estimated r was 0.245 (range: 0.139–0.431). Rank-ordering followed by iterative removal revealed that omission of the five worst loci was required to significantly reduce the mean r_cumulative value (t = 1.720, d.f. = 40, P = 0.047). The removal of these loci did not cause a concomitant reduction in mean information content, as measured by PI (t =  − 0.435, d.f. = 57, P = 0.332).

Homoplasy. For 31 loci, the available allele sequence trimmed at the 5′ ends of both primers yielded a predicted PCR product that was within the empirically determined size range for eastern United States D. frontalis. For two loci with unexpectedly short or long amplicons, an alternative priming site (i.e., a region of high sequence similarity located up or downstream of the original target site) was readily identifiable (Table 3). Three loci constrained two distinct microsatellite regions. These were either composed of different (i.e., compound) repeat motifs located immediately adjacent to one other (locus Dfr-16; Schrey et al., 2008), or the same (i.e., impure) repeat motif interrupted by a 2-bp (locus Dfr-24; Schrey et al., 2008) or 3-bp (locus SPB1534; this study) insertion. In all three cases, only the longest repeat region was used as a proxy for potential for homoplasy. Following re-scaling (see Methods), across 33 loci the average value of the median contiguous number of repeats was 7.6 (range: 4.945–14.901; Table 3). The five most concerning loci needed to be removed in order to significantly reduce the number of repeat units (t = 1.819, d.f. = 50, P = 0.037), but their omission did not also reduce mean PI (t =  − 0.780, d.f. = 58, P = 0.219).

Of the two sets of five omitted loci, one (Dfr-24) was common to both approaches used to reduce noise. Accordingly, the largest subset of loci with reduced null allele and homoplasy issues was made up of 24 loci. Reassessment of the impacts of our strategic removal of nine loci to create the low noise dataset showed that the reduction in mean r_cumulative remained close to significant (t = 1.520, d.f. = 42, P = 0.068), and so did reduction in mean number of repeats (t = 1.499, d.f. = 53, P = 0.070). As before, there was no inadvertent change in mean PI (t =  − 1.180, d.f. = 51, P = 0.122). Although our low noise dataset happened to contain the same number of loci originally analyzed by Havill et al. (2019), their compositions differed by eight loci.

Identification and retention of loci that may enhance overall signal

Based on a visual examination of the ΔF_ST plot, several inflection points marking pronounced decreases in F_ST among adjacent rank-ordered loci were apparent (Fig. 2). Aside from steep declines that would have set a threshold for retention of too few loci for meaningful population genetic analyses (i.e., ≤ 5 loci), there were two other inflection points: one at the transition from the 10th to 11th ranked locus, and another at the transition from 19th to 20th. Accordingly, we created two alternative versions of the high sensitivity dataset (HSDS). One contained 10 loci with the highest F_ST values, and the other was an expansion of this, where we retained the 19 highest ranked loci (Table 4). The 10-locus and 19-locus high signal datasets contained 3 and 6 loci, respectively, that were not shared with the 24-locus low noise dataset. Compared to Havill et al.’s (2019) original dataset, the 10-locus and 19-locus high signal datasets contained 2 and 5 novel loci, respectively. Ultimately, all datasets analyzed in the present paper differed from each other and from the previous study, thereby enabling meaningful comparisons.

Genotypic clustering analyses

STRUCTURE. Across all 12 datasets (i.e., ADS, LNDS, and HSDS with 10 or 19 loci; males and females combined, or separated) there were five instances where the two alternative approaches for identifying the best-fit value of K were in conflict. In these cases, Pritchard & Wen’s (2003) qualitative method inferred K = 1 whereas Evanno, Regnaut & Goudet’s (2005) ΔK method inferred K = 2, 3, or 4, depending in the dataset (Table 5; but note that ΔK is not capable of assessing K = 1). More generally, irrespective of whether there was agreement or conflict between methods, wherever K>1 was inferred, the optimal clustering solution almost always contained “ghost clusters” (sensu Guillot et al., 2005); i.e., those for which no individuals were strongly assigned with a membership coefficient of Q>0.50). In some cases, all clusters were nonsensical. For instance, based on ΔK, the best-fit value was K = 3 for LNDS males plus females, and LNDS females-only, yet not a single individual was strongly assigned; Table 5). Additionally, when K>1 there was also an absence of geographic localization of clusters. Ultimately, as with Havill et al.’s (2019) STRUCTURE analysis of a 24-locus D. frontalis dataset, our reanalysis using augmented and/or sub-setted datasets did not provide any indications of discrete population structure.

BAPS. The optimal partition of geo-referenced molecular data was K = 1 for all datasets, using both individual- and population-based clustering approaches (Table 5). Given this outcome, we performed post-hoc HWE exact tests in GENEPOP, treating all individuals as members of a single population. Each of the 12 datasets showed highly significant departures from HWE (all P < 0.0005), indicating that eastern United States D. frontalis cannot simply be characterized as panmictic.

DAPC. Across all datasets, there were only two instances where K = 1 was inferred (Table 5; but note that DAPC was not intended for this purpose; see Miller, Cullingham & Peery, 2020). However, wherever K ≥2, there were strong indications that the groups were artefacts rather true reflections of population structure. For example, there were no cases where a given local sampling site contained members of only a single inferred cluster. Indeed, in seven of the datasets with best-fit K ≥2, all sampling sites contained representatives of all inferred clusters (Table 5); such levels of coexistence among putatively distinct gene pools has no plausible biological explanation.

Isolation-by-distance analyses

Population-based metrics. Regardless of which microsatellite dataset was analyzed, or which genetic distance metric was used, no significant IBD was detected (all P ≥0.130; Table 6).

Individual-based metrics. All instances of significant IBD were limited to Lynch & Ritland’s (1999) relatedness (note that such relationships are negative). Furthermore, this significant IBD was almost always associated with datasets that included males, with the one exception being the female-only 19-locus high sensitivity dataset (Table 6). Despite relatively strong significance (all P = 0.001, except for the aforementioned female-only dataset, which had P = 0.031; Table 6), the nature of relationships between relatedness and geographic distance were consistently weak [augmented dataset with males plus females: slope = −2 ×10⁻⁶, correlation coefficient (r) = 0.036; augmented dataset with males only: slope = −2 ×10⁻⁶, r = 0.037; low noise dataset with males plus females: slope = −2 ×10⁻⁶, r = 0.020; low noise dataset with males only: slope = −2 ×10⁻⁶, r = 0.030; and even weaker for the 19-locus high signal dataset with females only: slope = −3 ×10⁻⁶, r = 0.001]. Overall, these outcomes are consistent with female-biased dispersal; however, males do not appear to be strongly philopatric.

Clusters

Regardless of whether we increased the number of microsatellite loci (augmented dataset) or analyzed subsets of loci with more favorable signal-to-noise ratio (low noise and high signal datasets), no spatially abrupt genetic subdivisions were detected. Using simulations, Miller, Cullingham & Peery (2020) showed that in the presence of even quite low levels of on-going gene flow (i.e., migration rate, m ≥ 0.005), the best-fit number of clusters inferred using STRUCTURE may be systematically underestimated, leading to erroneous inferences of K = 1. Following the recommendations of those authors and others, we compared outcomes across several analytical methods with different assumptions and found that the inferred absence of discrete clusters was robust (Table 5). We interpret this to indicate that either D. frontalis does not exhibit discrete population structure, or gene flow is sufficiently high that such structure is not detectable with our data. Other types of molecular markers, such as single nucleotide polymorphisms (SNPs) assayed via sub-genomic sampling methods that are applicable to non-model species (e.g., Elshire et al., 2011) could reveal subtle fine-scale spatial-genetic structure. For example, compared to inferences from microsatellites, SNPs have identified additional structure in the mountain pine beetle, D. ponderosae (Batista et al., 2016; Janes et al., 2018). That said, although only 4 × to 12 × more bi-allelic SNPs than multi-allelic microsatellite loci might offer comparable resolution of population structure (e.g., Liu et al., 2005; but see Haasl & Payseur, 2011 for a potentially much higher ratio), microsatellites have desirable properties for assignment tests, kinship analyses and estimating heterozygosity, and therefore remain valuable (Narum et al., 2008).

Clines

When using sampling sites as the unit of analysis, we found that microsatellite dataset composition did not impact conclusions about the absence of gradients of genetic differentiation. This outcome was also robust to the chosen genetic distance metric (Table 6). Range expansion is often considered in the context of a series of founder events that each give rise to a new population with reduced genetic variation along the moving wave of advance, leading to genetic differentiation from the source population. However, when range expansion is recent or on-going, multi-directional, and/or effective population sizes are consistently large, there can be a lag time or weakening of bottleneck effects (also see Roques et al., 2012 for an example of how Allee effects at a low-density wave front can prevent successive loss of genetic diversity). Likewise, if occasional long-distance dispersal is involved, new populations can establish in leaps and bounds, far from the parent population (Nichols & Hewitt, 1994; Hewitt, 1996; Hewitt, 1999; Ibrahim, Nichols & Hewitt, 1996). Under any of these scenarios, IBD is not expected to be immediately detectable. Indeed, the current expansion of D. frontalis to the northeastern United States is both recent and on-going: the species spread throughout New Jersey by the mid- to late-2000’s, to Long Island in New York by 2014, and was trapped in Connecticut, Rhode Island and Massachusetts by 2015–2016 (Dodds et al., 2018). Effective population sizes certainly have the potential to be very large, but little is known about non-outbreak populations, making it difficult to speculate about whether genetic drift would be a strong driver of allele frequency differences over increasing geographic distances. However, the capacity for long-distance dispersal by D. frontalis has been well-documented: this could occur via active flight over several kilometers, or via passive above-canopy wind-assisted dispersal over tens of kilometers (Jones et al. 2019, and references therein). Ultimately, lack of IBD is not entirely unexpected, at least not in the newly invaded portion of this species’ current geographic range. However, other explanations exist. For example, detection of IBD might be scale-dependent, whereby its signal fades at increasingly larger geographic distances owing to an upper bound on the maximum attainable genetic distance (Castric & Bernatchez 2003). Such non-linearity might be revealed by iterative reanalysis using different distance thresholds (Van Strien, Holderegger & Van Heck, 2015), provided that straight-line geographic distance (cf. isolation-by-environment; Wang & Summers, 2010) is the primary driver of genetic divergence.

Female-biased dispersal

Despite overall weak spatial-genetic structure among eastern United States D. fontalis, we did find significant individual-based IBD attributable to males (Table 6). Although this inference of female-biased dispersal was limited to Lynch & Ritland’s (1999) measure of relatedness, it was robust to different compositions of microsatellite loci across datasets. Indeed, this finding is concordant with life history traits of D. frontalis. Females are the pioneering sex, who locate a suitable host before recruiting males via release of aggregation pheromones (Jones et al. 2019). Furthermore, females are larger than males, and given that body mass is correlated with stored energy reserves, their greater capacity for flight (mean distance of 3.4 km vs. 2.7 km for males) in a flight mill experiment was attributed to sexual size dimorphism (Kinn et al., 1994). Notably, the notion that males are themselves not strongly philopatric is supported by capture-mark-recapture data, which showed no significant differences in dispersal between the sexes over distances up to 1 km (Turchin & Thoeny, 1993).

There are several reasons that could explain why individual-based IBD was detectable in our study, yet population-based IBD was not. For instance, this may simply reflect an issue of sample size (i.e., 2,415–32,385 pairwise comparisons among individuals vs. only 28–36 among sampling sites depending on the dataset; but note the elevated pseudo-replication associated with the former). Alternatively, a difference in the timescales over which spatial-genetic structure has evolved and/or is assessed may be responsible. Indeed, individual-based genetic distances may better reflect contemporary processes such as intra-generational dispersal, whereas population-based metrics primarily measure historical connectivity given that most signal comes from the accumulated multi-generational effects of dispersal and gene flow (Sunnucks, 2000; Garrick, Caccone & Sunnucks, 2010; Epps & Keyghobadi, 2015). In the context of this study, however, these explanations are plausible only if the other two individual-based genetic distance metrics that we assessed (i.e., D_kf and D_ps) had low power. One potential advantage of Lynch & Ritland’s (1999) r is that loci with more alleles—particularly those with rare alleles—provide more information about relatedness, such that highly polymorphic genetic markers such as microsatellites are very powerful (Ritland, 2000). That said, this remains speculative, as simulations have shown D_ps to be more sensitive than r under a variety of IBD scenarios, albeit under simplified conditions (Shirk, Landguth & Cushman, 2017).

Reconciliation with previous work

Our findings for D. frontalis are inconsistent with those of Schrey et al. (2011) who reported both clusters and clines, in the form of an east–west division approximately coincident with southern Appalachian Mountains and population-based IBD across the eastern United States. Those authors used eight microsatellite loci, all of which were included in Havill et al.’s (2019) 24-locus dataset, as well as the present study’s 33-locus augmented dataset. Generally, trade-offs between the number of loci and number of individuals favor adding loci (Landguth et al., 2012), but given the shared loci, it seems unlikely that the marker set alone can account for contrasting inferences. However, geographic and temporal sampling did differ. Whereas Schrey et al.’s (2011) study design included 19 sites with an average of 63 individuals per site (range: 26–100) collected between fall 2004 to spring 2005, ours included 9 sites with an average of 28 individuals (range: 21–31) collected >10 years later (i.e., summer 2016 to summer 2017, except for Ponte Vedra, Florida, collected in summer 2013). Thus, although the overall spatial scale of sampling was similar between studies (see Introduction), effort allocation was not. Indeed, Landguth & Schwartz (2014) have shown that the optimal allocation of individuals for accurate estimation of F_ST changes depending on the extent of gene flow limitation, with misallocations leading to either over- or under-estimation of F_ST, depending on the circumstances. Furthermore, direct comparison between studies may be impacted by the fact that Schrey et al. (2011) sampled D. frontalis during years of relatively low population densities (albeit following a significant outbreak in 2001–2002) whereas our specimens were collected at a time of considerable increase in trap captures in Louisiana and Alabama, and an active outbreak in Mississippi (B Sullivan, pers. comm., 2020).

During D. frontalis outbreak years, pulse irruptions have unusual population dynamics where beetles emerge from infested trees in several waves and inter-tree distance is an important determinant of attack success, with crowded healthy trees being at greater risk owing to the spill-over effects of aggregation pheromones (Hain et al., 2011). Conversely, in endemic (non-outbreak) situations, D. frontalis persists at low densities and generally must travel further after emergence to encounter trees already under stress (e.g., from lightning strike damage or disease). Indeed, seasonal differences within a given year can also impact dispersal behavior, whereby individuals have the greatest potential for long-distance dispersal in Fall and the lowest potential in mid-Summer, owing to temporal differences in fat content levels (Turchin & Thoeny, 1993). Ultimately, the impacts of inter- and intra-year differences in sampling across studies may at least partly account for contrasting conclusions about the nature and strength of spatial-genetic structure.

The idea that preexisting structure may have been “overwritten” either by a change in outbreak status (i.e., from endemic to epidemic) and/or on-going expansion processes underscores the importance of considering both and ecological and geographic context. For example, James et al. (2015) showed that while the spruce budworm (Choristoneura fumiferana), a cyclical irruptive pest in North America, exhibits weak spatial-genetic structure, there is considerable nuance to this general finding. For example, the legacy of high connectivity during past outbreaks can persist for some time afterwards, high admixture can be driven by migrant adults sourced from relatively few sites, and different life stages (i.e., adults vs. larvae) can concurrently exhibit different degrees of spatial genetic structure. Although the periodicity of spruce budworm outbreaks is considerably longer than that of D. frontalis (∼35 years vs. ∼5–10 years, respectively) and the synchrony of irruptive populations may occur at different spatial scales, there may be considerable similarity in the underlying processes that generate apparent weak spatial-genetic structure. Future work on D. frontalis could attempt to distinguish between migrant and resident individuals, and model range expansion dynamics and associated source–sink dynamics, explicitly incorporating temporal effects.