Introduction
The ability to correctly identify individual specimens to a species level is important for water quality and biodiversity evaluations of aquatic systems (Jones Reference Jones2008). Separating closely related species is typically completed through a comparison of diagnostic characteristics, yet these comparisons are often limited by the amount of morphological divergence visible by eye (Bickford et al. Reference Bickford, Lohman, Sodhi, Ng, Meier and Winker2007). Most descriptions of stoneflies are based primarily on adult males because their genitalia are distinctive enough among species to be diagnostic (Szczytko and Kondratieff Reference Szczytko and Kondratieff2015). Accurate morphological descriptions of females and larvae remain difficult for many species and, in some instances, are impossible with standard light microscopy (Zhou et al. Reference Zhou, Jacobus, DeWalt, Adamowicz and Hebert2010). Moreover, there are few stonefly taxonomists, and at times they are unable to identify specimens due to the lack of objective diagnostic traits (Beaty Reference Beaty2015).
Modern molecular techniques, coupled with open access data (e.g., GenBank), can help overcome taxonomic issues such as insufficient knowledge of morphological variability for traditional species descriptions (Hausmann et al. Reference Hausmann, Parisi and Sciarretta2016). Molecular data can provide a less labour-intensive and faster method for making male–female and adult–larvae associations of Plecoptera where only subtle morphological differences exist (Miller et al. Reference Miller, Alarie, Wolfe and Whiting2005; Mynott et al. Reference Mynott, Webb and Suter2011; Hausmann et al. Reference Hausmann, Godfray, Huemer, Mutanen, Rougerie and van Nieukerken2013). Molecular data of species for any combination of male, female, and larvae remain highly conserved within species (Calleja et al. Reference Calleja, Peña, Ugalde, Ferreiro, Marco and Garesse1993). Therefore, gene sequences from authoritatively identified adult males can be used as a reference for associating females and larvae. Integrative molecular and morphological approaches can delineate minor morphological differences that conform to diagnostic characters (McDaniel and Shaw Reference McDaniel and Shaw2003). Several recent studies have used molecular data for the purposes of evolutionary systematics of stoneflies (Mynott et al. Reference Mynott, Webb and Suter2011; Gill et al. Reference Gill, Kondratieff, Stark and Sandberg2015a, Reference Gill, Sandberg and Kondratieff2015b; Vitecek et al. Reference Vitecek, Vinçon, Graf and Pauls2017; South et al. Reference South, Skinner, Dewalt, Kondratieff, Johnson and Davis2021). Studies have also used DNA to associate different stonefly life stages or to associate males with females (Miller et al. Reference Miller, Alarie, Wolfe and Whiting2005; Mynott et al. Reference Mynott, Webb and Suter2011; Grubbs et al. Reference Grubbs, Metzger and Wei2020). The increasing use of gene sequence data for stoneflies has allowed for an integrated approach that combines morphology and molecular information (Vitecek et al. Reference Vitecek, Vinçon, Graf and Pauls2017). Much of the previously completed work focusing on DNA-based life-stage association relies on the mitochondrial cytochrome c oxidase subunit 1 (mtCO1) gene (Caterino and Tishechkin Reference Caterino and Tishechkin2006).
Despite decades of technological advances with light and electron microscopy, describing the adult female stage and associating females with males remain difficult with some stonefly genera. One example is with Nearctic Leuctra Stephens, 1836 (Plecoptera: Leuctridae). Thirty-one Leuctra species are currently recognised from the eastern Nearctic region, with most found in the Appalachian Mountains and unglaciated regions of the southeastern United States of America (DeWalt et al. Reference DeWalt, Maehr, Neu-Becker and Stueber2020). Harper and Harper (Reference Harper and Harper1997) placed most of the species known at that time (n = 26) into five proposed species groups: L. biloba Claassen, 1923 group, L. duplicata Claassen, 1923 group, L. ferruginea (Walker, 1852) group, L. grandis Banks, 1906 group, and L. tenuis (Pictet, 1841) group. Two species were left unassigned, but the five species described since 1997 have each been placed into a group (Grubbs and Sheldon Reference Grubbs and Sheldon2008; Grubbs Reference Grubbs2010; Harrison and Stark Reference Harrison and Stark2010). The most recent Nearctic taxonomic key (Hitchcock Reference Hitchcock1974) for this genus includes only 17 species, with descriptions of larva and females provided for only seven of these species. Taxonomic information in general on females is absent or minimal at best (e.g., Hitchcock Reference Hitchcock1974; Harper and Harper Reference Harper and Harper1997) due to conservative morphological characters of females and larvae characteristics. Although the females of some Leuctra species are described in taxonomic literature (e.g., Claassen Reference Claassen1923; Poulton and Stewart Reference Poulton and Stewart1991; Grubbs Reference Grubbs2010; Harrison and Stark Reference Harrison and Stark2010), they are often difficult to distinguish due to conservative morphology of the subgenital plate, setal patterns, and shape of the internal spermathecae (Harper and Harper Reference Harper and Harper1997). This is especially problematic in the more species-rich Appalachian Mountains, where as many as 7–10 species may inhabit the same drainage. Despite this, confident identification of females and larvae is critical for understanding the diversity, ecology, and conservation of species in the region.
Only one phylogenetic treatment of Nearctic Leuctra has been undertaken (Grubbs et al. Reference Grubbs, Metzger and Wei2020), and it is small in scope, with only 12 of the 31 known species. The two partial phylogenetic treatments of Palearctic Leuctra are comparatively larger in scope. Gattolliat et al. (Reference Gattolliat, Vinçon, Wyler, Pawlowski and Sartori2016) analysed an mtCO1 sequence library for 90 species from Switzerland, including 37 Leuctra species. Vitecek et al. (Reference Vitecek, Vinçon, Graf and Pauls2017) used mitochondrial and nuclear genes to assess the phylogenetic relationships of 18 species of the L. inermis Kempny, 1899 species group.
This study had two main objectives. The first objective was to incorporate mtCO1 sequence data with tree- and genetic distance–based phylogenetic approaches to associate putative females with males collected in 2019 for several species of Leuctra known from a southern Appalachian Highlands landscape, western North Carolina. The second objective was to build upon the phylogenetic framework established in Grubbs et al. (Reference Grubbs, Metzger and Wei2020) to a more robust assessment of phylogenetic relationships amongst eastern Nearctic Leuctra species and the proposed species groups.
Materials and methods
Study region
Fieldwork occurred in Mount Mitchell State Park and in the adjacent Appalachian Ranger District of Pisgah National Forest (Fig. 1) in North Carolina. Mount Mitchell is the high point in North America east of Hudson Bay and the Mississippi River Basin. The mountain is located within the ancient Blue Ridge Physiographic region that extends from Pennsylvania to northern Georgia, United States of America, and it varies from narrow ridges to high peaks in mountainous areas. Both landscapes are nested within the United States Environmental Protection Agency Level III Blue Ridge Ecoregion 65. A total of 44 sites were selected to represent a broad range of stream sizes, from small seeps to the South Toe River, and along an elevation gradient that spanned from 831 to 1983 m (from 2726 to 6506 ft) above sea level (Supplementary material, Table S1). Most sites were accessed along backcountry trails and United States Forest Service roads (Fig. 1). Coordinates for each site were taken with handheld geographic positioning system units and were subsequently projected to check for accuracy using Acme Mapper 2.2 (WGS-84; Acme Laboratories, Vancouver, British Columbia, Canada; http://mapper.acme.com). All geographic coordinate data are represented in decimal degrees. Dot distribution maps were prepared using ArcMap 10.7 (Esri, Redlands, California, United States of America).
Fig. 1. Collection sites during April–October 2019 from Mount Mitchell State Park and Pisgah National Forest, North Carolina, United States of America. Map inset of North Carolina in upper right shows the Appalachian Ranger District, Pisgah National Forest in dark green. Sites are represented by black stars.
Field sampling
Fresh adult specimens of Leuctra were collected in April, June, July, August, and October 2019. We used beating sheets to dislodge specimens from riparian vegetation during cool weather, sweep nets on warmer days, hand-picking from emergent leaf packs, rocks, trees, and bridges, and light trapping with a standard ultraviolet light. Leuctra specimens were preserved in 95% ethanol for molecular analyses.
Microscopy and digital imaging
Males were identified using an Olympus SZ61 stereomicroscope (Olympus, Shinjuku City, Tokyo, Japan) and several resources (Claassen Reference Claassen1923; Hanson Reference Hanson1941; Hitchcock Reference Hitchcock1974; Harper and Harper Reference Harper and Harper1997; Grubbs Reference Grubbs2015). Subgenital plates for all females sequenced in this study (n = 272) were digitally imaged using a Leica MZ16 stereomicroscope (Leica, Wetzler, Germany) equipped with a JVC KY-F75U camera (JVC, Yokohama, Japan). Images were layered using the Auto-Montage Pro programme (Syncroscopy, Cambridge, United Kingdom).
Molecular methods
Total DNA was extracted from two adult legs and associated thoracic tissue per individual from nine regional species collected in 2019 (L. ferruginea (Walker, 1852), L. grandis Banks, 1906, L. mitchellensis Hanson, 1941, L. sibleyi Claassen, 1923, L. tenella Provancher, 1878, L. tenuis (Pictet, 1841), L. triloba Claassen, 1923, L. truncata Claassen, 1923, and L. variabilis Hanson, 1941) (Supplementary material, Table S2). We extracted DNA using a DNeasy Blood and Tissue Kit (Qiagen, Inc., Hilden, Germany) following the manufacturer’s instructions. We amplified the mtCO1 gene through a polymerase chain reaction using the primers HCO2198 (5′-TAAACTTCAGGGTGACCAAAAAATCA-3′) and LCO1490 (5′-GGTCAACAAATCATAAAGATATTGG-3′; Folmer et al. Reference Folmer, Black, Hoeh, Lutz and Vrijenhoek1994). The polymerase chain reaction was performed with 10 µL total (5 µL GoTaq® Colorless Master Mix (Promega Corporation, Madison, Wisconsin, United States of America), 0.5 µL of each of 10-µM solutions of HCO2198 and LCO1490, and 4 µL DNA). The polymerase chain reaction thermal regime consisted of an initial denaturation at 95 °C (1 minute), 35 cycles at 95 °C (20 seconds each), annealing at 50 °C (20 seconds), and an extension period at 72 °C (1 minute), followed by a final extension period at 72 °C (5 minutes). Polymerase chain reaction products were verified on 1% agarose gel (0.35 g agarose and 35 mL 1 × Tris-acetate-EDTA buffer; Thermo Fisher Scientific, Waltham, Massachusetts, United States of America) using a FisherBiotech FB-SB-710 Electrophoresis System (Thermo Fisher Scientific) and subsequently visualised compared to a 1-kB DNA ladder with an IO Rodeo Large Blue LED transilluminator (IO Rodeo Inc., Pasadena, California, United States of America). Polymerase chain reaction products were sequenced at the North Carolina State University Genomics Lab (Raleigh, North Carolina). Full collection data, photographs, sequencing primers, and all sequences have been added to a private data set on the Barcode of Life Database (BOLD; University of Guelph, Guelph, Ontario, Canada). Barcode of Life Database IDs can be seen in Supplementary material, Table S1.
To add more information for phylogenetic analyses, six species were also collected in 2019 (L. alta James, 1974, L. carolinensis Claassen, 1923, L. moha Ricker, 1952, L. rickeri James, 1976, L. schusteri Grubbs, 2015, and L. tenuis). Cytochrome c oxidase subunit 1 sequences were treated as above, and high-quality mtCO1 sequences for L. alexanderi Hanson, 1941, L. duplicata Claassen, 1923, L. ferruginea, L. rickeri, L. sibleyi, and L. tenuis were downloaded from BOLD (Supplementary material, Table S2). Phred quality scores indicated the quality of each base in the sequences. All sequences were trimmed to a score of 30 using 4peaks, version 1.8 (Nucleobytes, Aalsmeer, The Netherlands). Edited sequences were assembled and aligned using Muscle in MEGA X (Kumar et al. Reference Kumar, Stecher and Tamura2016).
Molecular analyses
Molecular data were analysed with tree- and genetic distance–based phylogenetic methods. Overall tree topology and a preliminary assessment of putative female–male associations were visualised with a neighbour-joining tree constructed in MEGA X. Nodal support was assessed using 1000 bootstrap replicates (Felsenstein Reference Felsenstein1985). With one exception, males of each species clustered as monophyletic units. The neighbour-joining tree was used to determine the species identity of each female. To confirm the placement of females, two more robust phylogenetic approaches were used. A maximum-likelihood tree was generated using RAxML with the autoMRE setting. Although nodal support was again assessed using 1000 bootstrap replicates, the autoMRE setting allowed this analysis to finish earlier if 1000 replicates were unnecessary. A Bayesian-inference tree was also constructed using MrBayes MPI to further assess species relationships (Huelsenbeck et al. Reference Huelsenbeck, Larget, Miller and Ronquist2002) with the following settings: sequence data, GTR_I_G model, default Markov chain Monte Carlo commands (four chains, three heated), default number of runs (n = 2), stationary nucleotide frequency, set burn-in of 25%, random starting tree, sampling frequency 1000 generation, and a stop value between 0.01 and 0.05. Approximately 65 000 000 generations were required. A consensus tree was generated following the stop rule value.
Pairwise divergence distances were calculated in MEGA X using the Kimura two-parameter model of nucleotide substitution with the complete codon position deletion option. We compared pairwise divergence values among males and putative females of the same species (intraspecific) and between all species (interspecific). Minimum and maximum intraspecific and interspecific values were generated as more reliable indicators of genetic distinctiveness than referring only to mean values (Meier et al. Reference Meier, Shiyang, Vaidya and Ng2006; Mynott et al. Reference Mynott, Webb and Suter2011). Automatic barcode gap discovery (Puillandre et al. Reference Puillandre, Lambert, Brouillet and Achaz2012) was performed on the pairwise-distance data generated in MEGA X to complement the phylogenetic approach. Automatic barcode gap discovery creates operational taxonomic units based on comparisons of pairwise distances. Relative gap width was set at 1.5, the Kimura two-parameter model was applied, and all other settings were default (Puillandre et al. Reference Puillandre, Lambert, Brouillet and Achaz2012). This method has been used successfully for stoneflies with new species descriptions (Teslenko et al. Reference Teslenko, Palatov and Semenchenko2019) and for supporting operational taxonomic unit delineations based on mtCO1 sequences (Morinière et al. Reference Morinière, Hendrich, Balke, Beermann, König and Hess2017; Grubbs et al. Reference Grubbs, Metzger and Wei2020).
Results
In total, 360 Leuctra males were collected from Mount Mitchell streams in 2019, and all could be readily identified to species level using traditional morphological characteristics (i.e., dorsal processes and paraprocts). We collected 477 females, and most could not be identified to the species level with any confidence because subgenital plate features were similar across several species (Figs. 2–4).
Fig. 2. Leuctra ferruginea group females, abdominal terminalia, ventral view: A, L. ferruginea OTU 1, Setrock Creek; B, L. ferruginea OTU 1, Setrock Creek; C, L. ferruginea OTU 2, Lower Creek; D, L. ferruginea OTU 2, tributary Lower Creek; E, L. truncata, Big Lost Cove Creek; and F, L. truncata, Neals Creek. OTU, operational taxonomic unit.
Fig. 3. Leuctra tenuis group females, abdominal terminalia, ventral view: A, L. tenella, Balsam Spring; B, L. tenella; seep; C, L. tenuis, South Toe River; D, L. tenuis; South Toe River; E, L. triloba, tributary Setrock Creek; F, L. triloba, seep; G, L. variabilis, Lower Creek; and H, L. variabilis, seep.
Fig. 4. Leuctra biloba group (A–B), L. grandis group (C–F), and undetermined Leuctra (G–H) females, abdominal terminalia, ventral view: A, L. mitchellensis, seep; B, L. mitchellensis, seep; C, L. grandis, Big Lost Cove Creek; D, L. grandis, Lower Creek; E, L. sibleyi, Left Prong South Toe River; F, L. sibleyi, Big Lost Cove Creek; G, Undetermined Leuctra operational taxonomic unit, Lower Creek; and H, Undetermined Leuctra operational taxonomic unit, Lower Creek.
We sequenced or gathered from BOLD a total of 401 sequences. These consisted of 162 determined males and L. grandis females collected in 2019, as well as BOLD downloads, and 239 undetermined females collected in 2019 (Supplementary material, Table S1). With two exceptions (see below), all valid species were supported by high bootstrap values (80–100%) and posterior probabilities (0.80–1.00). The maximum-likelihood and Bayesian-inference phylogenies displayed similar topologies, but none strongly supported the species group hypotheses of Harper and Harper (Reference Harper and Harper1997).
Morphology-based L. grandis female determinations were confirmed through molecular methods with successful pairing with males of the same species. The maximum-likelihood and Bayesian-inference trees associated females with males with very high nodal support for at least 16 species, including two L. ferruginea operational taxonomic units and one undetermined Leuctra operational taxonomic unit. The two L. ferruginea operational taxonomic units (L. ferruginea OTU 1 and L. ferruginea OTU 2) overlapped in space (Fig. 5). In addition, six of the seven other species that were included to provide additional information also had well-supported clades.
Fig. 5. Distribution of the two Leuctra ferruginea operational taxonomic units from Mount Mitchell State Park and Pisgah National Forest, North Carolina, United States of America. Map inset of North Carolina in upper right shows the Blue Ridge Ecoregion in grey and the Appalachian Ranger District, Pisgah National Forest in dark green.
Of particular importance was the clade for L. carolinensis. This is a regional species with a “North Carolina, Black Mountains” type locality (Claassen Reference Claassen1923) that was collected in east Tennessee, United States of America, only approximately 40 km from Mount Mitchell but not during this study. No Mount Mitchell females grouped with the L. carolinensis individuals (Figs. 6–7).
Fig. 6. Maximum likelihood tree for 18 eastern Nearctic Leuctra species made in RAxML. Support greater than 70% is represented on nodes. Horizontal triangle width is a function of the number of sequences per species. Scale bar represents substitution rate. Canadian and United States of America postal codes: Canada, ON = Ontario, NS = Nova Scotia; United States of America, IL = Illinois, KY = Kentucky, NC = North Carolina, NY = New York, SC = South Carolina, and TN = Tennessee. Proposed Leuctra species groups are represented by colour: orange = L. grandis group, grey = L. duplicata group, green = L. biloba group, black = L. tenuis group, red = L. ferruginea group, and blue = ungrouped species.
Fig. 7. Bayesian-inference tree for 18 eastern Nearctic Leuctra species made in MrBayes. Support greater than 0.7 is represented on nodes. Horizontal triangle width is a function of the number of sequences per species. Scale bar represents substitution rate. Canadian and United States of America postal codes: Canada, ON = Ontario, NS = Nova Scotia; United States of America, IL = Illinois, KY = Kentucky, NC = North Carolina, NY = New York, SC = South Carolina, and TN = Tennessee. Scale bar represents estimated substitution rate. Proposed Leuctra species groups are represented by colour: orange = L. grandis group, grey = L. duplicata group, green = L. biloba group, black = L. tenuis group, red = L. ferruginea group, and blue = ungrouped species.
Genetic divergence and automatic barcode gap discovery
There was overlap between maximum intraspecific and minimum interspecific divergence values for several species. Intraspecific divergence values ranged from 0% to 18.6% (Table 1), whereas minimum interspecific values ranged from 2.2% to 15.5% (Table 2). Although intraspecific values of at most 3% occurred only with a combination of six species and operational taxonomic units (Table 1), there was a general, albeit variable, relationship with the number of individuals sequenced (Fig. 8). A larger number of individuals sequenced (e.g., L. alta, L. triloba; Fig. 8) were associated with a broader range of intraspecific divergence values. The notable exception was for L. sibleyi, with a high intraspecific divergence value (18.6%; Table 1) that was attributed to one individual (Fig. 8). Interestingly, the haplotype outlier was from Mount Mitchell and not the BOLD sequence of a Nova Scotia, Canada male.
Table 1. Mean and range of intraspecific divergence values of mtCO1 sequences for 16 eastern Nearctic Leuctra species, including two L. ferruginea operational taxonomic units, plus three undetermined females (L. operational taxonomic unit). Leuctra species are organised by groups as proposed in Harper and Harper (Reference Harper and Harper1997). “F” and “M” refer to numbers of female and male mtCO1 sequences included in this study, respectively. Species in bold type were collected from Mount Mitchell, North Carolina, United States of America in 2019. Standard Canadian and United States postal codes were used in the table. OTU, operational taxonomic unit
* All from BOLD sequences from the United States of America – Tennessee.
† All from BOLD sequences from the United States of America – New York.
‡ Includes one BOLD sequence from Canada – Ontario.
§ Includes one BOLD sequence from the United States of America – Illinois.
║ Includes one BOLD sequence from Canada – Nova Scotia.
# Includes one BOLD sequence each from Canada – Nova Scotia and the United States of America – Tennessee. ABGD, automatic barcode gap discovery.
Table 2. Minimum interspecific divergence values of mtCO1 sequences for at least 16 eastern Nearctic Leuctra species, including two L. ferruginea operational taxonomic units and three undetermined females (L. operational taxonomic unit), collected from Mount Mitchell, North Carolina, United States of America in 2019. Species are organised alphabetically by groups as proposed in Harper and Harper (Reference Harper and Harper1997). Species names are abbreviated: alex (L. alexanderi), mitc (L. mitchellensis), dupl (L. duplicata), ferr (L. ferruginea OTU 1 and OTU 2), rick (L. rickeri), trun (L. truncata), alta (L. alta), gran (L. grandis), sibl (L. sibleyi), caro (L. carolinensis), schu (L. schusteri), tene (L. tenella), tenu (L. tenuis), tril (L. triloba), vari (L. variabilis), and moha (L. moha). OTU, operational taxonomic unit
Fig. 8. Box and whisker jitter plot of intraspecific genetic distance values for each individual species unit. Species are arranged alphabetically within each proposed species group. Each dot represents a distance value between two unique sequences. The vertical line, rectangle, and horizontal line represent the median, interquartile range, and entire range of data, respectively.
Higher intraspecific divergence values were also associated with a higher number of automatic barcode gap discoveries (Table 1). Automatic barcode gap discovery groupings were concordant with species clades. There were 41 sequence operational taxonomic units (31 clusters and 10 singletons) recognised with automatic barcode gap discovery. All operational taxonomic units were generated by the initial partition with a prior intraspecific divergence of Pmax = 0.00774 (Table 1). Most groups represented species clades (24%), including the two L. ferruginea operational taxonomic units, and clades within species clades (71%). The remaining 5% were the undetermined Leuctra operational taxonomic unit and the outgroup.
There were low (≤ ca. 3%) minimum interspecific divergence values for regional members of the L. ferruginea group, and these were only most notably between L. truncata and the two L. ferruginea operational taxonomic units. The low-divergence value between the two L. ferruginea operational taxonomic units is concordant with the maximum-likelihood tree (Fig. 6). These two operational taxonomic units were nested as a single clade, albeit without strong nodal support (< 80%). The Bayesian-inference tree grouped both L. ferruginea operational taxonomic units with L. truncata as a three-taxon polytomy, also without strong nodal support (< 0.80). Overall, the general lack of resolution depicted within and between species-group relationships was a common theme with both the maximum-likelihood and Bayesian-inference analyses. Most notable was paraphyly shown for the L. tenuis group (Figs. 5–6). Whereas L. schusteri and L. tenuis appear to form a sister-species group with strong nodal support (maximum-likelihood: 92%; Bayesian-inference: 1.00), the remaining species (L. carolinensis, L. tenella, L. triloba, and L. variabilis) form their own monophyletic group (maximum-likelihood: 82%; Bayesian-inference: 1.00). Despite the possible paraphyletic nature of the L. tenuis group (Figs. 5–6), this group as a whole was most closely grouped to members of the L. ferruginea group. Minimum interspecific divergence values ranged from 5.4% to 9.0% (Table 1).
Discussion
General patterns
For both the Gattolliat et al. (Reference Gattolliat, Vinçon, Wyler, Pawlowski and Sartori2016) and Vitecek et al. (Reference Vitecek, Vinçon, Graf and Pauls2017) phylogenetic treatments of Palearctic Leuctra, the molecular approach generally supported morphology-based species groups (i.e., Harper and Harper Reference Harper and Harper2003). In this study, mtCO1 phylogenetic analyses showed strong nodal support for 18 monophyletic groups that represented at least 16 of the 31 Nearctic species. Males were grouped into distinct clades with high to very high nodal support. These values likewise provided strong support for inclusion of sequenced females, thus permitting confidence with mtCO1-based species determinations in lieu of subjective, subtle morphological differences.
Females of most Nearctic Leuctra species are difficult to identify to species with confidence. Notable exceptions include L. duplicata and L. maria Hanson, 1941 (Grubbs and Wei Reference Grubbs and Wei2017), L. hicksi Harrison and Stark, 2010 (Grubbs et al. Reference Grubbs, Metzger and Wei2020), and L. grandis. We found that the shape of the subgenital plate lobes and medial cleft for females studied here appear different with small changes in the viewing angle (e.g., Fig. 3E–F). Plates can appear different with age (e.g., teneral versus mature; Fig. 3A–B), and the space between lobes also varies between gravid and nongravid females. Using morphological features alone to identify females would have resulted in many incorrect determinations for all species except L. grandis. An integrative taxonomic approach has provided a more accurate and robust data set.
Many Leuctra species in our study showed a high amount of maximum intraspecific genetic variation. Broad ranges have been reported for several species of Plecoptera (Zhou et al. Reference Zhou, Adamowicz, Jacobus, DeWalt and Hebert2009; Mynott et al. Reference Mynott, Webb and Suter2011; Sweeney et al. Reference Sweeney, Battle, Jackson and Dapkey2011; Boumans and Baumann Reference Boumans and Baumann2012; Avelino-Capistrano et al. Reference Avelino-Capistrano, Nessimian, Santos-Mallet and Takiya2014; Gill et al. Reference Gill, Kondratieff, Stark and Sandberg2015a). For example, Avelino-Capistrano et al. (Reference Avelino-Capistrano, Nessimian, Santos-Mallet and Takiya2014) reported maximum intraspecific values ranging from 0% to 15% for Brazilian Kempnyia Klapálek, 1914. Mynott et al. (Reference Mynott, Webb and Suter2011) reported maximum intraspecific values ranging from 0% to 5.8% for Australian Riekoperla McLellan, 1971. Sweeney et al. (Reference Sweeney, Battle, Jackson and Dapkey2011) reported high intraspecific variation for Nearctic Perlesta Banks, 1906 but with no gaps between intraspecific and interspecific pairwise values.
The bi- and tri-model distributions reported here for intraspecific data are not due solely to a large proportion of individuals. Leuctra grandis, L. ferruginea OTU 1, L. sibleyi, and L. triloba each had maximum intraspecific values that were greater than 5%. High values for L. ferruginea OTU 1 and L. triloba, however, were both due to only one divergent individual. Leuctra sibleyi exhibited a trimodal divergence distribution, with one individual yielding the highest intraspecific variation value (18.6%) and three individuals representing an intermediate haplotype (10%). Males of L. sibleyi were morphologically inseparable even though the intraspecific variation was high. Considering the small geographic scope of this study, there are essentially no barriers separating these groups, and intermediate haplotypes may not have been collected. Leuctra grandis exhibited a bimodal divergence distribution, with 16 individuals showing maximum intraspecific divergence values that were greater than 5% (Fig. 8). Again, intermediate haplotypes may not have been collected.
There was overlap in distance values between individuals of species with high maximum intraspecific and low minimum interspecific values. Most overlaps were due to low numbers of divergent individuals, as noted above. For example, if those few individuals were removed, the following species or operational taxonomic unit pairs would no longer overlap – L. ferruginea OTU 1 and L. ferruginea OTU 2, and L. variabilis and L. triloba. The only species pairs that would still have overlapping values would be L. ferruginea OTU 2 and L. truncata, and L. ferruginea OTU 1 and L. truncata. These latter two species have been placed in the L. ferruginea group (Harper and Harper Reference Harper and Harper1997).
The overlap between intraspecific and interspecific genetic distance values suggests that other methods may better explain the data. The automatic barcode gap discovery algorithm separated all species into singletons or clusters. There were no shared clusters between morphologically defined species. In many cases, automatic barcode gap discovery overdivided species into multiple clusters or singletons due to variations in intraspecific values (e.g., L. sibleyi). The variation could be due to uncommon variant haplotypes, diverse population genetic structure, incomplete geographic sampling (i.e., intermediate elevation sites were most difficult to access), or a combination of the three.
Species groups
Members of all five Nearctic species groups (Harper and Harper Reference Harper and Harper1997) were included in this study, although none were fully represented. Group membership, however, was only partially supported by the mtCO1 phylogenetic analyses. Including L. moha, which remains ungrouped (Harper and Harper Reference Harper and Harper1997; Grubbs et al. Reference Grubbs, Metzger and Wei2020), the maximum-likelihood and Bayesian-inference phylogenetic trees had similar topology. The maximum-likelihood analyses recognised all individuals from Mount Mitchell as monophyletic units, including nine recognised species. The Bayesian-inference tree similarly showed similar patterns of monophyly, except for the polytomy of L. ferruginea OTU 1, L. ferruginea OTU 2, and L. truncata. The deeper relationships (e.g., species groups), however, were not well resolved in either the Bayesian-inference or maximum-likelihood analysis.
Leuctra alexanderi and L. mitchellensis (L. biloba group) consistently grouped together but only with L. atla (L. grandis group). Harper and Harper (Reference Harper and Harper1997) noted an apparent close relationship between the L. biloba and L. grandis groups. Interestingly, in the present study, L. duplicata always grouped with the three aforementioned species with strong nodal support. The inclusion of more sequences of L. duplicata, together with L. maria – the two species of the L. duplicata group – in subsequent analyses is needed to better confirm this relationship. Although the original inclusion of L. maria was tentative (Harper and Harper Reference Harper and Harper1997), morphological similarities between these two species in paraproct, vesicle, and subgenital plate features (Grubbs and Wei Reference Grubbs and Wei2017) suggest a sister-species relationship.
Leuctra alta was placed tentatively in the L. grandis group (Harper and Harper Reference Harper and Harper1997) due to common characteristics with other members (i.e., bilobed dorsal process, narrow subanal lobes). In the present study, the L. grandis group was represented by three of five species, but none grouped together. Despite the superficial similarities in the shape of the dorsal process of L. alta with that of L. sibleyi and in the paraproct characteristics of L. grandis with those of L. sibleyi (Grubbs, unpublished data), from a phylogenetic perspective, these three species were always recovered as completely independent of each other. Moreover, minimum interspecific divergence values ranged from 11.7% to 14.4%.
The L. tenuis group was paraphyletic. Harper and Harper (Reference Harper and Harper1997) split this group into two subgroups, based on emergence timing and subanal lobe characteristics. They placed L. carolinensis and L. tenella into one subgroup and the remaining species into a second. Although two monophyletic units were recovered in the present study, they do not support Harper and Harper (Reference Harper and Harper1997). Monophyly of L. carolinensis and L. tenella was not recovered in either analysis. In the present study, L. schusteri and L. tenuis were clustered as a sister-species pair with high nodal values. This was further supported by relatively low interspecific divergence (5.4%). This agrees with Grubbs (Reference Grubbs2015), who noted that these two species vary from the other members of the group by dorsal-process characteristics. Of additional importance was that the sequences of L. tenuis from Kentucky and Tennessee, United States of America, and from Nova Scotia, Canada, grouped together with the North Carolina specimens with good to strong nodal support. The remaining four group species (L. carolinensis, L. tenella, L. triloba, and L. variabilis) grouped together with good to strong nodal support. Similarly, interspecific divergence only ranged from 4.5% to 8.1%.
Although two subgroups were proposed for the L. ferruginea group (Harper and Harper Reference Harper and Harper1997), this taxon is an enigma. For example, Harrison and Stark (Reference Harrison and Stark2010) provided good evidence using scanning electron microscopy that L. rickeri is a junior synonym of L. alabama James, 1974. This is particularly relevant because these two species were placed into different subgroups. Although the Kentucky specimens of L. rickeri that were used in the present study grouped with L. ferruginea and L. truncata, the single sequence of L. rickeri from Illinois, United States of America, was consistently basal to the clade of L. ferruginea group + L. tenuis group + L. moha.
Unresolved species
There were two separate L. ferruginea operational taxonomic units and the undetermined Leuctra operational taxonomic unit. The low interspecific divergence value (3.3%) between the L. ferruginea operational taxonomic units and the close grouping on both the Bayesian-inference and maximum-likelihood trees suggest that they represent one distinct, albeit genetically variable, species. An important item to note is that the single BOLD sequence for a male from Ontario, Canada consistently grouped with L. ferruginea OTU 1 in all phylogenetic approaches.
Although there were only three undetermined Leuctra operational taxonomic unit females, the high minimum interspecific divergence values suggest that this may represent a distinct species. What are the options based on currently recognised species? No Mount Mitchell females grouped with the L. carolinensis individuals (Figs. 5−6). This is a regional species with a “North Carolina, Black Mountains” type locality (Claassen Reference Claassen1923) that was collected in east Tennessee, only approximately 40 km from Mount Mitchell but not during the present study. Leuctra biloba Claassen, 1923, L. nephophila Hanson, 1941, and L. monticola Hanson, 1941 were the only three remaining regional Nearctic Leuctra species not included in this analysis. Leuctra biloba is endemic to the southern Appalachian Highland region but was not collected in our study. Leuctra nephophila is also endemic to the southern Appalachian Highland region. Although not included here, we included BOLD sequences (SMSTO172, SMSTO173) of L. nephophila in preliminary analyses. This species did not group with the Leuctra operational taxonomic unit females, and interspecific divergence values exceeded 10%. Leuctra monticola is arguably the rarest Nearctic Leuctra, known only from the type series collected from Cades Cove (Great Smoky Mountain National Park, Tennessee) in 1938. Overall, males and subsequent molecular analyses are needed before a species name can be confidently applied to the Leuctra operational taxonomic unit individuals.
Conclusions
An integrative morphological and molecular approach added robust taxonomic information to the Leuctra faunal diversity of Mount Mitchell, North Carolina. Species names were assigned with confidence to an additional 235 females. The additional information included from other 2019 localities and from BOLD sequences served to better contextualise tree- and genetic distanced–based analyses. Morphological characteristics and gene sequence analyses may not align the same species into individual groups. Evidence presented in the present study suggests revisions to group memberships are plausible with future phylogenetic work. No changes, however, are proposed at this time to the Nearctic Leuctra species groups (Harper and Harper Reference Harper and Harper1997) due to the use of only one gene. Previous studies have shown that the use of a single mitochondrial gene is often insufficient to understand relationships amongst different species (Mynott et al. Reference Mynott, Webb and Suter2011; Gill et al. Reference Gill, Kondratieff, Stark and Sandberg2015a). Additional genes (e.g., nuclear 28S) or single nucleotide polymorphisms of whole genomes could help resolve these older relationships. Vitecek et al. (Reference Vitecek, Vinçon, Graf and Pauls2017) used a concatenated data set of two mitochondrial and two nuclear genes to assess the L. inermis group, but even a combination of faster- and slower-evolving genes did not resolve relationships between some species. The authors suggested that diversification within the species group was recent. A similar phenomenon could be operating on Nearctic Leuctra. Molecular characterisation of the entire genus Leuctra from both Nearctic and Palearctic realms should be considered as more DNA sequence data become available.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.4039/tce.2022.5.
Acknowledgements
Funding to support this research was provided by the Western Kentucky University (WKU) Center for Biodiversity Studies and the WKU Graduate School. Cris Metzger, Elizabeth Metzger, Margaret Metzger, and Garret Kratina kindly assisted with long hours of fieldwork. Garret Kratina also provided assistance with ArcMap. Three WKU faculty, John Andersland, Kevin Bilyk, and Albert Meier, helped with microscopy, workstation issues, and research ideas and editing, respectively. Dr. R. Edward DeWalt (Illinois Natural History Survey, Champaign-Urbana, Illinois) gave helpful advice for successful DNA extractions from stoneflies.
Competing interests
The authors declare no competing interests.
