Panmixia across elevation in thermally sensitive Andean dung beetles

Abstract Janzen's seasonality hypothesis predicts that organisms inhabiting environments with limited climatic variability will evolve a reduced thermal tolerance breadth compared with organisms experiencing greater climatic variability. In turn, narrow tolerance breadth may select against dispersal across strong temperature gradients, such as those found across elevation. This can result in narrow elevational ranges and generate a pattern of isolation by environment or neutral genetic differentiation correlated with environmental variables that are independent of geographic distance. We tested for signatures of isolation by environment across elevation using genome‐wide SNP data from five species of Andean dung beetles (subfamily Scarabaeinae) with well‐characterized, narrow thermal physiologies, and narrow elevational distributions. Contrary to our expectations, we found no evidence of population genetic structure associated with elevation and little signal of isolation by environment. Further, elevational ranges for four of five species appear to be at equilibrium and show no decay of genetic diversity at range limits. Taken together, these results suggest physiological constraints on dispersal may primarily operate outside of a stable realized niche and point to a lower bound on the spatial scale of local adaptation.

A theory that is a notable exception in attempting to predict dispersal rates of organisms is Dan Janzen's seasonality hypothesis, which mechanistically links temperature variation across latitude with its consequences for dispersal and biogeographic patterns (Janzen, 1967). The hypothesis begins with the observation that (a) tropical ecosystems generally have less seasonal variation in temperature than temperate ecosystems and then assumes that (b) populations and species are adapted to the range of climates they experience. As a result, Janzen predicted that a tropical organism climbing a mountainside is more likely to encounter a physiologically challenging climate than a temperate organism would be, leading to selection against dispersal across elevational gradients (the principle also applies to downslope movement and the physiological challenges of crossing lowland valleys). A suite of downstream predictions naturally follow, for example, that tropical organisms should also have narrower elevational ranges, show increased population subdivision and a pattern of isolation by environment, and ultimately have higher speciation rates (Gadek et al., 2018;Ghalambor, Huey, Martin, Tewksbury, & Wang, 2006;Sheldon, Huey, Kaspari, & Sanders, 2018;Wang & Bradburd, 2014).
Despite reduced thermal tolerance and generally smaller elevational ranges of tropical species, surprisingly few researchers have examined dispersal across tropical gradients, though it is the mechanism in the seasonality hypothesis that links physiology with elevational ranges. Evidence from Andean sparrows (Cheviron & Brumfield, 2009) and subtropical forest trees in China (Shi, Michalski, Chen, & Durka, 2011) indicates gene flow across mountain sides can be sufficiently reduced to generate population genetic structure. Evidence for speciation across elevational gradients-a possible long-term consequence of population genetic structure-has been reported in tropical kingfishers (Linck, Freeman, & Dumbacher, 2019) and butterflies (Elias et al., 2009), though it appears rare in general (Caro, Caycedo-Rosales, Bowie, Slabbekoorn, & Cadena, 2013). To our knowledge, the only paper explicitly estimating effective migration rates across elevational gradients found reduced gene flow and greater population subdivision in tropical stream insects compared to their temperate relatives (Polato et al., 2018). However, the authors did not control for the possibility that dispersal is systematically reduced in the tropics for reasons other than physiological tolerance, for example, by comparing the relative contributions of isolation by environment and isolation by distance across latitude.
We used a densely sampled population genomic dataset to ask whether dispersal is reduced across elevation within the ranges in five species of dung beetles (Scarabaeinae) with well-characterized thermal physiologies and elevational distributions that conform to predictions of Janzen's seasonality hypothesis. We tested for reduced dispersal across elevational gradients relative to within elevational bands by estimating fine-scale population genetic structure, Wright's neighborhood size, and by using a Bayesian approach for describing the relative contribution of geographic distance and environmental distance to neutral genetic differentiation. To understand whether observed patterns were instead produced by recent population expansion (Gadek et al., 2018), we further asked whether elevational ranges are at equilibrium and tested for evidence of declining genetic diversity at elevational range limits.
Ectotherms with a global distribution, they are useful taxa for comparative studies of natural history, physiology, and population genetics. Our previous work using phylogenetically matched dung beetles from locations spanning 60° of latitude found thermal tolerance of species in the tribes Canthonini and Dichotomini generally increased with seasonality and was positively correlated with elevational range width-key predictions of Janzen's seasonality hypothesis (Sheldon & Tewksbury, 2014). We focused our current study on two species in the tribe Canthonini (Deltochilum speciosissimum, Deltochilum tessellatum), two species in the tribe Dichotomini (Dichotomius podalirius, Dichotomius satanas), and a single species (Eurysternus affin. flocossus) from a third tribe, Oniticellini. Species in the genus Deltochilum are ball-rolling dung beetles that feed and breed using dung or carrion, though at least one tropical species of Deltochilum is predatory and kills millipedes (Larsen, Lopera, Forsyth, & Génier, 2009). Beetles in the tribe Dichotomini excavate tunnels near or below dung deposits and then transport it underground to create brood balls for laying and incubating eggs, often closing off the tunnel's entrance (Hanski & Cambefort, 2014). Finally, Eurysternus spp. are unique among dung beetles in several aspects of their reproductive biology, including a "nuptial feast," or aggregation and consumption of dung balls prior to nesting, an inability to roll balls with their legs, multiple nests consisting of a shallow crater with several brood balls, nest care, and pair bond behavior in some species (Halffter, Halffter, & Huerta, 1980). We sampled beetles using pitfall traps baited with human dung (Davis, Scholtz, & Chown, 1999) along two elevational and two horizontal transects on the eastern slope of the Andes Mountains in Napo Province, Ecuador (Figure 1). Our two elevational transects spanned 730 to 1,175 m and 1,500 to 1,950 m in largely undisturbed humid forest, with four sampling localities spaced as close to 125 m apart as possible given local constraints of soil and topography. Our two horizontal transects were each approximately 2 km long, with one transect having four points and the other having two points.
Both horizontal transects were located in replicate corridors of montane forest at 2,150 m asl separated by a natural barrier, the Cosanga River. No species were shared between the two elevational transects, reflecting high beta diversity in Andean scarabs (Sheldon & Tewksbury, 2014). In all cases, we sampled nearly the complete elevational range of our focal species (K. Sheldon, unpublished data).
Following capture, all beetles were euthanized and stored in 95% EtOH prior to processing.

