Range-restricted upland taxa are prone to population bottlenecks and thus typically have low genetic diversity, making them particularly vulnerable to environmental change. In this study, we used a combination of genotyping-by-sequencing (10,419 SNPs) and mitochondrial COI sequencing to test for population genetic structure within the narrow-range flightless sub-alpine stonefly Zelandoperla maungatuaensis Foster. This species is restricted to only a handful of upland streams along a 4 km stretch of the isolated Maungatua range in southeast New Zealand. We identified striking genetic structure across the narrow range of Z. maungatuaensis, with three deeply divergent allopatric lineages detected. These distinct lineages likely diverged in the early-mid Pleistocene, apparently persisting in separate microrefugia throughout subsequent glacial cycles. Our results illustrate how secondary flight loss can facilitate insect diversification across fine spatial scales, and demonstrate that intraspecific phylogenetic diversity cannot necessarily be predicted from range-size alone. Additional demographic analyses are required to better understand the conservation status of these divergent Z. maungatuaensis lineages, and to assess their potential susceptibility to climate change and other anthropogenic impacts.
Understanding the forces driving biological diversification remains one of the fundamental goals of evolutionary biology (Mayr 1942, Sobel et al. 2010). The evolution of insect flight, some 400 million years ago (Ma) (Misof et al. 2014) was likely a key driver of their diversification, allowing species to disperse, and exploit novel habitats (Dudley 2000). Recent studies, however, indicate that secondary flight loss may further enhance these rates of diversification (Salces-Castellano et al. 2021), primarily by limiting gene-flow across populations (Ortego et al. 2021a). These high levels of intraspecific diversification may also be correlated with increased speciation rates (Ikeda et al. 2012, Waters et al. 2020, Fenker et al. 2021, Gálvez-Reyes et al. 2021).
The evolutionary trajectories of low-dispersal taxa can be significantly impacted by ‘vicariant’ earth-history processes, such as mountain building (Hoorn et al. 2013, Craw et al. 2016) and glaciation (McCulloch et al. 2010, Wallis et al. 2016, Ortego et al., 2021b). Nondispersive freshwater species may be particularly vulnerable to such abiotic processes, as such taxa have limited ability to disperse beyond their catchments (Avise et al. 1987, Craw et al. 2016, McCulloch et al. 2019a). While recent studies have highlighted the biological effects of physical processes in shaping the genetic structuring of such lineages over broad geographic scales (e.g., among major drainage systems; Mayden 1988; Waters et al. 2001; Finn et al. 2016; and among mountain ranges; McCulloch et al. 2010, Ortego et al. 2021b), relatively little is known about how such processes operate across finer spatial scales (e.g., among adjacent river tributaries within mountain ranges; Wishart and Hughes 2003).
Range-restricted taxa (e.g., endemic to a single mountain range) represent challenging systems for conservation and biogeographic research. Such narrow-range endemics can potentially arise via a combination of intrinsic (e.g., limited dispersal ability; strict habitat requirements) and/or extrinsic (geographic barriers; competitive and density-dependent ecological processes; Waters et al. 2020) factors. In many cases, range-restricted species have experienced major population bottlenecks and associated rapid genetic drift, and are thus expected to exhibit only low levels of genetic diversity (e.g., Frankham, 1998; but see Binks et al. 2015). Such low diversity may be particularly pervasive for isolated upland species (Finn et al. 2014, Jordan et al. 2016, Hinojosa et al. 2019), especially those vulnerable to climate-driven reductions (Muhlfeld et al. 2011, Giersch et al. 2015). However, some narrow-range upland taxa present notable exceptions to these expectations, with substantial divergence between neighbouring populations (e.g., Wishart and Hughes 2003).
Many of New Zealand's upland stoneflies (order Plecoptera) are flightless, and have restricted geographic distributions (McLellan 2006, McCulloch et al. 2017). Such lineages potentially present strong model systems for assessing the effects of flightlessness on genetic diversity (Waters et al. 2020). As a case in point, the recently described flightless Maungatua stonefly, Zelandoperla maungatuaensis Foster (Plecoptera: Gripopterygidae), is restricted to a handful of adjacent streams (within just a few km of one another) along a single mountain range in southeast New Zealand (Fig. 1; Foster et al. 2020). Like many of New Zealand's flightless stoneflies, this species is found almost exclusively above the alpine treeline (McCulloch et al. 2019a, McCulloch et al. 2019b). As Z. maungatuaensis diverged phylogenetically from its closest relative (the lowland, full-winged Z. denticulata McLellan; western South Island) some 2 million years ago, a single (ancient) loss of flight ability has been suggested for this upland lineage (Foster et al. 2020).
In this study we test for intraspecific genetic diversity within the narrow-range ‘relict’ Z. maungatuaensis. Specifically, we use population genomic analyses to assess phylogenetic diversity and structuring in this endemic species. We hypothesize that, due to strict habitat requirements and limited dispersal capacity, there will be little gene flow among neighbouring subpopulations, leading to the development of genetic structuring among the different stream populations. We assess the distribution of this intraspecific genetic diversity under both evolutionary and conservation frameworks.
Methods
Sampling
Z. maungatuaensis was originally described from two unnamed subalpine streams on the eastern face of the Maungatua range. We conducted additional sampling across the Maungatua range, from February to November 2020, to better assess the geographic distribution of the species. Nymphs were collected from under small stones in riffle zones, and preserved in 100% ethanol.
Mitochondrial Sequencing
DNA was extracted from stonefly head tissue using DNeasy kits (Qiagen). We amplified a 651-bp portion of the mitochondrial cytochrome oxidase I (COI) gene from 33 individuals across 8 streams (3–6 specimens per stream), following the protocols of McCulloch et al. (2009). We tested for evidence of pseudogenes by examining whether there were any ambiguous sites or stop codons in the sequence (see Zhang and Hewitt 1996).
All new COI sequences obtained in this study were deposited in GenBank (GenBank accession numbers OM802728-OM802761).
The newly amplified COI sequences were aligned with two additional Z. maungatuaensis COI sequences (accession numbers MK942108 and MK942109) in Geneious 11.1.5 ( https://www.geneious.com/; Kearse et al. 2012), using the MUSCLE plugin (Edgar 2004). Phylogenetic relationships were reconstructed in MrBayes 3.2.7a (Ronquist et al. 2012), with Zelandoperla agnetis McLellan (GQ414594) and Zelandoperla denticulata (GQ414593) used as outgroups. We ran four Markov chains for five million generations, with chains sampled every 200 generations. The first 2,000 trees were discarded as burn-in. Tracer v1.7.0 (Rambaut et al., 2018) was used to confirm that the parameters had all converged, and to ensure that the effective sample size was greater than 200 for each of the priors. We estimated the divergence among the mitochondrial clades using standard insect calibrations (3.54% per million years; Papadopoulou et al. 2010).
Genotyping-by Sequencing
We conducted genotyping-by-sequencing (GBS) to investigate fine-scale genetic structuring among Z. maungatuaensis samples. The GBS library preparation and sequencing was conducted at AgResearch, Invermay, New Zealand following the methods of Elshire et al. (2011), with modifications as outlined in Dodds et al. (2015). We used the same 35 Z. maungatuaensis specimens that were used for COI sequencing. Briefly, DNA was fragmented using the restriction enzyme PstI. Sample specific barcodes were added, and the samples were pooled into a single library. The library was column purified, then size selected (220–520 bp) using a Pippen Prep (Sage Science). We used single end sequencing (1 × 101bp), with the library sequenced as 2/3 of a lane of Illumina NovaSeq6000, utilizing v1.5 chemistry.
We evaluated the quality of the data using fastqc 0.11.9 (Andrews 2010). Variants were called using Stacks 2.58 (Rochette et al. 2019). The detailed SNP calling procedure can be found in Appendix S1. Briefly, we removed adapters using cutadapt 3.5 before demultiplexing the raw reads using the process_radtags command in Stacks. All reads were trimmed to a common length (60 bp). Low-quality reads, and reads that did not contain the enzyme recognition site, were removed. As there is no Z. maungatuaensis reference genome, we de novo assembled the cleaned reads, using the denovo_ map.pl pipeline with default parameters. Following genotyping, we used the populations module of Stacks to further filter the data. Specifically, we only retained loci that were found in 80% of samples, and discarded any loci containing SNPs with a heterozygosity above 0.65 (as these loci are potentially paralogs [see O'Leary et al. 2018]).
We assessed population structuring using principal component analysis (PCA) in the R package ADEGENET version 2.1.1 (Jombart 2008; R Core Team 2020). We examined the phylogenetic relationships among populations by constructing a maximum likelihood phylogeny (using IQ-Tree v2.1.1; Minh et al. 2020). ModelFinder (Kalyaanamoorthy et al. 2017) was used to select the best fit model for our data, with an ascertainment bias applied (as our data set lacks invariant sites). We assessed node support with 5,000 ultrafast bootstrap approximations (Hoang et al. 2018). The phylogeny was rooted using the midpoint rooting method.
Pairwise FST values among streams were calculated in arlequin 3.5 (Excoffier and Lischer 2010). The significance of these values was assessed with 10,000 permutations, with Bonferroni corrections applied to account for multiple tests. We tested for an association between genetic distance (linearized FST) and geographic distance (the overland distance between sites) using a Mantel test (Mantel 1967) in Genepop 4.7.5 (Rousset 2008).
Results
Z. maungatuaensis nymphs were collected from eight streams along the eastern face of the Maungatua range ( Supp Table S1 [online only]; Fig. 2a). We did not detect the species in any of the streams on the southern or western faces of the range (Fig. 2a). Our surveys thus confirm that Z. maungatuaensis has a very narrow geographic distribution (less than 4 km across).
Sequencing revealed four discrete COI haplotypes across the eight streams. No ambiguous sites or stop codons were evident, suggesting all sequences were of mitochondrial origin. Phylogenetic analysis revealed three well-supported mitochondrial clades corresponding to distinct regions: a southern clade, a central clade, and a northern clade (Fig. 3a). The three clades did not overlap geographically (Fig. 2a). The northern and central clades were recovered as sisters, and differed by 1.0%. A significantly larger divergence was evident between the central and southern clade (4.3%), and the northern and southern clade (4.2%). Under standard insect calibrations rate, our results suggest that the three clades began to diverge about 1.2 Ma. Almost no genetic diversity was evident within each of the clades, with all but the northern clade containing only a single haplotype.
GBS yielded average of 3,353,899 reads per individual. Following SNP calling, 35 Z. maungatuaensis samples were compared across 10,429 SNPs. The proportion of missing data in the final dataset was 6.7%, with a mean coverage of 51× across all loci and individuals. Population genomic analyses revealed significant population structure across the Maungatua range, with individuals again clustering into the distinct northern, central, and southern lineages (Fig. 2b), and significant pairwise FST values were evident among populations from the three geographic regions (Table 1). In contrast to the COI analysis, the central and southern clades were recovered as sister taxa (Fig. 3b). Limited population genetic structure was also evident within the northern and central clades ( Supp Fig. S1 [online only]), though pairwise FST values among different streams within clades were relatively small, and not statistically significant (Table 1). We detected weak, but significant isolation-by-distance across the streams (r2 = 0.167; P = 0.031; Supp Fig. S2 [online only]), though this pattern was primarily driven by the differentiation among the distinct geographic lineages.
Discussion
Evolution of Fine-scale Diversity
In the current study we detect striking genetic substructuring within a narrow-range endemic insect over just a few kilometres of upland habitat in southern New Zealand. Specifically, strong genome-wide (GBS) differentiation among three allopatric lineages of Z. maungatuaensis on the Maungatua range is mirrored by deep mtDNA divergences (up to 4.3%), with concordant groupings detected for both marker sets. These results illustrate how wing loss can drive rapid diversification across very small spatial scales. Our results imply that the neighbouring lineages diverged in the early-mid Pleistocene (ca. 1 Ma), and likely persisted locally in separate microrefugia throughout repeated glacial cycles. While such microrefugial dynamics might be anticipated over broad geographic scales (e.g., separate mountain ranges, hundreds or thousands of kilometres apart; Wallis et al. 2016, King et al. 2020); to discover them operating locally over kilometre scales is intriguing. Along similar lines, while occasional contact zones involving divergent lineages have been reported for widespread upland taxa (Weston and Robertson 2015, King et al. 2020, Ortego et al. 2021b), finding such deep genetic structure within a narrow-range endemic taxon is unexpected.
Several recent phylogeographic studies of upland insect lineages over fine spatial scales have revealed substantial genome-wide structure, sometimes linked to disjunct stream habitats (Dussex et al. 2016, McCulloch et al. 2019a, McCulloch et al. 2021). However, these previous local examples typically involved recently-diverged lineages lacking corresponding mitochondrial divergence (Dussex et al. 2016, McCulloch et al. 2021). Additionally, in these previous cases, the local genetic structure stems from repeated flight loss events, with neighbouring flightless populations having been independently derived (McCulloch et al. 2019a, Suzuki et al. 2019). However, in the current study there are no closely-related flighted lineages likely to have seeded such independent dispersal-reduction events (see Waters et al. 2020), with flightless Z. maungatuaensis clearly representing a monophyletic assemblage, with its closest relative (the full winged Z. denticulata) found more than 200 km away (Foster et al. 2020).
While both GBS and mtDNA yielded three congruent groupings of Z. maungatuaensis samples, the relationships among these key groupings varied across marker sets. Such mitonuclear discordance is a common feature of isolated freshwater populations (Wallis et al. 2017), and could stem from a variety of factors, including past introgression or stochastic lineage sorting (see Toews and Brelsford 2012). Additional possible factors contributing to mitonuclear discordance could include sex-biased dispersal (e.g., Yannic et al. 2012), Wolbachia infection (Gompert et al. 2008), positive selection (e.g., Sunnucks et al. 2017), and genetic incompatibilities (e.g., Arntzen et al. 2009). Regardless, the evolution and preservation of deep genetic divergence among allopatric alpine stream populations apparently attests to long histories of population isolation and stability across different portions of the Maungatua range. In the context of comparable studies of range-restricted upland invertebrates (e.g., Finn et al. 2014, Jordan et al. 2016, Hinojosa et al. 2019), the current data seem exceptional in terms of both the timeframe (early-mid Pleistocene), and the very narrow spatial scale involved (a few km). Broadly, these findings suggest that phylogeographic diversity is not necessarily predictable from geographic range size alone (Binks et al. 2015). Indeed, additional factors such as habitat stability and connectivity may play crucial roles in facilitating genetic diversification over even the narrowest of geographic scales (Wishart and Hughes 2003).
Geological Controls on Stonefly Evolution
The Maungatua Range that hosts Z. maungatuaensis is a relatively isolated high topographic feature, and the nearest mountains of similar altitude are more than 30 km away (Fig. 1). This topographic isolation has almost certainly facilitated initial speciation of these flightless insects. The topography is relatively young, in that the whole area has risen from below sea level since the Miocene (last 20 million years), and is dominated by tectonic uplift rather than erosion (Youngson and Craw 1996). Within this region, the Maungatua Range has risen as an elongated domal feature above the surrounding landscape in even more recent times, through the Pleistocene. This more localised uplift has been associated with Pleistocene–Holocene fault activity along the margins of what is now the Taieri Basin, a Pleistocene subsidence feature ( Supp Fig. S3 [online only]; Litchfield 2001). Pleistocene uplift occurred along the northern side of North Taieri Fault and the Maungatua Range rose above, e.g., the adjacent Barewood Plateau because of locally enhanced uplift rates. The elongated domal shape of the Maungatua Range uplands has evolved progressively towards the northeast and southwest through the Pleistocene and this may be on-going. The three allopatric lineages of Z. maungatuaensis have developed separately on this dome as it evolved.
Fig. 2.
Striking genetic substructuring across the narrow geographic range of Zelandoperla maungatuaensis. (a) Collection localities (coloured circles) for Z. maungatuaensis across the Maungatua range. White crosses indicate sites where Z. maungatuaensis has not been found. (b) Principal component analysis, based on 10,429 genome-wide SNPs, illustrating the genetic differentiation among Northern, Southern, and Central Z. maungatuaensis lineages.
Fig. 3.
Phylogenetic reconstruction reveals three well-supported Zelandoperla maungatuaensis clades. (a) Bayesian maximum likelihood consensus phylogeny illustrating the relationships among Z. maungatuaensis lineages based on the mitochondrial COI gene. Posterior probability values are noted above each node. Outgroups (Zelandoperla agnetis and Zelandoperla denticulata) are excluded for diagrammatic clarity. (b) Midpoint rooted maximum likelihood phylogeny illustrating genome-wide relationships among Z. maungatuaensis lineages, based on 10,429 SNP markers. IQTREE ultrafast bootstrap values ≥ 90 are indicated at each node.
The approximately 4 km range of Z. maungatuaensis incorporates several parallel, high-gradient (unnamed) streams draining the southeast slopes of the Maungatua range, at and above the alpine treeline (550–700 m a.s.l.). These streams are typically connected to each other only via the downstream swampy Taieri plain, a habitat highly unsuitable for stoneflies, rendering many of these Z. maungatuaensis populations effectively isolated from one another. These distinctive landscape features thus appear to explain the evolution of multiple isolated stonefly lineages over a narrow spatial scale. In a few cases, the detection of genetically-related lineages among adjacent tributary populations (e.g., F, G & H; C & D; Fig. 2a) likely represents gene flow facilitated by continuous riverine connectivity (i.e., where tributaries merge at intermediate elevations). Such gene flow may have been enhanced during glacial periods, as geographic ranges of the lineages likely extended to lower elevations in response to associated cooling and loss of forest cover. By contrast, some apparent genetic connectivity among adjacent drainages that have no high-elevation wet connections might reflect occasional small-scale (e.g., 50–100 m) overland dispersal of upland adults, particularly in adjacent headwater regions (e.g., D & E; B & A; see Fig. 2a; Craw et al. 2007).
Table 1.
Pairwise FST values among collection localities for Zelandoperla maungatuaensis (see Fig. 2a for locality details)
Conservation Implications
The detection of substantial genome-wide and mitochondrial divergence within the narrow geographic distribution of a range-limited taxon raises potentially challenging questions for conservation, especially when management resources are already stretched (Lester et al. 2014, Drinan et al. 2021). Indeed, the levels of mtDNA divergence detected here are comparable to those often found between distinct insect species pairs (e.g., McCulloch et al. 2010, Allegrucci et al. 2014, Wallis et al. 2016, Alfaro et al. 2018). As such, future systematic analyses should reassess the taxonomic status of these genetically isolated subpopulations of Z. maungatuaensis. While these populations are partly protected as a result of their distributional overlap with the Maungatua Scenic Reserve, additional demographic analyses are required to better assess their conservation status and potential susceptibility to climate change and other potential anthropogenic impacts.
More broadly, our data imply that caution is needed when making conservation and management (e.g., threat-ranking; Drinan et al. 2021) decisions that focus primarily on the geographic area occupied by rare taxa (see Rattis et al. 2018). Specifically, we suggest that geographic range is not always a useful surrogate for genetic diversity (Hanson et al. 2017). Indeed, our study highlights the potentially complex microrefugial dynamics of low-dispersal species over even small geographic scales (Binks et al. 2015), where local genetic diversification can substantially exceed that found in more dispersive taxa over relatively vast spatial scales (e.g., Allibone and Wallis 1993, Waters et al. 2020).
Acknowledgments
We thank Brodie Foster for assisting with early fieldwork. We also thank Tracy van Stijn, Rudiger Brauning, and John McEwan for conducting the GBS analysis, supported by the MBIE program C10X1306 “Genomics for Production and Security in a Biological Economy”. This research was supported by the Royal Society of New Zealand (Marsden contract UOO2016). Specimens were collected under DOC permit 59871-RES. We acknowledge the use of New Zealand eScience Infrastructure (NeSI) high performance computing facilities. New Zealand's national facilities are provided by NeSI and funded jointly by NeSI's collaborator institutions and through the Ministry of Business, Innovation & Employment's Research Infrastructure programme. URL https://www.nesi.org.nz.
Author Contributions
GM and JW conceptualised the study and carried out the sampling. GM, GK, and LD conducted the laboratory work, and analysed the data. GM and JW wrote the manuscript, with significant input from all authors.
Data Availability
All GBS data used in this study can be found in the NCBI Sequence Read Archive PRJNA809479.
