Introgression and species demarcation in western European Leuctra fusca (Linnaeus, 1758) and L. digitata Kempny, 1899 (Plecoptera: Leuctridae)

ABSTRACT The compilation of a DNA barcoding library of Norwegian stonefly (Plecoptera) species revealed that Leuctra fusca (Linnaeus, 1758) and Leuctra digitata Kempny, 1899 (Leuctridae) share haplotypes in northernmost Scandinavia. Phylogenetic analyses of the mitochondrial (mt) DNA barcode marker COI and the nuclear marker 28S show that the shared haplotypes must result from the introgression of a L. fusca mitochondrion into a L. digitata population on at least two occasions. Although mt introgression is widespread in animals, this represents the first documented case in Plecoptera. This study also included specimens of L. cf. fusca from the Sierra Nevada massif in Spain, a population previously known as L. carpentieri Despax, 1945. Their mt haplotypes are ca. 13% different from other European L. fusca. However, their 28S alleles are compatible with their morphological identification as L. fusca. In view of the possibility of mt introgression, the taxonomic status of this population remains undecided.


Introduction
Inventories of macroinvertebrate fauna are routinely used in the assessment of freshwater ecosystem quality. In the near future, morphological identifications will likely become complemented with assay platforms based on DNA signatures that will increase the precision and ease of quantifying biodiversity (Pfrender et al. 2010;Hajibabaei, Shokralla, Zhou, Singer, and Baird 2011;Sweeney, Battle, Jackson, and Dapkey 2011). For animals, cytochrome c oxidase subunit I (COI), also known as the barcode gene, is the most widely used marker for species identification (Hebert, Cywinska, and Ball 2003). During the assembly of a DNA barcode library of Norwegian stoneflies (project 'Norwegian Barcoding of Life [NorBOL] À Freshwater Insects'), it was revealed that Leuctra digitata Kempny, 1899 specimens from the northernmost province of Finnmark share haploclades with Leuctra fusca (Linnaeus, 1758), a closely related species belonging to the same species-group (DeWalt, Maehr, Neu-Becker, and Stueber 2015), suggesting mitochondrial (mt) introgression. Both are autumn-emerging species that often occur syntopically in Scandinavia (records in the Norwegian database Artskart (http://artskart.artsdatabanken. no/); in central Europe they are sympatric at the water basin level (Sold an, Zahr adkov a, Hele sic, Du sek, and Landa 1998, pp. 281, 284, 288) or even in the same section of a stream (Mur anyi 2006, p. 88). Adults of these species can be distinguished in both sexes by means of the genital morphology (Lillehammer 1988, pp. 141À143;Tierno de Figueroa, S anchez-Ortega, Membiela Iglesia, and Luz on-Ortega 2003, pp. 294À295, 300À301).
In order to distinguish introgression from recent hybridisation, we verified our morphological identification of Norwegian specimens by sequencing the nuclear marker 28S. We examine the geographical extent of this mt introgression and evaluate the alternative explanation of incomplete lineage sorting (ILS) with the help of additional specimens and DNA barcode sequences of both species from Finnmark, Finland and central and southern Europe.
Our data-set also includes brachypterous and normal-winged L. cf. fusca specimens from the Sierra Nevada mountain range in southern Spain. Despax (1945) originally described this population as a distinct species, Leuctra carpentieri Despax, 1945, from brachypterous individuals from high mountain localities. Aubert (1963) lowered this population to a subspecies of L. fusca indicating that individuals were macropterous until 1800 m a.s.l. S anchez-Ortega and Alba-Tercedor (1985), finding no defining morphological characters other than the brachyptery found in some, but not all, individuals from high altitude, synonymised the taxa. The taxonomic position of the Sierra Nevada population is re-evaluated in the light of mt and nuclear sequence data.

Specimens
Leuctra Stephens, 1836 specimens from Finnmark were collected in the framework of large-scale inventory of aquatic insects in the Norwegian county of Finnmark, using Malaise traps and manual collecting methods (Boumans and Brittain 2012;Ekrem et al. 2012).
Further specimens were collected in southern and central Norway as part of the DNA barcoding project 'NorBOL À Freshwater Insects' (Norwegian Ephemeroptera, Plecoptera and Trichoptera, NOEPT) and supplemented with specimens from central and southern Europe (barcoding project West Palaearctic Plecoptera, WPPLE). We included six specimens of L. cf. fusca from the Sierra Nevada mountain range in southern Spain, two of which are brachypterous. Adult specimens were identified on the basis of tergites VI and VII in the males, and sternite VIII in the females.

Sequences
Part of the COI sequences discussed in this paper were produced at the sequencing facility of the Canadian Centre for DNA Barcoding in Guelph; for another part, the sequence amplicons were produced at the Natural History Museum in Oslo using the primers LCO1490-L and HCO2198-L (Nelson, Wallman, and Dowton 2007) and the protocol described in Boumans and Baumann (2012), while the final sequencing was outsourced to Macrogen Europe in Amsterdam.
Further COI sequences were downloaded from the Barcode of Life Data System (BOLD) database and added to our data matrix, notably one sequence from Latium, Italy (Fochetti, Sezzi, Tierno de Figueroa, Modica, and Oliverio 2009), six sequences from the Barcoding Fauna Bavarica project (Hendrich et al. 2010), three from the Swiss DNA barcoding initiative (Gattolliat et al. 2015; http://www.swissbol.ch) and 29 sequences originating from one Malaise trap site in central Finland, belonging to the Global Malaise Trap program (http://globalmalaise.org). We borrowed these Finnish specimens from the Biodiversity Institute of Ontario, University of Guelph and identified them morphologically. A sequence of Leuctra hippopus Kempny, 1899 was included as outgroup.
In order to verify our morphological identifications of L. fusca and L. digitata, we amplified the D2 region of nuclear ribosomal gene 28S in Oslo with the primers Road1a and Road4b (Crandall, Harris, and Fetzner 2000;Whiting 2002), following the protocol described in Boumans and Baumann (2012); these amplicons were likewise sequenced at Macrogen.
All sequences were uploaded to the database of the BOLD (Ratnasingham and Hebert 2007), projects NOEPT and WPPLE. All specimens and sequences are listed in Appendix 1 (see Supplemental Material online). Additional metadata on these specimens as well as photographs are available online from the BOLD database.

Sequence analyses
Phylogenetic trees of the COI sequences were constructed with distance (neighbour joining, NJ) and maximal parsimony (MP) methods implemented in Paup Ã , and with Bayesian inference (BI) in MrBayes version 3.2.2 (Ronquist and Huelsenbeck 2003). The distance measure for the neighbour-joining method in Paup Ã was set as 'dist D GTR rates D equal Pinvar D 0.6972', based on the evolution model selected under the Akaike information criterion with MrModeltest 2.2. (Nylander 2004). We carried out heuristic searches under the optimality criteria distance and parsimony with tree bisec-tionÀreconnection branch swapping and 100 random addition sequence replicates. Bootstrapping (2000 replicates) was performed to obtain support values for branches. For BI, the COI data were divided in a partition for the third codon position, and another for the first and second position. Based on the Akaike information criterion, the evolution model Generalised time-reversible with a proportion of invariable sites (GTR C I) was selected for the first and second codon position and Generalised time-reversible with gammadistributed rate variation (GTR C G) for the third position. We ran two independent analyses, each consisting of four Markov chains that ran for 40 £ 10 6 generations and were sampled every 1000 generations, with default priors and partition-specific substitution rates (setting prset ratepr D variable). After discarding the first 10 million generations, remaining trees from both analyses were combined and a 50% majority rule consensus tree was calculated. MrBayes output files and Tracer v1.6.0 (Drummond and Rambaut 2007) were used to inspect trace plots and convergence diagnostics (average standard deviation of split frequencies < 0.01, effective sample size > 200) in order to ensure that the Markov chains had reached statistical stationarity and converged on the parameter estimates and tree topology after the burn-in phase, which was set at 25%. Calculations in MrModeltest, Paup Ã and MrBayes were performed at the Lifeportal computer facility (www.lifeportal.uio.no) at the University of Oslo, Norway.
The 28S sequence divergence within and between both species was visualised with the median-joining network method (Bandelt, Forster, and R€ ohl 1999) implemented in the program Network 4.613 (www.fluxus-engineering.com).

Species identification by means of morphology and 28S
The partial 28S nuclear marker (783 bp) distinguishes Norwegian L. fusca and L. digitata, concurring with morphological identification in all studied specimens (Figures 1, 3). The 28S alleles of L. digitata and (L. fusca C L. cf. fusca) are separated by 4À6 substitutions (Figure 2). Specimens from Norway showed no intraspecific variation. L. fusca specimens from Spain show intraspecific variation, with some specimens bearing the allele found in Norway, and others showing up to three substitutions relative to this allele. Three different alleles were found in the Sierra Nevada.
Twenty-three of the specimens from central Finland, borrowed from the Biodiversity Institute of Ontario, were identified morphologically as L. digitata ( Figure 1GÀH; additional photographs in the BOLD database). Two female specimens had an atypically shaped subgenital plate, with a deep incurvature between the lateral lobes of sternite 8 ( Figure 1H). Typical L. digitata has a straight posterior border between slightly outward pointing lateral lobes ( Figure 1F (Figure 174 in Aubert 1959, p. 58). Six further DNA-barcoded Finnish specimens from the same Malaise trap were not identified because they were incomplete or damaged.

Mitochondrial phylogeny
Phylogenetic tree reconstruction of COI with BI, MP and NJ yields largely the same topology. Figure 3 shows the tree diagram resulting from BI with support values from the three methods; incompatible topologies all involve relationships without statistical support and are indicated below. Four major groups are retrieved: I. Sierra Nevada, a well-supported haploclade. II. L. digitata specimens from Central Europe, southern Norway and central Finland. These sequences form a supported haploclade in the MP and NJ trees (brace in Figure 3). The BI tree diagram, however, suggests this L. digitata cluster is paraphyletic with respect to the next two clades, but statistical support for this paraphyly is negligible (posterior probability 0.60). III. The single sequence from Latium, Italy. IV. A well-supported haploclade that includes all other L. fusca from southern and northern Europe plus L. digitata from northernmost Norway and part of the L. digitata samples from central Finland. The relationships between these four groups are not resolved. Notably, the specimen from Latium is shown as sister to clade IV in the Bayesian tree and as sister to clade II in the MP and NJ trees, but without relevant statistical support in any of the analyses. Clades II and IV are 6.7%À7.7% different (uncorrected distance). Two well-supported subclades of group IV include all L. fusca and L. digitata specimens from northern Norway and half of the L. digitata from Finland, with identical and near-identical (1 substitution) haplotypes shared by both species. Groups II and IV each include one of the two Finnish L. digitata females that have an atypical subgenital plate.
Haplotypes of cf. L. fusca from the Sierra Nevada are 12.9%À13.7% different from those of unambiguous L. fusca. Neither COI nor 28S distinguishes brachypterous from normal-winged Sierra Nevada specimens (Figures 2, 3). The single sequence from Italy is ca. 10% different from the sequences of group IV.

Hybridisation and introgression in L. digitata
Since two subclades of L. fusca include sequences of L. digitata, the COI phylogenetic tree is clearly not congruent with the species identification based on morphology and the nuclear marker 28S. Non-congruence of mitochondrial and nuclear gene trees can be due to either introgression or ILS.
As for ILS, the fact that ancestral polymorphism sorts randomly into descendent lineages can lead to two ancestral pre-speciation mt lineages occurring in one of the descendant species and only one in another descendant species. However, when two species share (near) identical haplotypes, ILS can only explain this finding if speciation occurred so recently that there has not been time for genetic divergence (McGuire et al. 2007). In the case of L. digitata, identical or near-identical haplotypes are shared with L. fusca in northern Fennoscandia, whereas uncorrected COI sequence divergence is ca. 7% elsewhere. The shared mt subclades are embedded in a larger, diverse fusca haploclade ( Figure 3) and must be younger than the diversification of the two species. Furthermore, if the shared haplotypes and subclades were already present in the common ancestor of both species, we expect them to be randomly distributed in the descendant populations of both species (McGuire et al. 2007 and references therein). L. digitata with L. fusca-type mitochondrion have thus far only been encountered in northern Fennoscandia. Here, the species occur sympatrically at many sites and the populations must be rather young, postdating the end of the last glacial period (10,000 years BP). In sum, based on the available evidence, L. digitata represents a typical example of recent mt introgression.
Hybridisation and introgression of the maternally inherited chloroplast or mitochondrial DNA is common in plants and animals respectively (e.g., Chan and Levin 2005;Linnen and Farrell 2007;Kindler et al. 2013). Introgression is often asymmetric, with all hybrid individuals bearing the cytoplasmic DNA of one species. A plausible explanation is frequency-dependant assortative mating, as when a (choosy) female of the rare species accepts mating with a heterospecific male (Chan and Levin 2005). So far, we have not encountered L. fusca specimens with a L. digitata mitochondrion. However, at sampling sites where both species occur, either species can be preponderant, as we found in a study of Malaise trap samples from Finnmark (Boumans and Brittain 2012).
Another hypothesis explaining the asymmetric mt introgression and the preponderance of the introgressed mitochondrion in one region is selective advantage associated with the L. fusca haplotypes. The L. fusca mitochondrion may be more adapted to the climatic conditions of the high North. Alternatively, a selective sweep on the mt genome can be caused by maternally inherited symbionts like Wolbachia that bestow a selective advantage on infected females either by causing cytoplasmic incompatibility (Ballard 2000) or killing infected male offspring (Johnstone and Hurst 1996;Graham and Wilson 2012). In other insects, strains of symbionts have been found associated with mt haploclades (e.g., Kvie, Hogner, Aarvik, Lifjeld, and Johnsen 2012).
Finally, the asymmetric introgression may reflect an asymmetry in pre-zygotic or postzygotic reproductive barriers. For example, L. digitata females may less readily accept heterospecific mating, or F1 hybrids with a L. digitata mother suffer more fitness loss. At present, nothing is known about the reproductive barriers between L. fusca and L. digitata. What is known about their biology is very similar: both are rheophilic detritus feeders with adult emergence in late summer and autumn. Habitat preferences are similar, although L. fusca is more eurytopic (Graf et al. 2009), and they often occur in sympatry. Vibrational sexual communication can constitute a pre-mating reproductive barrier in stoneflies (Boumans and Johnsen 2014). However, mating signals are not known from either species, though attested in some other Leuctra species (Rupprecht 1977;Boumans and Johnsen 2015, Suppl. 1).
Introgressed mitochondria were found in all nine L. digitata specimens tested from Finnmark (69 À70 N), half of the specimens from central Finland (65 N) and none of the specimens from more southern sites in Norway or Europe (south of 62 N). This geographical pattern may result from either of two different processes. First, it may simply reflect the dispersal history of the species, as when northernmost Scandinavia was accidentally populated by an introgressed population from the northeast, whereas the remainder of Scandinavia was colonised from the South. Mitochondrial phylogenies of other stonefly species indicate that these colonised the Scandinavian Peninsula from both sides (Boumans and Baumann 2012;Boumans and Brittain 2012). However, our data do not show in which direction the introgressed population dispersed. The introgressed population may also have arisen in the upper North and dispersed southward into central Finland. Second, the geographical distribution can result from a selective sweep replacing the ancestral mitochondrion with the introgressed one. Testing these alternative hypotheses would require additional sampling, extensive genotyping and testing for symbiont infection along NorthÀSouth transects in Scandinavia and Finland/Karelia. Genetic studies are also necessary to reveal whether the atypical morphology of the subgenital plate found in some Finnish L. digitata females (Figure 1h) results from hybridisation with L. fusca.

The populations from the Sierra Nevada and Latium
The very divergent COI sequences of L. cf. fusca from the Sierra Nevada revive the question whether this population belongs to a different species, in which case the name L. carpentieri should be resurrected. However, introgression can lead to similarly distinct haplotypes occurring within a single species, even at a single site (cf. the Finnish Malaise trap site). The 28S alleles do support the identification as L. fusca. For these reasons the taxonomic assessment of these L. cf. fusca populations requires further research, including sampling of L. fusca from southern Spain outside the Sierra Nevada and study of further nuclear or morphological markers. Likewise, the Italian COI sequence originating from Fochetti et al.'s (2009) study may represent a different species. Interestingly, the BOLD database includes a sequence from Iraqi Kurdistan that is very similar to this sequence. Based on its sampling site, this Iraqi specimen could be L. fusca latior Berthelemy and Dia, 1982. Finally, further examination of other L. fusca populations, particularly from the Far Eastern Palaearctic, as well as closely related species like L. tergostyla Wu, 1973 from China, would also contribute to the taxonomic assessment of these taxa.

Practical implications
Mt introgression introduces error in DNA barcode-based species inventories (e.g., Hurst and Jiggins 2005). A practical consequence of introgression is that DNA barcoding by means of COI cannot reliably distinguish between L. fusca and L. digitata in Fennoscandia, at least when the attested haplotype belongs to the fusca clade. At the same time, this case may raise awareness of possible mt introgression in other geographical regions and in other aquatic taxa.

Disclosure statement
The authors are not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as constituting a conflict of interest.

Funding
The Natural History Museum in Oslo (NHMO) financed lab work and other research facilities.