| Library preparation and DNA sequencing
We extracted genomic DNA from all samples using Qiagen DNeasy kits and the manufacturer's recommended protocol for insects, homogenizing a small amount of wing muscle tissue in pH 7.2 PBS prior to lysis (Qiagen). To prepare libraries for reduced representation high-throughput sequencing, we used a double restriction enzyme digest approach adapted from Gompert et al. (2014) and described below.
We placed 6 μl DNA from each sample (with a minimum concentration of 20 ng/μl) in separate wells of a chilled 384-well plates, filling the remainder with samples from a separate but related study.
We then added 2.6 μl of a restriction digest master mix consisting of 10× T4 buffer, 1M NaCL, 1 mg/ml BSA, water, and the MseI and EcoRI enzymes in a 1:0.52:0.52:0.2125:0.10:0.25 ratio. After sealing, vortexing, and centrifuging the plate, we incubated it for 2 hr at 37°C on a thermocycler with a heated lid, followed by 20 min at 65°C to denature the active enzymes. We next added 1.44 μl of an adaptor ligation master mix to each well. This master mix consisted of a

| Sequence assembly and variant calling
After initially demultiplexing our libraries using our unique barcodes and a Python script implemented in the first step of ipyrad (Eaton, 2014), we assembled sequencing reads and called and filtered variants for each species independently using the dDocent pipeline (v. 2.7.8), a set of bash wrapper scripts optimized for population genomics of nonmodel organisms requiring de novo assemblies (Puritz, Hollenbeck, & Gold, 2014 We then used VCFtools, vcffilter (1.0.0; https://github.com/vcfli b/vcflib), and the dDocent script dDocent_filters to filter this file to high-quality SNPs alone. We first dropped any individuals with missing data at more than 30% of loci, any loci missing data at more than 25% of individuals, all SNPs with a minimum minor allele frequency of 0.05 and a minimum minor allele count of 3, all SNPs with a quality (Phred) score of <30, and any genotypes with fewer than 5 reads across all individuals. We subsequently dropped sites with an allelic balance of >0.75 or <0.25, where the proportion indicates the ratio of reference allele reads to all reads; sites with reads from both strands; sites with a ratio of mapping qualities between reference and alternate alleles that were >0.9 or <1.05, and sites with a quality score >¼ below its depth. Finally, we identified sites with a depth exceeding 3× the square of the mean and removed any from this subset that lacked quality scores exceeding 2× their depth or were outliers, which we identified qualitatively based on the approximate point at which a histogram of mean depth across loci began to asymptote.
For the remainder of the paper, we refer to genotypes filtered in this way as our "primary" SNP datasets; we also generated a "secondary" SNP dataset without minor allele frequency filtering.

| Describing population genetic structure
To estimate population genetic structure within species across our gradients, we used two approaches appropriate for large SNP datasets. First, we performed principal component analysis (PCA) of genotypes using our primary SNP datasets using adegenet v. 2.1.1.
We plotted samples on principal component axes 1 and 2 and visually identified outliers representing misidentified samples, cryptic species, or highly divergent subpopulations. Second, we tested for fine-scale population subdivision within each species using fine-RADstructure v. 0.3.2, which uses haplotype data to identify a closest relative for each allele and then sums these data into coancestry similarity matrix among all individuals. We again used our primary SNP datasets, but dropped individuals identified as potentially belonging to different species through our PCA.

| Testing for isolation by environment
We quantified the relative contributions of geographic distance and environmental distance to observed genetic differentiation within all five species using BEDASSLE v. 1.5 (Bradburd, Ralph, & Coop, 2013). BEDASSLE is a Bayesian method that models allele frequencies in unlinked loci in a set of populations as a spatially correlated Gaussian process, in contrast to correlation-based Mantel and multiple matrix regression methods (Bradburd et al., 2013). Covariance is a decreasing function of both ecological and geographic distance, and parameters are estimated with a Markov Chain Monte Carlo simulation. We first filtered our primary SNP datasets for linkage disequilibrium using bcftools v. 1.9 and a maximum r 2 value of .1 and converted these files into allele counts per population (sampling locality) using adegenet v. 2.1.1 and custom R code (https:// github.com/elinc k/scarab_migra tion/blob/maste r/scrip ts/bedas sle.R). We calculated geographic distances among all sampling localities using the earth.dist() function in the R package fossil (Vavrek, 2011) and calculated environmental distances using the R package raster (Hijmans et al., 2019) and data for elevation, mean annual precipitation, and mean annual temperature from WorldClim2 (Fick & Hijmans, 2017). Given WorldClim2's 1 × 1 km resolution and reliance values interpolated from climate stations which may be few and far between in the tropics, we acknowledge this is an imperfect solution.
To partly counteract these issues and promote more efficient chain

| Demographic modeling
To test whether elevational ranges were in equilibrium or the product of recent population expansion, we compared simple demographic models of drift-mutation equilibrium and exponential population growth using an approximate Bayesian computation approach (Beaumont, Zhang, & Balding, 2002). We first calculated the site frequency spectrum (SFS) of a separate SNP dataset treated identically to our primary SNP dataset except for filters based on minimum minor allele count and minor allele frequency using the R package adegenet's glSum() function (Jombart, 2008). We next defined demographic models of a single population with either a single population-scaled nucleotide diversity parameter (our "null") model or both theta and an exponential population growth parameter α (our "growth") model using the coalescent simulator framework coala v. 0.5.3 in R (Staab & Metzler, 2016). Under the "growth" model, the population size changes by a factor e (− t) , where t is the time in generations since the growth has started. To approximate our empirical SNP data, we used scrm to simulate 50 three-nucleotide diploid loci for a sample size equivalent to each of our five species after filtering using the sequential coalescent with recombination model (Staab, Zhu, Metzler, & Lunter, 2015). We ran 100,000 simulations for each and calculated the resulting SFS. We then used the R package abc v. 2.1 to estimate parameters and perform model selection for each species (Csilléry, François, & Blum, 2012). We first performed leave-one-out cross-validation to evaluate the ability of ABC to distinguish among our models using tolerance rates of 0.01, 0.05, and 100 simulations. We then estimated parameters and performed model selection using the abc() and postpr() functions, implementing the rejection algorithm with a tolerance rate of 0.05 for both.

| Spatial patterns of genetic diversity and Wright's neighborhood size
To describe genetic diversity across the elevational distribution of our focal species, we calculated the population-scaled nucleotide diversity of each sampling locality using an estimator of theta ( = 4N e ) based on the mean homozygosity of gene frequencies in the R package pegas v. 0.1.1 (Paradis, 2010). We used our primary SNP dataset and calculated theta values for each RAD locus independently, plotting their full distribution. We examined the relationship between theta and the absolute distance in meters from mean sampling elevation (as a proxy for proximity to putative range limits) using linear mixed effects models with population as a fixed effect and tested for statistical significance using a likelihood ratio test. We additionally calculated Wright's neighborhood size for each species, a metric proportional to the average number of potential mates for an individual given its dispersal ability and defined as N w = 4 2 , where is mean parent-offspring distance, and is population density (Battey, Ralph, & Kern, 2020;Wright, 1946). To do so, we used Rousset's finding that the reciprocal of the slope of a linear regression of the natural log of geographic distance against F ST /(1−F ST ) is an estimator of 4 2 (Rousset, 1997). We analyzed pairwise F ST values from BEDASSLE's calculate.all.pairwise.Fst() function, and pairwise geographic distances calculated as described above (Bradburd et al., 2013). We transformed variables and ran simple linear regressions using custom R code (https://github.com/elinc k/scarab_migra tion/ blob/maste r/stats.R).

| Sampling and DNA sequencing
In total, we extracted DNA from 230 individuals of our five focal    (Figure 2b). Despite efforts to disentangle these variables in our sampling scheme, elevational distance was strongly correlated with geographic distance in our data (Pearson's r = .9129).

| Isolation by environment
Bayesian analysis of the relative contribution of isolation by environment (IBE) and isolation by distance (IBD) to genetic differentiation in BEDASSLE found little evidence of IBE in any focal species   [Weir and Hill, 2002]) likely impeded MCMC efficiency: In all but one species, acceptance rates fell to near 0 following burnin.

| Demographic modeling
Demographic model testing with approximate Bayesian computation suggests four out of five focal taxa have not experienced recent range expansion (Figure 4a) with moderate-to-strong support (Kass & Raftery, 1995). In contrast, the Bayes factor value for a comparison of population growth to the null model for E. affin. flocossus was 1.631, indicating low-tomoderate support. These patterns were reflected in the posterior probabilities of the population growth rate parameter in coala ("r" in Figure 4b).

| Spatial patterns of genetic diversity and wright's neighborhood size
Estimates of theta were uncorrelated to distance from mean sam-

| D ISCUSS I ON
Janzen's seasonality hypothesis predicts that organisms distributed in regions with relatively limited climatic variability-as is true in much of the tropics-will evolve a reduced physiological tolerance breadth for temperature (hereafter thermal tolerance; Janzen, 1967). In turn, narrow thermal tolerances should bias dispersal to occur most frequently between environments with similar temperature regimes, leading to neutral genetic differentiation associated with temperature and independent of geographic distance (Wang & Bradburd, 2014). Given recent evidence for the apparent ubiquity of isolation by environment in nature (Herrera, Medrano, & Bazaga, 2017;Manthey & Moyle, 2015;Sexton, Hangartner, & Hoffmann, 2014;Sexton et al., 2014;Shi et al., 2011;Wang & Bradburd, 2014;Weber, Bradburd, Stuart, Stutz, & Bolnick, 2017) and previous work demonstrating elevational ranges and thermal physiologies of our focal taxa broadly conform to predictions of the seasonality hypothesis (Sheldon & Tewksbury, 2014), we expected to see evidence of genetic differentiation across elevational ranges in the present study. We were therefore surprised to find no evidence of population genetic structure associated with elevation and little evidence of isolation by environment in any form (Figures 2 and 3). Reflecting these patterns, Wright's neighborhood sizes for all five species-proportional to the average number of potential mates for an individual given its dispersal ability-were effectively infinite relative to the geographic and environmental scale of our sampling (Figure 2c), indicating panmixia.
Importantly, this result appears to reflect long-term population dynamics rather than a temporary artifact of when we sampled.
Over evolutionary time scales, demographic processes may obscure the signal of selection against maladaptive physiological phenotypes, biasing inferences about the drivers of genetic differentiation.
For example, rapid population expansion across an environmental gradient may temporarily generate a signal of panmixia, before subsequent range contraction due to fitness costs, or the emergence Dichotomini, and Oniticellini are strong fliers with highly developed olfactory systems (Hanski & Cambefort, 2014), and they rely on an ephemeral, randomly distributed resource for feeding and reproduction. With low mammal population density in the tropical Andes (Jiménez et al., 2010;Lizcano, Pizarro, Cavelier, & Carmona, 2002;Ríos-Uzeda, Gómez, & Wallace, 2007) and putative variation in the nutritional content of mammalian dung (Hanski & Cambefort, 2014  to the validity of Janzen's hypothesis, as speciation is one of its potential consequences. Given previous work in this system demonstrating a greater number of cryptic species at lower latitudes (Gill et al., 2016), it seems likely intraspecific reproductive isolation is a major driver of their observed reduction in gene flow across elevation. However, it is at least as plausible that cryptic speciation was driven by the evolution of reproductive isolation through processes unrelated to physiology, such as genetic incompatibilities arising during a period of allopatry prior to secondary contact and range displacement. More detailed characterization of performance curves across temperature in these species could inform a mechanistic niche model (Kearney & Porter, 2009) and establish whether their current distributions approach abiotic tolerance limits.
Though Janzen (1967) did not directly address the consequences of the seasonality hypothesis for speciation, our findings raise two points related to the field that merit discussion. First, as discussed above, populations in the tropics separated by climatically challenging regions should be more isolated and diverge more rapidly than similarly separated populations of temperate congeners. While the absence of robust taxonomic, phylogenetic, and distributional data have largely precluded tests of diversification rates in tropical dung beetles (but see Davis & Scholtz, 2001;Davis et al., 1999), well-designed comparative phylogeographic studies across latitude targeting common species could help evaluate this hypothesis.
Second, strong constraints on thermal niche might trigger divergent selection and ecological speciation (Nosil, 2012;Schluter, 2001) in the event of niche expansion. However, the apparent stability of elevational ranges over evolutionary time scales suggests unusual circumstances might be required to facilitate this process, such as ecological release from a competitor, or a change in standing genetic variation enabling adaptation.
Our results suggest tropical dung beetles may have limited capacity to respond to climate change. Anthropogenic global warming may have a dramatic impact on tropical dung beetle communities assuming species cannot keep pace with their thermal niche through range shifts (Sheldon et al., 2011), though the magnitude and direction of range shifts remain unpredictable given variation in responses among closely related taxa and uncertainty in climate forecasts for tropical regions (Sheldon, 2019). Regardless, these predictions of community change ignore the potential for adaptation. Though experimental approaches to local adaptation remain the gold standard, data on the evolutionary ecology of range limits can inform predictions of the likelihood species can respond to climate changes at evolutionary timescales. Our data suggest this may be difficult for the species in the present study. In addition to evidence that most species' ranges are at equilibrium, the absence of any decay in genetic diversity away from the elevational range center ( Figure 5) suggests elevational ranges are not constrained by low population growth rates or limited standing genetic variation-that is, sourcesink dynamics driven by fitness costs or poor habitat quality at range limits (Moeller, Geber, & Tiffin, 2011).
While it is beyond the scope of our work to identify the ultimate drivers of elevational distributions, a pattern of high genetic diversity range-wide is consistent with range limits formed by either a biotic interaction or a swamping of maladaptive alleles (Angert & Schmeske, 2005;Bridle & Vines, 2007;Holt & Gomulkiewicz, 1997;Kirkpatrick & Barton, 1997;Moeller et al., 2011;Sexton et al., 2009).
The hypothesis of gene swamping seems particularly likely given steep environmental gradients in our study and the small spatial scales involved. If this process is indeed occurring, it points to an intriguing negative interaction between physiological traits, dispersal, and the scale of local adaptation. As broad thermal tolerances are selected against in populations in aseasonal environments, dispersal across elevation gradients is also reduced, leading to narrow elevational ranges. A high ratio of dispersal distances to the distance upper and lower range limits then exposes populations at different points of the gradient to gene swamping-further constraining the species from evolving into new niche space. In this scenario, where niche constraints are primarily evolutionary, a complex change in patterns of gene flow might be required to permit adaptation to novel climates. We encourage future workers to explore this rich area of inquiry.

ACK N OWLED G M ENTS
We thank Claire Winfrey for assistance in the field, C.C. Nice at Texas State for help with library preparation, and the staff of Universidad Regional Amazónica Ikiam molecular laboratory for hosting our work. This research was conducted under permit no.
MAE-DNB-CM-2017-0062 from El Ministerio de Ambiente to J. Celi, and was supported in part by an NSF IOS-1930829 to K.S.S.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
All code used in this study and an annotated digital notebook are available at https://github.com/elinc k/scarab_migra tion.