Flight loss is a common feature of upland insect assemblages, with recent studies detecting parallel wing reduction events across independent alpine lineages. However, the geographic scale over which such repeated evolution can operate remains unclear. In this study, we use genotyping-by-sequencing to assess the genomic relationships among vestigial-winged and full-winged populations of the widespread New Zealand stonefly Nesoperla fulvescens, to test for repeated wing loss events over small spatial scales. Biogeographic analyses indicate that alpine wing loss in this widespread species is restricted to a single, narrow mountain range. Intriguingly, our coalescent analyses indicate that upland vestigial-winged N. fulvescens populations are not sister to one another, suggesting wings have been lost independently in disjunct populations of this species, over a <30 km scale. Our results suggest that selection against flight above the alpine treeline can drive rapid and repeated adaptation even across narrow spatial scales. We propose that such repetitive processes may represent a far more pervasive feature of alpine insect adaptation than is currently recognized.
The extent to which similar selective forces can underpin predictable and repeated evolutionary shifts across different regions and lineages remains a key question in evolutionary biology (Lässig et al. 2017, Blount et al. 2018, Nosil et al. 2018, Reznick and Travis 2018, Waters and McCulloch 2021). As a case in point, the repeated reductive evolution of secondarily flightless lineages from flighted ancestors has long fascinated biological researchers (Darwin 1859, Darlington 1943). Despite the numerous evolutionary and ecological benefits conferred by flight, this dispersal capacity has been lost repeatedly across diverse insect and bird lineages (Wagner and Liebherr 1992, Roff 1994, Trautwein et al. 2012, Misof et al. 2014, Waters et al. 2020). Secondary flight loss is a particularly common feature of ‘insular’ insect assemblages on oceanic islands and alpine habitats (Roff 1990, Hodkinson 2005), and this dispersal reduction is thought to have played a substantial role in associated rapid divergence and speciation (Mitterboeck and Adamowicz 2013, Waters et al. 2020, Ortego et al. 2021, Salces-Castellano et al. 2021). Recent analyses further suggest that frequent exposure to strong winds in these habitats may play a key role in the loss of flight in weak-flying insects (McCulloch et al. 2019b, 2019a; Leihy and Chown 2020), supporting an evolutionary hypothesis originally proposed by Darwin (1859).
In addition to interspecific differentiation linked to flight loss, some individual insect taxa exhibit intraspecific dispersal polymorphisms, often involving energetic trade-offs between flight and other life-history traits (Harrison 1980, Roff 1990, Wagner and Liebherr 1992, Langellotto et al. 2000, Guerra 2011). Such dispersal-polymorphic lineages represent strong systems for assessing the evolutionary and ecological dynamics of insect flight loss (Van Belleghem et al. 2018, Suzuki et al. 2019, McCulloch et al. 2021b). For instance, the detection of multiple dispersal-limited populations within the broader geographic range of an otherwise dispersive taxon raises intriguing questions for evolutionary genomic analysis. Such systems promise to shed light on both the causes and consequences of dispersal reduction (Waters et al. 2020), and provide opportunities to test for the role of repeated selective processes in driving similar evolutionary shifts across independent populations and regions. Indeed, recent genomic studies of polymorphic taxa are beginning to uncover the genomic basis of repeated evolution, and are revealing the relative contributions of ‘repeated sorting’ (selection on standing variation) versus ‘molecular convergence’ (independent de novo mutations) to this process (Waters and McCulloch 2021).
New Zealand's (NZ's) species-rich plecopteran fauna (>100 endemic taxa) comprises a broad variety of full-winged, wing-reduced, and wingless taxa (McLellan 2006, McCulloch et al. 2016). Wing reduction is a particularly prominent feature of NZ's high-altitude grassland ecosystems (McCulloch et al. 2019b), which are dominated by flightless stonefly populations. While the vast majority of NZ's stonefly species are essentially monomorphic in their flight capability, two genera (Zelandoperla; Nesoperla) include taxa with wing polymorphisms (McLellan 2006), potentially providing strong systems for assessing the causes and consequences of insect flight loss (Waters et al. 2020). Interestingly, while the Zelandoperla fenestrata McLellan (Plecoptera: Gripopterygidae) complex includes wing-reduced alpine populations over a wide region of southern New Zealand (McLellan 1999, McCulloch et al. 2009, 2019b, 2021b; Dussex et al. 2016), wing-reduced forms of the widespread species Nesoperla fulvescens are relatively limited in their distribution, being apparently confined to a few alpine habitats of southern North Island (Fig. 1A and B; McLellan 1977, 2006). While independent evolution of flight loss might be anticipated over broad spatial scales (e.g., among widely-separated mountain ranges, 100s of kilometers apart; Suzuki et al. 2019; McCulloch et al. 2021b), the extent to which such parallelism can evolve over fine spatial scales (e.g., 10s of kilometers; disjunct alpine grassland habitats within single mountain ranges; McCulloch et al. 2019b) remains less clear.
Here we undertake biogeographic and genomic analyses of full-winged versus wing-reduced populations of the widespread N. fulvescens to test for repeated evolution of flightless alpine populations. Initially, we examine broad-scale distributional records to pinpoint the geographic distribution of wing-reduced populations. We then use genotyping-by-sequencing (GBS) to assess the genomic relationships among full-winged and vestigial-winged populations. If flight has been lost just once in this taxon, then disjunct vestigial-winged populations should together form a monophyletic flightless clade. Alternatively, if multiple wing loss events have occurred across independent upland populations (either through repeated sorting or molecular convergence; see Waters and McCulloch 2021), then flightless lineages should be polyphyletic. Based on emerging data from comparable systems (Hendrickx et al. 2015; McCulloch et al. 2019b, 2021b; Suzuki et al. 2019), we hypothesize that local flight loss may have evolved repeatedly across disjunct montane grassland habitats within this widespread wing-dimorphic species.
Materials and methods
Ecotype Distribution
To assess the geographic distribution of wing reduction in N. fulvescens, we measured hind tibia length (HTL) and forewing length for 141 adult specimens (across 59 localities) from the New Zealand Arthropod Collection (Auckland), Canterbury Museum (Christchurch), and from the personal collection of Ian Henderson (Massey University). We classified each specimen as either full-winged (wing:HTL > 1.8) or vestigial-winged (wing:HTL < 0.3), and used NZ Topo Map ( http://www.topomap.co.nz) to confirm whether specimens were collected above or below the alpine treeline.
Stonefly Sampling and DNA Extraction
We sampled for full-winged and vestigial-winged N. fulvescens from three populations across the Tararua Range from 21st to 25th November 2017 (Field, Dundas, Ruapae; Fig. 2A). Additional samples of full-winged N. fulvescens were collected from two other lower North Island localities from 22nd to 25th October 2020 ( Supp Table S1 (ixac027_suppl_supplementary_tables.docx) [online only]). We note that this species tends to be rare in lowland sites, with patchy distributions, though large populations have been recorded at some sites in the South Island (Winterbourn 2010). Both nymphs and adult specimens were collected by hand from underneath stones in small alpine streams, and from immediately adjacent stony/muddy habitat. Additionally, two full-winged, flighted specimens were collected opportunistically from the surfaces of high-elevation pools. Samples were preserved in 80% ethanol. Late instar nymphs were classified as either full-winged or vestigial-winged based on clear differences in wing bud development (see Veale et al. 2018). We extracted DNA from the head and pronotum of stoneflies using DNeasy kits (Qiagen), following the manufacturer's protocols.
Fig. 1.
(A) An alpine grassland tributary of the Hector River (elevation 1350 m above sea-level) on Mt Field, Tararua Range (North Island, New Zealand), where vestigial-winged Nesoperla fulvescens were collected; (B) A vestigial-winged N. fulvescens male from Mt Field; (C) National collection records of N. fulvescens, illustrating the narrow distribution of vestigial-winged ecotypes (Tararua Range).
Mitochondrial Sequencing
We amplified a 634-bp portion of the mitochondrial cytochrome oxidase I (COI) gene from 26 individuals, 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). To explore relationships among closely related COI haplotypes we constructed TCS haplotype networks, using PopART 1.7 (Leigh and Bryant 2015).
Genotyping-by-sequencing
We conducted GBS on 51 N. fulvescens samples from eight locations across the lower North Island ( Supp Table S1 (ixac027_suppl_supplementary_tables.docx) [online only]). The GBS library preparation and sequencing were conducted at AgResearch, Invermay, New Zealand following the methods of Elshire et al. (2011), with modifications as outlined in Dodds et al. (2015). Briefly, samples were indexed, and the restriction enzyme PstI was used to fragment the DNA. Samples were pooled into a single library, which was column purified. Size selection (220–520 bp) was conducted using a Pippen Prep (Sage Science). The library was sequenced as 2/3 of a lane of an Illumina NovaSeq6000 (1 × 101bp), utilizing v1.5 chemistry.
We used fastqc 0.11.9 (Andrews 2010) to evaluate the quality of the data. Variants were called using Stacks 2.58 (Rochette et al. 2019). The detailed SNP calling procedure can be found in Appendix S1 (ixac027_suppl_supplementary_appendix.docx). Initially, adapters were removed using cutadapt 3.5. We then demultiplexed the raw reads using the process_radtags command in Stacks. All reads were then trimmed to a common length (60 bp). We removed low-quality reads, and reads that did not contain the enzyme recognition site. As there is no reference genome available for N. fulvescens, we de novo assembled the cleaned reads, using the denovo_map.pl pipeline. The parameters used in de novo assembly can significantly impact downstream analyses (Mastretta-Yanes et al. 2015). We therefore trialed a range of different –M and –n parameters (–M = 2 and –n = 2; –M = 3 and –n = 3; –M = 4 and –n = 4; –M = 5 and –n = 5), and tested how they impacted our results.
We further filtered these data using the populations module of Stacks. We only retained biallelic loci that were found in at least 70% of samples, and discarded loci with a heterozygosity above 0.65 (as these loci are potentially paralogs [see O'Leary et al. 2018]).
Genetic Structure
When assessing population structure and demographic history we excluded outlier loci (see the Outlier Analysis section, below), as these analyses rely on selectively neutral markers. We assessed genetic similarity among samples using principal component analysis (PCA) in the R package ADEGENET 2.1.1 (Jombart 2008). Pairwise FST values were calculated among the three neighboring populations from the Tararua Range using ARLEQUIN 3.5 (Excoffier and Lischer 2010). We assessed the significance of the FST values using 20,000 permutations, with Bonferroni corrections applied to account for multiple tests.
Fig. 2.
(A) Map of Nesoperla fulvescens sampling sites from in southern North Island, New Zealand. (B) Principal component analyses, based on 39,642 SNP markers, demonstrating the genomic divergence among different N. fulvescens populations. Individuals are colored by geographic location (from A). Vestigial-winged ecotypes are marked with an asterisk.
Population structure across the Tararua Range was further assessed in STRUCTURE v2.3.4 (Pritchard et al. 2000), using the admixture model with alleles frequencies correlated. The analysis was run for 500,000 generations, after an initial burn-in of 50,000 iterations. We ran five independent runs from K = 1 to K = 6. These runs were combined in CLUMPAK (Kopelman et al. 2015), and the most likely number of clusters assessed using the Delta K approach (Evanno et al. 2005).
We used BA3-SNPs (Wilson and Rannala 2003, Mussmann et al. 2019) to assess the rates of contemporary migration among the three upland N. fulvescens populations, and to test for evidence of asymmetric gene flow. The –a and –f parameters were adjusted to 0.40 and 0.03, so that acceptance rates for proposed changes for each mixing parameter were between 0.2 and 0.6 (following the BA3-SNPs user manual). We ran the program for 10,000,000 generations, with the first 1,000,000 generations discarded as burn-in. To assess chain convergence, we ran the analysis ten times, with each run started with a different random seed.
Demographic History
To explore the evolutionary relationships among the three Tararua populations in more detail—and to test for repeated wing loss events—we used coalescent analysis with DIYABC-RF 1.0.14 (Collin et al. 2021). This program implements a model-based inference to simulate coalescence of observed populations given different demographic scenarios. Simulated datasets are then compared to the observed data set using the machine learning Random Forest approach (Breiman 2001), to determine the best supported model and to estimate posterior probability. This approach thus allowed us to compare several different historical models, and to estimate divergence timing among the populations. We compared three different potential demographic scenarios. In Scenario One, the full-winged Ruapae population is sister to vestigial-winged populations from Field and Dundas, which diverged from one another relatively recently. In Scenario Two, the vestigial-winged Field population is sister to the Ruapae and Dundas populations. In Scenario Three, the vestigial-winged Dundas population is placed sister to the Field and Ruapae populations. Assuming the ancestral N. fulvescens population was full-winged (which seems likely, given that secondary flight loss is a common feature of alpine insect evolution; McCulloch et al. 2021b), Scenario One would imply a single loss of flight, whereas Scenarios Two and Three would both imply multiple wing reduction events.
We removed SNPs that were not present in at least one individual from each group, and set the minimum minor allele frequency to 0.10 (as recommended by Cornuet et al. 2008). We used broad priors with a uniform distribution (102–107 for each parameter), as the specific divergence times and population sizes are unknown. We performed 500,000 simulations per scenario to produce posterior distributions, with each scenario considered equally probable at the outset. All summary statistics were used in the Random Forest analyses. To assess the reliability of the observed summary statistics, we ran a principal component analysis to test whether observed summary statistics of the most likely scenario were surrounded by the simulated summary statistics (see Collin et al. 2021). To select the most likely scenario, and to estimate the timing of divergence among lineages, we ran 5,000 Random Forest trees, using the default number of noise variables (five). The divergence timing of populations was estimated based on a one generation per year lifecycle, as N. fulvescens is univoltine (Scarsbrook 2000).
Outlier Analysis
We used BayeScan 2.1 (Foll and Gaggiotti 2008) to search loci potentially associated with wing reduction. Specifically, we categorized each stonefly as either full-winged or vestigial-winged, and tested for loci that were significantly differentiated among these two groups. We used the default parameters, but with prior odds of 100 to minimize the number of false positives (because of the large number of markers in our data set). Loci with q-values of less than 0.05 were considered as potential outliers.
In addition to excluding outlier loci when assessing population structure and demographic history, we also tested whether any of these outliers were linked to genes with known roles in wing reduction or other alpine adaptations (e.g., delayed emergence (McCulloch and Waters 2018a), increased body size (McCulloch and Waters 2018b), longevity (McLellan 1999), and fecundity (Guerra 2011), and olfactory shifts (Neupert et al. 2022)). We applied two different approaches: 1) using the BLASTx function on the UniProt Knowledgebase ( https://www.uniprot.org/blast/), with an E-value threshold of <1e–3, and 2) aligning outlier loci to the annotated Z. fenestrata reference genome ( https://www.genomics-aotearoa.org.nz/data; McCulloch et al. 2021a), and searching for genes of known function within 20 kb of the outlier.
Results
Ecotypic Distribution
Distributional and morphological analyses based on specimens retrieved from national invertebrate collections (incorporating 141 adult N. fulvescens specimens across 59 localities) confirmed that vestigial-winged ecotypes of this species have been detected only at high elevation grassland sites across a narrow region of the North Island's Tararua Range (Fig. 1C; Supp Table S2 (ixac027_suppl_supplementary_tables.docx) [online only]). These alpine wing-reduction localities also represent the only sites where N. fulvescens has been detected above the alpine treeline ( Table S2 (ixac027_suppl_supplementary_tables.docx) [online only]). Similarly, our field sampling detected vestigial-winged N. fulvescens ecotypes from just two high elevation grassland sites (>1,300 m) within the Tararua Range (Fig. 2; Field, Dundas). All samples from these two sites were vestigial-winged, apart from a single full-winged specimen from each locality, which were collected from the surface of small alpine pools as adults ( Supp Table S1 (ixac027_suppl_supplementary_tables.docx) [online only]). Samples from all other sites were full-winged, including those from the forested, subalpine (1,200 m) Ruapae site.
Genetic Structure
Mitochondrial sequencing revealed 13 discrete COI haplotypes across the N. fulvescens samples. No ambiguous sites or stop codons were evident, suggesting all sequences were of mitochondrial origin. Phylogeographic structuring was evident, with samples broadly grouping by site (Fig. 3). However, several samples from Field (including the full-winged adult) were more closely related to samples from Ruapae than the remaining Field samples, and two haplotypes were shared by samples from Field and Dundas. The vestigial-winged and full-winged samples were broadly represented by distinct haplogroups, although some vestigial-winged samples from Field and Dundas were mitochondrially more similar to full-winged samples (Fig. 3).
Genotyping-by-sequencing yielded an average of 1,480,265 reads per individual. Different SNP filtering schemes resulted in very similar PCAs ( Supp Fig. S1 (ixac027_suppl_supplementary_figure_s1.jpeg) [online only]), so all additional analyses were based on a single filtering scheme (–M = 2, –n = 2). Following SNP calling (and removal of outlier SNPs), we compared 51 N. fulvescens samples across 39,642 SNPs. Samples had a mean coverage of 29× across all individuals and loci, and the proportion of missing data in the final dataset was 9.4%.
Fig. 3.
COI haplotype network of Nesoperla fulvescens from five locations. Each circle represents a haplotype, and circles are scaled by size according to the number of sequenced individuals per haplotype. Haplotypes are colored by locality (see key). Uninterrupted lines in the network represent single step mutations. Small open circles indicate hypothetical intermediate (unsampled) haplotypes.
Principal component analysis revealed strong population structuring across the Tararua Range, with samples generally clustering by locality (Fig. 2), and significant pairwise FST values among all three Tararua localities (P < 0.0001). However, samples from the vestigial-winged population from Dundas were genomically closely related to the geographically adjacent (5 km northeast) full-winged population from Ruapae (Fig. 2), rather than to the morphologically-similar (but geographically remote) vestigial-winged population from Field Peak (approximately 30 km southwest). Similarly, the pairwise FST value between Ruapae and Dundas samples (0.076) was much lower than that between Ruapae and Field (0.153), or between Dundas and Field (0.135). Interestingly, the two anomalous flighted specimens opportunistically sampled from Field and Dundas localities did not cluster with the vestigial-winged samples from their respective sites, but were closely related to one another (Fig. 2).
STRUCTURE results (Fig. 4A) were broadly congruent with the results of the PCA. The delta K method indicated that K = 3 was the most likely number of population clusters, though K = 2 also received some support ( Supp Fig. S2 (ixac027_suppl_supplementary_figure_s2.jpeg) [online only]). K = 2 separated the Field population from the Dundas and Ruapae populations, whereas K = 3 broadly differentiated all three populations from one another (Fig. 4A).
Estimates of contemporary migration rates were identical across all runs. Migration rate estimates were very low across all population pairs, with estimated rates of migration from full-winged to vestigial-winged populations only slightly larger than the corresponding reverse estimates (Fig. 4B).
Demographic Analysis
Coalescent analyses suggest that the vestigial-winged Field population diverged initially from the Ruapae and Dundas populations, with Ruapae and Dundas populations diverging from one another more recently (Scenario 2; Fig 5). This scenario received 85.3% of the classification votes, and had a very high posterior probability (0.931; Fig. 5). By contrast, Scenario 1 received fewer than 4% of the classification votes (Fig. 5). A PCA confirmed that the observed data were contained within the simulated summary statistics for scenario 2 ( Supp Fig. S3 (ixac027_suppl_supplementary_figure_s3.jpeg) [online only]). Coalescent analyses of the favored scenario suggest that the Field N. fulvescens population diverged from Ruapae and Dundas populations ca. 153 kya (95% CI: 278 kya–71 kya). The Ruapae and Dundas populations, by contrast, are estimated to have diverged substantially more recently, around 81 kya (95% CI: 164 kya–36 kya).
Outlier Analysis
We identified 64 outlier SNPs that were significantly differentiated between full-winged and vestigial-winged N. fulvescens ( Supp Table S3 (ixac027_suppl_supplementary_tables.docx) [online only]). Three of these outlier SNPs mapped to genes of known function in the NCBI database, but none of these genes have previously reported roles in insect wing development ( Supp Table S3 (ixac027_suppl_supplementary_tables.docx) [online only]). Of the 21 outlier SNPs that successfully mapped to the Z. fenestrata reference genome, many were located within 20 kb of annotated genes, although none of these identified genes were identified as having documented roles in wing development ( Supp Table S3 (ixac027_suppl_supplementary_tables.docx) [online only]).
Discussion
Our distributional analyses confirm that flight loss is rare across the range of N. fulvescens, having been detected only in a few alpine localities within the Tararua Range. Intriguingly, our new demographic analysis indicates that the two vestigial-winged populations from this range are not each other's closest relatives. Assuming the common ancestor of these populations was full-winged, this finding suggests that wing reduction has occurred independently in two adjacent but disjunct populations.
While mtDNA analysis might suggest that the two vestigial-winged populations are more similar to one another than either is to the full-winged Ruapae population (Fig. 3), this suggestion is not supported by more powerful GBS analysis (Figs. 2, 4, 5). We note that phylogenetic resolution of recently-diverged lineages provided by single locus approaches is typically limited relative to that provided by multilocus genome-wide approaches (see McCulloch et al. 2019b, 2021b). Regardless, any potential mito-nuclear discordance (between mtDNA and GBS datasets) could stem from a variety of factors including selection, incomplete lineage sorting, or perhaps introgression among populations (see Toews and Brelsford 2012). The potential for such introgression is indicated in the current study by our migration analysis which suggests low levels of gene flow among lineages (Fig. 4). Additional sampling and screening (particularly of full-winged lowland population) promise to shed further light on the evolutionary history of these lineages.
Although our demographic analyses suggest that wings have been lost independently in the Field and Dundas N. fulvescens populations, it remains possible that this apparently repeated flight loss has been underpinned by selection on standing genetic variation (see Van Belleghem et al. 2018, Haenel et al. 2019, Lai et al. 2019), rather than due to independent de novo mutations in each population (see Hart et al. 2018, Zhang et al. 2021). This proposition may be circumstantially supported by the fact that the vestigial-winged phenotype is apparently restricted to the Tararua Range. The proposed localized ‘spread’ of wing reduction among neighboring alpine Tararua populations could perhaps have involved transportation of ‘flight-loss’ alleles via flighted individuals (i.e., transporter hypothesis; Schluter and Conte 2009). However, further sampling and additional genomic data (e.g., whole-genome sequencing) are required to better assess the role that standing genetic variation has played in repeated flight loss in N. fulvescens.
Recent studies have shown that the repeated evolution of flight loss can be common across broad spatial scales (e.g., among populations on distinct mountain ranges over 100s of kilometers; Suzuki et al. 2019, McCulloch et al. 2021b). Our new results for Nesoperla, however, suggest that such repeated processes may operate even over fine spatial scales (<30 km; see also McCulloch et al. 2019b). These findings of apparently replicated alpine adaptation over even small geographic scales raise intriguing questions about the potential extent of deterministic, repeated evolutionary processes in upland ecosystems more broadly, particularly in tectonically active, young mountain ranges (Hörandl et al. 2005, Dunning et al. 2013, Wallis et al. 2016, Koot et al. 2020, Waters et al. 2020). We propose that such repeated evolution may be a far more pervasive alpine evolutionary process than is currently recognized.
Fig. 4.
(A) STRUCTURE plots (K = 2 and K = 3) of Nesoperla fulvescens from three sites across theTararua Range, based on 39,642 SNP markers. Each vertical bar represents an individual, with color indicating the inferred genomic cluster.The anomalous full-winged individuals from Field and Dundas are indicated with stars (B)The estimated proportion of migrants (+95% HPD intervals) to and from each N. fulvescens population on theTararua Range, derived from BA3-SNPs.
Fig. 5.
Coalescent analyses support parallel flight loss events in Nesoperla fulvescens across different portions of theTararua Range. Competing demographic scenarios (A–C) were assessed with DIYABC-RF, with the level of support noted under each scenario (time is not to scale). Stars indicate inferred wing-reduction events, assuming a winged ancestor. Divergence estimates for the scenario with the highest support (outlined) were calculated with DIYABC-RF.
The apparent rarity of full-winged N. fulvescens from locations above the alpine treeline suggests that harsh environmental conditions (e.g., strong winds), coupled with limited habitat size, may render these grassland ecosystems unsuitable for insect flight, and lead to strong selection for wing reduction (McCulloch et al. 2019b, 2019a; Leihy and Chown 2020; Foster et al. 2021). Interestingly, the two anomalous full-winged N. fulvescens adults collected from ponds above the alpine treeline are genomically-distinct from the corresponding vestigial-winged samples, but closely related to one another. These findings suggest that these flighted adults had likely dispersed (post-metamorphosis) into alpine habitats from neighboring (lower-elevation) forested source populations. If these individuals did disperse from lower-elevation populations, this implies they are strong flyers, a proposition also tentatively supported by our migration analysis (Fig. 4). Conversely, the apparent absence of vestigial-winged N. fulvescens populations from below the alpine treeline suggests strong selection against wing-reduced individuals in forested habitats. Our results thus further demonstrate how the alpine treeline can shape the distributions of winged versus wing-reduced insects (Foster et al. 2021).
The Tararua Range is geologically young, with recent research suggesting that upland alpine habitats in this region have existed for fewer than 500 ka (Craw et al. 2019). The young age of this range implies that the upland N. fulvescens populations likely evolved flight loss in the late Pleistocene, a proposition supported by our temporal genomic analyses, which suggest that the two vestigial-winged populations sampled from the Tararua Range diverged from one another less than 200 kya. Additionally, the detection of multiple mtDNA haplotypes (Fig. 3) within the flightless populations of N. fulvescens (Field: 6 vestigial haplotypes; Dundas: 4 vestigial haplotypes), and the apparent absence of shared haplotypes between ecotypes, is broadly consistent with ecotype divergence prior to the last glacial maximum (ca. 20,000 years ago; Fig. 5).
Intriguingly, the genetic variant/s underpinning flight loss do not appear to have spread beyond the Tararua Range. As a consequence, N. fulvescens is apparently restricted to habitats below the alpine treeline throughout the rest of its broad distribution, in spite of the availability of ample upland habitat (e.g., Ruahine Range, central North Island). The relative rarity of wing reduction in N. fulvescens contrasts with the situation for the similarly polymorphic Z. fenestrata, which comprises numerous wing-reduced upland populations (McCulloch et al. 2019b, 2021b). The contrast among these two taxa may suggest that wing reduction has evolved relatively recently in the case of N. fulvescens, and has yet to spread across its geographic distribution. Alternatively, gene-flow among N. fulvescens populations may be particularly limited—perhaps reflecting the relatively strict habitat requirements of this species, and its relatively fragmented distribution (see Phillipsen et al. 2015).
We identified more than 50 outlier SNPs that were significantly differentiated between full-winged and vestigial-winged ecotypes. While none of these SNPs were close to recognized insect wing development genes, it should be noted that the molecular mechanisms underpinning wing loss among different insect lineages remain poorly understood (Saastamoinen et al. 2018). One outlier SNP was closely linked to muskelin, a gene that plays a role in fecundity, and we thus speculate that this locus might potentially be linked to life history divergence among ecotypes (McLellan 1999, McCulloch et al. 2021a). Additionally, some of the outliers detected may potentially be artifacts of geographic structure, rather than being linked to genes involved in alpine adaptation. While we tested for outliers across a large suite of SNPs, the markers we used covered only a fraction of the genome (see Hoban et al. 2016). Thus, additional sampling and genomic analysis are clearly required to better elucidate the key gene/s involved in wing reduction and other alpine adaptations.
Overall, our findings suggest that natural selection can drive repeated flight loss in alpine insects (McCulloch et al. 2019b, 2021a), over even small spatial scales. Future analyses promise to further enhance our understanding of the ecological genomic bases, and phylogeographic dynamics of such repetitive local evolutionary processes (Waters et al. 2020).
Acknowledgments
We thank Marianne Vogel for assisting with 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’. We wish to acknowledge the use of New Zealand eScience Infrastructure (NeSI) high performance computing facilities ( https://www.nesi.org.nz). 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. This research was supported by Marsden contracts UOO2016 and UOO1412 (Royal Society of New Zealand). Specimens were collected under DOC permit 63798-FAU.
Author Contributions
JW and BF conceptualised the study and carried out the sampling. GM conducted the laboratory work. GM and LD analysed the data. GM 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 PRJNA825078. Cytochrome oxidase I sequences have been deposited in GenBank (Accession numbers OP565413-OP565436). Demultiplexed GBS data and aligned COI sequences have been added to Datadryad ( doi:10.5061/dryad.w9ghx3fsn).
