Taxonomist's Nightmare … Evolutionist's Delight : An Integrative Approach Resolves Species Limits in Jumping Bristletails Despite Widespread Hybridization and Parthenogenesis.

Accurate species delimitation is fundamental to biology. Traditionally, species were delimited based on morphological characters, sometimes leading to taxonomic uncertainty in morphologically conserved taxa. Recently, multiple taxonomically challenging cases have benefited from integrative taxonomy-an approach that highlights congruence among different disciplines and invokes evolutionary explanations for incongruence, acknowledging that different methods can mirror different stages of the speciation continuum. Here, we used a cohesive protocol for integrative taxonomy to revise species limits in 20 nominal species and 4 morphospecies of an ancestrally wingless insect group, the jumping bristletail genus Machilis from the European Eastern Alps. Even though morphologically conserved, several small-scale endemic species have been described from the Eastern Alps based on variation in hypodermal pigmentation patterns-a highly questionable character. As valuable as these endemics are for conservation, they have never been verified by alternative methods. Using traditional morphometrics, mitochondrial DNA, ribosomal DNA, and amplified fragment-length polymorphism markers, we identify six nominal species as taxonomic junior synonyms (Machilis alpicola Janetschek, 1953 syn. n. under M. vagans Wygodzinsky, 1941; M. ladensis Janetschek, 1950 syn. n., M. robusta Wygodzinsky, 1941 syn. n., and M. vicina Wygodzinsky, 1941 syn. n. under M. inermis Wygodzinsky, 1941; M. aleamaculata Wygodzinsky, 1941 syn. n. under M. montana Wygodzinsky, 1941; M. pulchra Janetschek, 1950 syn. n. under M. helleri Verhoeff, 1910) and describe two new species (Machilis cryptoglacialis sp. n. and Machilis albida sp. n.), one uncovered from morphological crypsis and one never sampled before. Building on numerous cases of incongruence among data sources, we further shed light on complex evolutionary histories including hybrid speciation, historical and recent hybridization, and ongoing speciation. We hypothesize that an inherent affinity to hybridization, combined with parallel switches to parthenogenesis and repeated postglacial colonization events may have boosted endemicity in Eastern Alpine Machilis We thus emphasize the importance of integrative taxonomy for rigorous species delimitation and its implication for evolutionary research and conservation in taxonomically challenging taxa.

Recently, species delimitation has experienced a renaissance among systematic and evolutionary biologists. With the unified species concept (De Queiroz 2007), a promising synthesis has been achieved which relaxes the rivalry among different species concepts and fosters multidisciplinary reasoning. Methodological advances regarding DNA-based species delimitation (Pons et al. 2006; Knowles and Carstens 2007; Yang and Rannala 2010; Zhang et al. 2013) and increasing popularity of an integrative approach to taxonomic problems (Dayrat 2005; Will et al. 2005; Padial et al. 2010; Schlick-Steiner et al. 2010) have triggered a growing number of comprehensive species delimitation studies (animal examples: Raxworthy et al. 2007; Leaché et al. 2009; Glaw et al. 2010; Ross et al. 2010; Pinzón and LaJeunesse 2011; Seppä et al. 2011; Gebiola et al. 2012; Puillandre et al. 2012; Satler et al. 2013; Sistrom et al. 2013; Derkarabetian and Hedin 2014; Wachter et al. 2015). By combining results from analyses covering independent evolutionary aspects of the same set of specimens, integrative taxonomy can provide robust species hypotheses even in taxonomically intricate instances and uncover complex evolutionary histories. Unfortunately, these strengths are often impaired by the reluctance of authors or editors to publish taxonomic consequences such as formally describing new species or synonymizing species resulting from species delimitation studies (Pante et al. 2014). This reluctance has presumably been based, among other reasons, on the common belief that taxonomic studies are poorly cited. However, Steiner et al. (2015) found no evidence for a citation impediment in the taxonomic literature.

Accurate species delimitation remains a challenging task in morphologically poorly differentiated animal taxa (e.g., Satler et al. 2013; Andújar et al. 2014; Wachter et al. 2015) or in cases where only larval specimens are available (e.g., Zhou et al. 2007), where appropriate morphological and molecular characters to quantify variation among proposed species hypotheses are lacking. Often, mitochondrial genes are the first markers of choice due to their fast rate of evolution and ease of amplification. However, mitochondrial phylogenies are especially prone to deviations from the species tree due to hybridization, introgression, incomplete sorting of ancestral polymorphism, and infection with reproduction-manipulating endosymbionts (Funk and Omland 2003). Moreover, nuclear pseudogenes can distort supposedly mitochondrial gene trees and are difficult to identify (Song et al. 2008). Thus, a modern species delimitation approach ideally includes multiple nuclear loci. Inference of allele variants in nuclear sequences is a minor issue in diploids, but higher ploidy levels pose a serious challenge to amplification, sequencing, and even more to the phylogenetic analysis of orthologous sequences. As a consequence, amplified fragment-length polymorphisms (AFLPs) have been frequently used to characterize nuclear genomic variation among species when polyploids are involved (Hedrén et al. 2001; Guo et al. 2005; Koopman et al. 2008; Meudt et al. 2009). These and other studies (Dasmahapatra et al. 2009; Fink et al. 2010; Khan et al. 2014; Kirchberger et al. 2014) have shown that AFLPs are suitable to reconstruct complex evolutionary histories when frequent hybridization, incomplete lineage sorting, or polyploidization complicate the use of nuclear markers.

In this study, we apply an integrative taxonomic approach to survey species diversity in the jumping bristletail genus Machilis, which is poor in morphological characters suitable for species delimitation. Jumping bristletails (Archaeognatha or Microcoryphia) are ancestrally wingless insects that present many morphological features plesiomorphic for Hexapoda (Sturm and Machida 2001) and are thus often called “primitive”. Phylogenetically, they are sister to the Dicondylia, the group comprising silverfish and winged insects (Hennig 1953; Misof et al. 2014). Worldwide, about 450 species have been described from ecosystems as different as tropical rainforests, marine supralittorals, dry steppes, and glaciated mountains. At 94 species, the genus Machilis (Mendes and de Jong 2013) is the most species-rich group of jumping bristletails, and at 55 species, its diversity is highest in the European Alps; “Alps” and “Alpine” henceforth refer to this mountain system. The radical drop in Machilis species numbers in European lowlands might be explained by their petrophilous lifestyle (Sturm and Machida 2001) rather than collection bias. As far as known, most species are confined to either the Western (32 spp.) or the Eastern Alps (21 spp.). Just two Machilis species occur across this geological (Grimm and Mattmüller 2004) and biogeographical (Schönswetter et al. 2005; Schmitt 2009) border, the Rhine–Splügen line.

Due to their winglessness and low vagility, jumping bristletails are thought to be slow dispersers. Accordingly, the genus Machilis includes a high portion of small-scale endemics compared with other arthropod groups in the Eastern Alps (Rabitsch and Essl 2009). Eleven (48%) of the twenty-three Machilis species known from the Eastern Alps have exclusively been recorded from their type locality and are thus classified as local endemics sensu Cowling (2001; Table 1). Another six species (26%) have a known distribution range below 10,000 km2 and are termed narrow-range endemics sensu Harvey (2002). The genus Machilis thus represents an outstanding example of Alpine biodiversity. However, many of the original species descriptions were based on just a few individuals. Moreover, due to the absence of prominent morphological characters, many species were delimited on the basis of hypodermal pigmentation patterns, which have also been reported to vary considerably within some Machilis species (Janetschek 1955). As a consequence, the real number of species may differ from the number of described species, and thus also the number of true local and narrow-range endemics remains unclear. Accurate species delimitation is therefore key to making this group accessible to evolutionary research, biodiversity studies, and conservation measures. A classical taxonomic revision of these species is, however, hampered by the loss or poor condition of type material (detailed information compiled in Table 1) and by the absence of additional samples in museums or universities—a consequence of the longstanding lack of taxonomic expertise in this animal group.

Table 1.

The Machilis species of the Eastern Alps

Genus	Species	Taxon authority	Endemic status	Geographic distribution	Reproductive mode	Number of sampled individuals	Information on type specimens
Machilis	aciliata	Janetschek, 1955	L	EA	Sexual	—	Holotype not at LMFI, untraceable (this study).
Machilis	aleamaculata	Wygodzinsky, 1941a	L	EA	Parthenogenetic	7	Holotype at NMB; uninformative due to dissolved tissue.
Machilis	alpicola	Janetschek, 1953	A	WA + EA	Parthenogenetic	19	Unclear if any type was deposited at MNHNP (no response).
Machilis	alpina	Riezler, 1941	L	EA	Parthenogenetic	—	No type specified by Riezler (1941); synonymized under M. lehnhoferi by Janetschek (1955).
Machilis	anderlani	Riezler, 1941	L	EA	Sexual	—	No type specified by Riezler (1941).
Machilis	engiadina	Wygodzinsky, 1941a	A	EA	Parthenogenetic	59	Syntypes uninformative due to oxidized embedding medium; holotype (incomplete) and paratypea of synonymized M. distincta Janetschek (1970) at LMFI.
Machilis	fuscistylis	Riezler, 1941	N	EA	Parthenogenetic	29	Lectotypea (Janetschek, 1954) at LMFI.
Machilis	gepatschi	Riezler, 1941	L	EA	Parthenogenetic	—	No type specified by Riezler (1941).
Machilis	glacialis	Verhoeff, 1910	L	EA	Sexual	18(34b)	Type material “almost certainly lost” according to Mendes (1990).
Machilis	helleri	Verhoeff, 1910	none	CE	Sexual	40	Type material “almost certainly lost” according to Mendes (1990); holotype of M. helleri ssp. styriaca lost (this study); paratypea of M. helleri ssp. styriaca at LMFI.
Machilis	hrabei	Kratochvil, 1945	none	CE	Sexual	28	Type material lost according to K. Bezděčková (personal communication; MHJ).
Machilis	inermis	Wygodzinsky, 1941a	N	WA	Sexual	16 (14c)	All but one (male) syntypes uninformative due to dissolved tissue or oxidized embedding medium.
Machilis	ladensis	Janetschek, 1950a	L	EA	Sexual	10	Type material untraceable (this study).
Machilis	lehnhoferi	Riezler, 1941	N	EA	Sexual	47	Holotype lost; lectotype (Janetschek 1955) untraceable (this study).
Machilis	longiseta	Janetschek, 1949c	L	EA	Parthenogenetic	—	Holotype and paratypes at LMFI.
Machilis	mesolcinensis	Wygodzinsky, 1941a	L	WA	Sexual	20	Syntypes (NMB) uninformative due to dissolved tissue.
Machilis	montana	Wygodzinsky, 1941a	N	WA	Sexual	5	Syntypes (NMB) largely uninformative due to incomplete specimens or dissolved tissue; one female (Pr. IV/46 and 47)a intact.
Machilis	nigrifrons	Wygodzinsky, 1941a	L	WA	Sexual	—	Syntypes at NMB.
Machilis	oblitterata	Janetschek, 1970	L	EA	Sexual	—	Type material untraceable (this study).
Machilis	pallida	Janetschek, 1949a	N	EA	Parthenogenetic	43	Holotype untraceable (this study).
Machilis	pasubiensis	Bach de Roca, 1982	L	EA	Sexual	—	Type material at MSNV.
Machilis	pulchra	Janetschek, 1950b	N	EA	Parthenogenetic	11	Holotypes of M. pulchraa and M. pulchra var. silvestrisa at LMFI.
Machilis	robusta	Wygodzinsky, 1941a	L	WA	Sexual	11	Syntypes (NMB) uninformative due to dissolved tissue or oxidized embedding medium.
Machilis	rubrofusca	Janetschek, 1950a	N	EA	Parthenogenetic	37	Holotypea at LMFI.
Machilis	simplex	Wygodzinsky, 1941a	N	EA	Parthenogenetic	—	Syntypes at NMB.
Machilis	strebeli	Wygodzinsky, 1941b	L	EA	Parthenogenetic	—	Holotype at NMB.
Machilis	ticinensis	Wygodzinsky, 1941a	A	WA	Sexual	54	Syntypes (NMB) uninformative due to dissolved tissue or oxidized embedding medium.
Machilis	tirolensis	Verhoeff, 1910	A	WA + EA	Sexuald	42	Holotype “almost certainly lost” according to Mendes (1990).
Machilis	vagans	Wygodzinsky, 1941a	A	WA	Parthenogenetic	3	Syntypes not at NMB, untraceable (this study).
Machilis	vallicola	Wygodzinsky, 1941a	L	WA	Sexual	—	Syntypes not at NMB, untraceable (this study).
Machilis	vicina	Wygodzinsky, 1941a	L	EA	Sexual	9	Syntypes not at NMB, untraceable (this study).
Machilis	sp. B	Undescribed; this study	A	EA	Sexual	13	—
Machilis	sp. C	Undescribed; this study	Unknown	EA	Unknown	2	—
Machilis	sp. D	Undescribed; this study	Unknown	EA	Unknown	1	—
Machilis	sp. E	Undescribed; this study	Unknown	EA	Unknown	2	—

Notes: Eight species from the Western Alps were suspected to potentially occur in the Eastern Alps and were therefore included. Additionally, four potentially unknown species that could not be assigned to any nominal species are listed. Endemic status was estimated based on known distribution ranges following criteria commonly used for terrestrial arthropods (Rabitsch and Essl 2009). A = Alpine endemic; N = narrow-range endemic (< 10,000 km2); L = local endemic (< 5 km2); EA = Eastern Alps; WA = Western Alps; CE = Central European. Collections are abbreviated as follows: LMFI = Landesmuseum Ferdinandeum Innsbruck, Austria; MHJ = Museum of the Highlands Jihlava, Czech Republic; MSNV = Museo di Storia Naturale di Verona, Italy; MNHNP = Muséum National d’Histoire Naturelle Paris, France; NMB = Naturhistorisches Museum Basel, Switzerland.

aThese type specimens were examined by means of traditional morphometrics.

bEighteen specimens as M. glacialis, 34 as M. cryptoglacialis sp. n. (M. sp. A).

cFourteen specimens as M. cf. inermis.

dEven though sexual individuals of M. tirolensis are reported in the literature (Verhoeff, 1910; Wygodzinsky, 1941a), we exclusively found parthenogenetic populations in our study area. Information about reproductive modes refers to explicit statements in the literature.

Recently, morphological and mitochondrial markers were used to facilitate species delimitation in Machilis species (Dejaco et al. 2012). In the same study, several anonymous nuclear sequence markers failed to resolve relationships, probably due to highly admixed or duplicated genomes. More recently, Gassner et al. (2014) presented evidence for polyploidy and highly variable genome sizes in several Eastern Alpine Machilis species.

Here, we aim at enhanced species delimitation in Eastern Alpine Machilis jumping bristletails using a cohesive protocol of integrative taxonomy (Schlick-Steiner et al. 2010). Based on the insights from Dejaco et al. (2012) and Gassner et al. (2014), we apply traditional morphometrics, sequencing of one mitochondrial (cytochrome oxidase subunit 1; cox1) and one nuclear (internal transcribed spacer 2; ITS2) sequence marker, and AFLP genotyping to geographically representative samples of nominal Machilis species from throughout the study area. We explicitly do not include hypodermal pigmentation in our integrative protocol because it is not reliably quantifiable due to fuzziness and differences across molting stages (jumping bristletails keep on molting as adults, Sturm and Machida 2001). Eventually, we link our insights to existing names via original species descriptions, topotypic samples, and— whenever possible—via inclusion of type specimens into traditional morphometric analyses and describe two new species.

Materials and Methods

Sampling

We defined our study area, the Eastern Alps, based on the line connecting Lake Constance in the north and Lake Como in the south (the Rhine–Splügen line), which represents an established geological and biogeographical barrier between the Western and Eastern Alps (Grimm and Mattmüller 2004; Schönswetter et al. 2005; Schmitt 2009). Our goal was to sample the type localities of the 23 nominal species known from the Eastern Alps and 8 additional species from the Western Alps with a potential distribution across the Rhine–Splügen line (Table 1). Because of the absence of males in parthenogenetic species, we limited morphometric analyses to females but included males in molecular analyses. We aimed at a minimum of 15 females from three different populations per nominal species. In total, 105 populations were sampled (Fig. 1). For details on individual samples, including GPS coordinates and GenBank accessions, see Supplementary Table S1, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr. All specimens were sampled between 2010 and 2013 and stored at −20°C in 96% ethanol at the Institute of Ecology, University of Innsbruck, Austria.

F

Sampling sites. Species are color- and shape-coded, implemented identically in all other figures. The type localities of M. alpicola and M. helleri are outside our study area (see map inlay) but were nonetheless sampled. AT, Austria; CH, Switzerland; CZ, Czech Republic; DE, Germany; FR, France; IT, Italy. The dashed line roughly marks the Rhine–Splügen line that separates the Eastern from the Western European Alps.

Examination of Type Material

Type material of jumping-bristletail species is typically composed of dissected body parts permanently mounted on one or more glass slides per individual. All relevant type specimens that could be traced (see Table 1 for comments on the availability of type material) are kept in Landesmuseum Ferdinandeum Innsbruck (Austria; LMFI) or Naturhistorisches Museum Basel (Switzerland; NMB). The condition of holotypes, lectotypes, paratypes, or syntype series was determined for most of the species listed in Table 1. Female types in good condition were used to cross-check species identifications of collected samples and were included in the reanalysis of traditional morphometric characters.

Traditional Morphometrics

Of 574 collected specimens, 331 females were dissected and mounted on glass slides using water-soluble Marc-André 2 medium (Christiansen 1990). Eight meristic characters were determined using a Nikon Eclipse E600 bright field microscope. Character definitions and the resulting data matrix are given in Supplementary Table S2 and Supplementary Material 1, respectively, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr. In short, these characters represent count data from articulation and chaetotaxy of the ovipositor, coxites, and legs of females and had been identified as most informative for discriminating Machilis species in Dejaco et al. (2012). The principal component analysis (PCA) based cluster identification method of Ezard et al. (2010) was applied to the 331 mounted females to find significant clusters in our data set without a priori hypotheses. Discriminant analyses (original and leave-one-out cross-validated) based on final species limits were performed in SPSS v.21 (IBM, New York, USA) using data from the 331 mounted females and seven female type specimens. The latter were not included for calculating the discriminant function but were treated as unknowns in the reanalysis and were assigned to final species based on the discriminant function.

PCR Amplification, Sequencing, and Alignment

DNA was extracted from ethanol-preserved muscle tissue using the GenElute Mammalian Genomic DNA Miniprep Kit (Sigma–Aldrich, St. Louis, USA). Amplification of the partial cox1 gene necessitated multiple primer pairs: MachF1/MachR3 (Dejaco et al. 2012), MachF5/MachR7 (Gassner et al. 2014), and MachF4 (5′-ATTCGAGCTGAACTAGGNC-3′)/UEA10 (this study/Lunt et al. 1996). Ten-microliter reaction volumes contained ca. 50 ng template DNA, 1 × MyTaq Buffer, 0.2 μM of each primer, and 0.5 U MyTaq polymerase (Bioline Inc., London, UK). PCR conditions were 95°C for 2 min, 35 cycles of 94°C for 30 s, 50°C for 45 s, and 72°C for 90 s, and a final elongation at 72°C for 10 min.

For the nuclear sequence marker, ITS2, primers Bel28S and revBel28s (Belshaw and Quicke 1997) were used to amplify a DNA fragment spanning the 3′-end of the 5.8S domain, complete ITS2, and the 5′end of the 28S domain of the ribosomal DNA (rDNA) tandem repeat. A new forward primer (ITS-MachF2, 5′-GGGTCGATGAAGAACGCAGCTA-3′) was designed to improve cross amplification among species. The primer combination ITSMachF2/revBel28s was then used throughout this study. Ten-microliter reaction volumes contained ca. 50 ng template-DNA, 1 × PCR Buffer, 0.2 μM of each primer, and 0.5 U MyTaq polymerase. PCR conditions were 95°C for 2 min, 35 cycles of 94°C for 30 s, 62°C for 45 s, and 72°C for 105 s, and a final elongation at 72°C for 10 min.

Amplicons were checked by gel electrophoresis and purified by incubating 8 μL PCR product, 1 U Exo1, and 0.05 U FastAP (both Thermo Fisher Scientific Inc., Waltham, USA) at 37°C for 15 min followed by inactivation at 80°C for 15 min. Sanger sequencing was conducted by a commercial sequencing facility (Eurofins MWG Operon, Ebersberg, Germany) using the amplification primers. Electropherograms were visually inspected for quality and absence of double peaks. cox1 sequences were aligned manually and checked for correct amino acid translation. Two additional cox1 sequences from Trigoniophthalmus alternatus (GenBank accession: EU016193.1) and Petrobius brevistylis (GenBank accession: AY956355.1) were included as outgroups. ITS2 sequences were aligned using the X-INS-i algorithm in Mafft v.7.215 (Katoh and Standley 2013) and adjusted manually. Gaps were coded using the simple indel-coding algorithm (Simmons and Ochoterena 2000) implemented in GapCoder (Young and Healy 2003). Both the X-INS-i alignment algorithm, which incorporates information from pairwise structural alignments, and the simple indel-coding algorithm, which is known for its “biological realism” (Young and Healy 2003), maximize homology of gap characters. The resulting 1/0 matrix was appended to the ITS2 alignment. Identical sequences were removed from the final alignments before phylogenetic analyses, thus giving unequal sample numbers in the two alignments. Since the information contained in the two alignments was partly incongruent, we decided to build two independent gene trees rather than building one poorly resolved phylogeny. Both alignments are available (Supplementary Materials 2 and 3, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). All sequences were deposited in GenBank (Supplementary Table S1, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr).

AFLP Fingerprinting

A detailed description of the AFLP protocol applied is given in Wachter et al. (2012). Briefly, DNA samples of 516 individuals were digested with the enzymes MseI and EcoRI, and AFLP adaptors were ligated simultaneously. Preselective amplification was carried out with the primers Eco-N and Mse-C. Selective amplification used four primer combinations (tEco-ACA/Mse-CAA; tEco-ACA/Mse-CTT; tEco-ACC/Mse-CAA; and tEco-TA/Mse-CTC), with each forward primer having a 5′-M13 tail. Additionally, FAM/HEX/NED/PET labeled M13 primers were added in a ratio M13:F:R = 10:1:10 to enable subsequent pooling of amplicons from different primer combinations. At three different points in time, DNA of two to four individuals was replicated on each of the six 96-well plates used to verify identical conditions throughout the AFLP workflow: 1) before the digestion/ligation (four individuals), 2) before the preselective PCR (two), and 3) before the selective PCR (four), giving a total of 60 replicates. Fragment analysis was performed by a commercial provider (Comprehensive Cancer Center DNA Sequencing & Genotyping Facility, University of Chicago, USA) on an ABI 3730 sequencer (Applied Biosystems, Waltham, USA). Peak data were converted into .csv format with Peakscanner v.1.0 (Applied Biosystems). optiFLP v.1.51 (Arthofer et al. 2011) was used in unsupervised mode to identify parameters for automatic scoring (Supplementary Table S3, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). Final scoring with optimized parameters was performed in tinyFLP v.1.40 (Arthofer 2010). Samples with a peak frequency in the 10th percentile of the overall peak distribution in at least two primer sets were considered as PCR failures and excluded from further analysis (n=61, of which 12 were replicates). The four neighbor joining (NJ) trees (one for each primer combination; Supplementary Fig. S1, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr) produced by optiFLP based on Jaccard distances (Jaccard 1908) were used to identify potentially cross-contaminated profiles under the following rationale: To account for signals of hybridization and incomplete sorting, we allowed individual profiles to cluster with at most two different species across the four NJ trees, but excluded profiles that clustered with three or four different species (n=27). The refined matrix consisted of 488 unique profiles including 49 replicates and one outgroup. To estimate error rates, pairwise Euclidean distances between replicates of the same sample were calculated in Splitstree v.4 (Huson and Bryant 2006). Replicates were excluded in the final matrix (Supplementary material 4, available on Dryad at http://datadryad.org/review?doi=doi:10.5061/dryad.tf7qr). The outgroup was excluded for Bayesian clustering methods because—in contrast to phylogenetic trees—no root is needed for correct interpretation.

Phylogenetic Analyses of cox1, ITS2, and AFLPs

Phylogenetic trees were inferred using both Bayesian (MrBayes v.3.2; Ronquist et al. 2012) and Maximum-Likelihood (ML; Garli v.2.0; Zwickl 2006) approaches. In the Bayesian analyses, nucleotide substitution models were estimated by sampling over the entire GTR+G model space using the reversible-jump Markov chain Monte Carlo (rjMCMC) algorithm implemented in MrBayes (nst = mixed). By directly estimating substitution models, uncertainties in model selection can be mitigated according to the posterior probability of each model (Huelsenbeck et al. 2004). The cox1 alignment was partitioned into three codon positions as suggested by PartitionFinder v.1.1.1 (Lanfear et al. 2012), and model parameters were unlinked between partitions. The ITS2 alignment was partitioned into nucleotide and gap information. Gap characters were treated as standard data, and coding bias was corrected for (coding = variable) to account for absence of invariable characters. For both loci, two parallel runs, each consisting of four chains, were iterated for 3×106 generations. Trees were sampled every 1000 generations, and the first 25% were discarded as burn-in. Convergence was assured by an average standard deviation of split frequencies below 0.01 and accurate parameter estimates as indicated by estimated sample sizes above 200 and potential scale reduction factor values close to 1.

For the ML analyses, the cox1 alignment was partitioned into three subsets (1st, 2nd, and 3rd codon position) with three different substitution models (TrNef+G, HKY+I+G, and GTR+G, respectively) as identified via Bayesian Information Criterion using PartitionFinder. The non-coding ITS2 data set was partitioned into nucleotide and gap characters. For the nucleotide partition, the GTR+G model was identified as best-fitting substitution model using jModeltest v.2.1.4 (Darriba et al. 2012). The gap partition used the Mkv model implemented in Garli. Twenty independent ML searches were iterated in Garli, and the tree showing the highest likelihood was defined as best topology. To estimate node support, 100 bootstrap replicates were iterated in a separate Garli run, and bootstrap values were then summarized and mapped onto the best topology using SumTrees, included in the DendroPy software package v.3.12 (Sukumaran and Holder 2010).

Presence/absence data from AFLPs were coded as restriction data (1/0 matrix). In MrBayes, a correction for coding bias was applied using the “noabsencesites” option. After preliminary runs, the temperature of the heated chains was lowered to 0.01 to increase chain swapping frequency. For the state frequency parameter, a Dirichlet prior was used that resembled the actual state frequencies in our matrix (2 and 1 for “0” and “1”, respectively). Substitution rate was allowed to vary across sites. Two parallel runs, each consisting of one cold and five heated chains, were iterated for 30 × 106 generations. Trees were sampled every 5000 generations and the first 50% were discarded as burn-in.

Admixture Analyses

Bayesian clustering as implemented in the programs STRUCTURE v.2.3.4 (Hubisz et al. 2009) and BAPS v.6.0 (Corander and Marttinen 2006) was used to infer the number of species and the presence of interspecific admixture in the final AFLP data set. STRUCTURE is frequently used for species delimitation (Meudt et al. 2009; Leavitt et al. 2011; Satler et al. 2013; Hedin 2015) but can produce misleading results when sample sizes vary heavily among groups (Kalinowski 2011). In contrast, BAPS is less sensitive to sample-size variation and tends to detect finer genetic differentiation (Wilkinson et al. 2011). STRUCTURE runs were conducted for K=2 to K=27 with 10 replicates for each K using the admixture model with independent allele frequencies among populations. Each MCMC consisted of 3 × 105 generations, of which the first 105 were discarded as burn-in. Graphical representations of likelihoods across runs L(K) and K (Evanno et al. 2005) were generated using Structure Harvester (Earl and von Holdt 2012). The results for K=16 to K=19 were summarized using the Greedy algorithm in CLUMPP v.1.1.2 (Jakobsson and Rosenberg 2007) and visualized with DISTRUCT v.1.1 (Rosenberg 2004). In BAPS, clustering of individuals was performed with five replicates for K=1 to K=30. K=21 was identified as the optimum number of clusters and used in a subsequent admixture analysis using default parameters.

STRUCTURE was used to reanalyze admixture in a subset of the final AFLP matrix, including profiles of M. alpicola, M. engiadina, M. rubrofusca, M. ticinensis, M. tirolensis, M. vagans, and M. sp. E. For K=1 to K=7, each 10 replicate runs were iterated using the admixture model with independent allele frequencies among populations. Each MCMC consisted of 105 generations, of which 2×104 were discarded as burn-in. The same tools as mentioned above were used to generate admixture plots. The software NewHybrids v.1.1beta (Anderson and Thompson 2002) was used to test for hybrid genomes. Principal coordinates analysis (PCoA) was calculated in PAST v.3.06 (Hammer et al. 2001).

Species Delimitation Protocol

Defining the species concept and the species delimitation criteria used under each methodological approach is essential for reproducibility in integrative taxonomy. We used the unified species concept of De Queiroz (2007), which defines species as independently evolving metapopulations and acknowledges that different methods can delimit species at different stages of the speciation continuum. Nominal species as determined using the key of Palissa (1964) and putative morphospecies (see results) were used as primary species hypotheses. Unsupervised (i.e., hypothesis-free) analyses were then used to explore the variation present in the various data sets, thus generating secondary species hypotheses. Cases of incongruence among primary and secondary species hypotheses were later resolved and final species hypotheses formulated. The following criteria were used to build secondary species hypotheses:

Traditional morphometrics.—The delimitation criterion of phenotypic distinctness was implemented using the cluster identification method of Ezard et al. (2010). This approach accounts for the absence of multivariate normal distribution typically found in morphometric data sets containing multiple species. More specifically, the traditional morphometrics data matrix was centered on the median and scaled to the mean absolute deviation. The b-stick model was used as stopping criterion for dimension reduction.

Sequence markers.—Two different delimitation criteria were applied: The first, reciprocal monophyly (Donoghue 1985), was used despite the problems known to afflict it (Wiens and Servedio 2000; Kizirian and Donnelly 2004). Specifically, only nodes supported by Bayesian posterior probabilities >0.95 and ML bootstrap values >75 (Pattengale et al. 2009) were accepted as monophyletic. Reciprocity was interpreted on the level of primary species hypotheses to avoid instances of excessive nested diversity (Kizirian and Donnelly 2004), but nested entities were recognized in case of congruence with primary species hypotheses or with results from other data sets. The second criterion was the transition from species-level to population-level branching patterns identified using posterior delimitation probabilities as inferred by the Bayesian Poisson tree processes algorithm (Zhang et al. 2013). In contrast to the Generalized Mixed Yule Coalescent approach of Pons et al. (2006), the Bayesian Poisson tree processes algorithm models speciation in terms of number of substitutions instead of time units. Therefore, this model can be used when time-calibrated (ultrametric) trees are not available. The Markov chain was iterated for 105 generations using default settings for thinning and burn-in. Convergence was verified on the likelihood trace plot. Out of several possible species partitions, the partition with the highest likelihood found by a simple heuristic search was accepted as the delimitation result, and nodes supported by posterior delimitation probabilities >0.95 were accepted as species. For recently diverged lineages, however, the Bayesian Poisson tree processes model may return low posterior delimitation probabilities for both hypotheses (one vs. two species). In such instances, the one-species hypothesis was always favored as the delimitation result.

AFLP genotyping.—We are unaware of any algorithm that specifically identifies species limits on trees generated from AFLPs; two different delimitation criteria were used: First, reciprocal monophyly on the Bayesian consensus tree (see details above), and second, genotypic clustering (Mallet 1995) as inferred by BAPS and STRUCTURE. BAPS exclusively outputs the number of clusters with the highest marginal likelihood. For the STRUCTURE results, Evanno’s ∆K and the change in marginal likelihoods across runs [L(K)] were used to determine the optimum number of clusters.

Once secondary species hypotheses had been established, we sought evolutionary explanations for incongruence among data sets and reached agreement on final species hypotheses. Eventually, we used these final hypotheses to reanalyze the traditional morphometrics data set; in more detail, character combinations were sought that can be used for fast species identification. Specifically, we split our data set into four subgroups (Groups 1–4 in Fig. 6) according to ovipositor length (i.e., number of articles on gonapophysis VIII). This corresponded to the first steps in the Machilis identification key of Palissa (1964). For each group separately, we then used an exhaustive combination analysis that sequentially performed discriminant analyses using random characters from our data set and calculated overall classification rates for each possible character combination (Moder et al. 2007). We used the combinations with the lowest classification error and the lowest number of characters and calculated pairwise original and cross-validated classification rates using SPPS Statistics v.21 for each species group.

F

Summary of species delimitation. Numbers on the right are the sums of delimited species in each line. a) Primary species hypotheses. Each colored box represents 1 out of 20 sampled nominal species. White boxes with solid outlines represent additional four, potentially unknown species as defined in Table 1. For definition of M. cf. inermis, refer to section “Morphology-Based Identification (Primary Species Hypotheses)” in the Results and Discussion. b) Species delimitation. Each row summarizes secondary species hypotheses delimited by one data set. Colored boxes depict congruence with primary species hypotheses. White boxes represent secondary species hypotheses that were incongruent with primary species hypotheses. Horizontal splits specifically indicate populations in which individuals belonged to either one or another secondary species hypothesis; they do not indicate admixed individuals. The results based on AFLPs are given independently for analyses conducted in MrBayes, BAPS, and STRUCTURE. TM, traditional morphometrics; bPTP, Bayesian posterior tree processes; AFLP, amplified fragment-length polymorphism. c) Evolutionary explanations for cases of incongruence among primary and secondary species hypotheses of different data sets. d) The first row indicates final species hypotheses. Boxes representing previously unrecognized species have solid outlines. The second row lists classification success (in %) based on discriminant functions calculated separately from optimal combinations of traditional morphometric characters for four groups of morphologically similar species. LOOCV, leave-one-out cross-validation.

Results and Discussion

Morphology-Based Identification (Primary Species Hypotheses)

In total, 574 specimens were collected of which 496 were assigned to 20 nominal species following the identification key of Palissa (1964); see Table 1 for the taxon authorities. The other 78 specimens could not be assigned and were thus classified as follows: 1) In four species belonging to what is here termed the M. inermis species group (M. inermis, M. ladensis, M. robusta, and M. vicina), identification was hampered by the gradual nature of otherwise diagnostic characters (i.e., patterns of pigmentation). Even though we hypothesized synonymy already at this stage (see section “Morphological Oversplitting”), we assigned individuals from the respective type localities to nominal species (M. inermis: n=16, M. ladensis: n=10, M. robusta: n=11, M. vicina: n=9) and 14 specimens sampled from additional sites to a fifth category (M. cf. inermis). 2) Eighteen specimens sampled from the Southeastern Alps could not be determined at all and were hence assigned to four previously unknown morphospecies: M. sp. B (n=13), M. sp. C (n=2), M. sp. D (n=1), M. sp. E (n=2). In total, we defined 25 primary species hypotheses (see Fig. 6a).

Type Specimens Examined

We compiled up-to-date information on type material of 31 nominal species listed in Table 1. In 13 species, the types are lost or untraceable. For three species, no types were specified at all in the original descriptions. The type material of 15 species was examined. Of these, five species were missing in our collection of fresh samples, and the types were therefore not included in further analyses. In six species, the types were uninformative due to oxygenized embedding medium, dissolved tissues, or missing body parts. Only two holotypes (M. pulchra and M. rubrofusca), one syntype (M. montana), and one lectotype (M. fuscistylis) were informative. We additionally included the types of two subspecies (holotype of M. pulchra ssp. silvestris and paratype of M. helleri ssp. styriaca) and one paratype of an already synonymized species (M. distincta) in the reanalysis of traditional morphometric characters. Machilis distincta had been synonymized under M. engiadina by Janetschek (1970), and therefore the available M. distincta paratype was included instead of the uninformative syntypes of M. engiadina.

Species Delimitation (Secondary Species Hypotheses)

Traditional morphometrics.—Meristic characters were quantified in 331 females of 20 nominal species and M. sp. B. Nine distinct clusters were recognized without any a priori assumptions by the method of Ezard et al. (2010; Fig. 2a). Color-coding nominal species on the same PCA plot (Fig. 2b) revealed extensive overlap among them, with most clusters lumping two or more nominal species. Since most information on PC1 and PC2 stemmed from length and chaetotaxy of the ovipositor, these clusters mainly reflected differences in ovipositor morphology. PCA scatterplots from the reanalysis of traditional morphometric characters based on final species limits and including the type specimens listed in Table 1 are given in Supplementary Figure S2, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr.

F

Scatterplot of the first two principal components (PCs) extracted from eight traditional morphometric characters from 331 females. a) Individuals coded according to nine groups as predicted by the best model (ellipsoidal, equal volume, shape, and orientation) following the method of Ezard et al. (2010). b) Individuals coded according to nominal species except (*) M. sp. A (see results of molecular data).

Mitochondrial DNA.—Approximately 800 bp from the 5′-region of the cox1 gene were successfully amplified and sequenced from 563 individuals representing the 20 nominal species and M. sp. B to M. sp. E. The final alignment included 703 bp and was reduced to 111 unique haplotypes for phylogenetic inference. The Bayesian majority-rule consensus and the ML tree returned the same topology (Fig. 3). Basal nodes were weakly supported and therefore collapsed. According to the reciprocal monophyly criterion, 16 secondary species hypotheses were recovered, of which only six were congruent with primary species hypotheses. Ten secondary species hypotheses were incongruent with respect to 18 primary species hypotheses, rendering them paraphyletic either due to splitting (M. glacialis, where three clades were recovered, hereafter named M. glacialis, M. sp. A_sex., and M. sp. A_parth.) or sharing of haplotypes among primary species hypotheses (M. alpicola, M. engiadina, M. rubrofusca, and M. vagans; M. helleri and M. pulchra; M. helleri and M. lehnhoferi; M. helleri and M. hrabei; M. inermis, M. cf. inermis, M. ladensis, M. robusta, and M. vicina).

F

Bayesian majority-rule consensus trees based on cox1 and ITS2 sequences (111 and 73 unique sequences, respectively). Node support values indicate Bayesian posterior probabilities and ML bootstrap values. Species identified by the Bayesian Poisson tree processes algorithm are highlighted by red terminal branches. Colored bands visualize congruence/incongruence among gene trees.

Applying the Bayesian Poisson tree processes algorithm to the Bayesian consensus tree resulted in 20 secondary species hypotheses (Fig. 3 and Supplementary Fig. S3, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr), of which 9 were congruent with primary species hypotheses. Eleven secondary species hypotheses were incongruent with respect to 15 primary species hypotheses, rendering them paraphyletic either due to splitting (M. glacialis, M. hrabei, and M. sp. B) or lumping (M. helleri and M. pulchra; M. glacialis and M. montana; M. alpicola, M. engiadina, M. rubrofusca, and M. vagans; M. inermis, M. ladensis, M. robusta, and M. vicina). Machilis sp. C, M. sp. E, and M. sp. D were recovered as separate species.

Nuclear ribosomal DNA.—Between 600 and 800 bp of ITS2 from 562 individuals were successfully amplified and sequenced. The edited sequence alignment included 73 unique sequences with 761 characters and numerous gaps of variable length (1–98 bp). Consequently, 219 gaps were coded, and the resulting 1/0 matrix was appended to the final alignment (Supplementary Material 4, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). The Bayesian majority-rule consensus and the ML tree returned the same topology (Fig. 3) and resolved more nodes compared with the cox1 tree. According to the reciprocal monophyly criterion, 20 secondary species hypotheses were recovered, of which 8 were congruent with primary species hypotheses. Twelve secondary species hypotheses were incongruent with respect to 16 primary species hypotheses, rendering them paraphyletic either due to splitting (M. alpicola, M. glacialis, M. helleri, and M. ticinensis) or lumping (sexual M. ticinensis and M. engiadina; M. alpicola, M. engiadina, M. rubrofusca, and M. vagans; M. hrabei and M. lehnhoferi; M. helleri and M. pulchra; M. aleamaculata and M. montana; M. inermis, M. cf. inermis, M. ladensis, M. robusta, and M. vicina). Machilis sp. B, M. sp. C, M. sp. D, and M. sp. E were recovered as separate species.

The Bayesian Poisson tree processes algorithm delimited 21 secondary species hypotheses (Fig. 3 and Supplementary Fig. S3, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr), of which 7 were congruent with primary species hypotheses. Fourteen secondary species hypotheses were incongruent with respect to 17 primary species hypotheses rendering them paraphyletic either due to splitting (M. alpicola, M. glacialis, M. helleri, and M. lehnhoferi) or lumping (M. alpicola, M. engiadina, M. rubrofusca, and M. vagans; M. aleamaculata and M. montana; M. engiadina and M. ticinensis; M. helleri and M. pulchra; M. hrabei and 12 individuals of M. lehnhoferi). Machilis sp. B, M. sp. C, M. sp. D, and M. sp. E were recovered as separate secondary species hypotheses.

AFLP genotyping.—The final alignment comprised 488 profiles including 49 replicates and one outgroup. An error rate of 23.4% (SD = 7.2%) was calculated, corresponding to the mean pairwise Euclidean distance among replicates of four different individuals. Similar error rates were reported in a study investigating the effect of automated scoring algorithms on AFLP marker quality (Holland et al. 2008). In the same study, Holland et al. showed that there is a tradeoff between scoring few high-quality markers and more low-quality markers and that the latter strategy potentially increases phylogenetic resolution. In line with this, using more stringent scoring parameters resulted in fewer markers but reduced the phylogenetic resolution in our data set (data not shown).

The Bayesian majority-rule consensus tree was almost fully resolved, but several inner nodes were weakly supported (Fig. 4). Resolution at the species level was good, though. Sixteen secondary species hypotheses were recovered following the reciprocal monophyly criterion, six of which corresponded to primary species hypotheses. Ten secondary species hypotheses were incongruent with respect to 13 primary species hypotheses, rendering them paraphyletic: Machilis sp. A was split from M. glacialis, and parthenogenetic populations of M. sp. A were nested within the sexual population of M. sp. A. Two reciprocally monophyletic secondary species hypotheses were recovered within M. helleri, corresponding to western and eastern populations, in line with entities recovered by the ITS2 data and AFLPs analyzed with BAPS (see paragraph below). A second nominal species, M. pulchra, was nested within the western clade of M. helleri. Machilis aleamaculata and M. montana formed a monophyletic group. Machilis ticinensis was split in two clades corresponding to one southern population (hereafter sexual M. ticinensis), and four northern populations (hereafter parthenogenetic M. ticinensis). The samples of M. inermis, M. ladensis, M. robusta, and M. vicina were not recovered as separate monophyla but clustered in one group that was not supported. Machilis engiadina individuals clustered on the branch in between sexual M. ticinensis and the species pair M. rubrofusca/M. alpicola.

F

Bayesian majority-rule consensus tree based on 466 AFLP markers from 438 individual profiles. Node support (posterior probability >0.95) is indicated by filled circles. Reciprocally monophyletic groups (secondary species hypotheses) that correspond to primary species hypotheses are displayed with colored branches. Secondary species hypotheses incongruent with primary species hypotheses have colored backgrounds. The single terminal branch basal to M. helleri (marked by *) corresponds to the hybrid individual from population “Leopoldsteiner See” (M. helleri × M. hrabei) as inferred via the NewHybrids software.

BAPS returned 21 secondary species hypotheses and almost no admixture (Fig. 5). Eight secondary species hypotheses corresponded to primary species hypotheses (Fig. 6b). In three instances, primary species hypotheses were split into northern and southern (M. ticinensis) or into eastern and western populations (M. hrabei and M. helleri). The western lineage of M. helleri included all samples of M. pulchra. Machilis glacialis was split in three secondary species hypotheses, corresponding to M. glacialis s. str., sexual M. sp. A, and parthenogenetic M. sp. A. The samples of M. inermis, M. ladensis, M. robusta, and M. vicina were assigned to two separate clusters, not reflecting any apparent geographical or altitudinal pattern.

F

BAPS and STRUCTURE admixture plots based on 466 AFLP markers from 438 individual profiles. Colors represent secondary species hypotheses recovered for each value of K. Black bars and names above the upper plot indicate primary species hypotheses (i.e., nominal species). Additional black bars and names below the plots label instances of incongruence among primary and secondary species hypotheses.

STRUCTURE runs first converged inconsistently, but this problem disappeared when the locprior model (Hubisz et al. 2009) was applied, using nominal species affiliations as prior population information. Both the K statistic and the mean likelihood across runs of the same K indicated 19 separate clusters (Supplementary Fig. S4, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). Background noise was prominent in the admixture plot, resulting in two biologically arbitrary groups (gray and dark gray in Fig. 5). We considered this artifact to be a result of our strategy to call a high number of low-quality peaks instead of few high-quality peaks. Thus, we ignored the two arbitrary groups and accepted the remaining 17 groups as the most conclusive result of STRUCTURE. Nine groups corresponded to primary species hypotheses (Fig. 6b). Splits of primary species hypotheses were supported in M. ticinensis (sexual M. ticinensis and parthenogenetic M. ticinensis) and in M. glacialis (M. glacialis and M. sp. A). The following primary species hypotheses were lumped: M. helleri and M. pulchra, M. aleamaculata and M. montana, and one group including four nominal species (M. inermis, M. ladensis, M. robusta, and M. vicina). Admixture was found in M. alpicola (from M. rubrofusca) and in sexual M. ticinensis (from M. tirolensis, M. engiadina, and M. alpicola). Moreover, there was a weak but consistent signal of admixture from M. engiadina into M. rubrofusca. An inconsistent signal of admixture from several different groups into the sexual population of M. sp. A might indicate genetic differentiation from parthenogenetic populations of M. sp. A (as recovered by cox1, the AFLP tree, and BAPS) rather than true admixture events.

Evolutionary Explanations for Incongruence, Final Species Hypotheses, and Taxonomic Consequences

We uncovered multiple cases of incongruence among species limits inferred with different methodological approaches (Fig. 6), possibly reflecting complex evolutionary histories in several Machilis species. Not a single nominal species was recovered by unsupervised analysis of the traditional morphometrics data set, and only 3 out of 20 primary species hypotheses were congruent across molecular data sets. In the following, we discuss instances of incongruence, present evolutionary explanations for the observed patterns, and suggest taxonomic changes including six new synonymies and the description of two new species.

Hybridization and parthenogenesis in the M. engiadina species group.—The six species M. alpicola, M. engiadina, M. rubrofusca, M. ticinensis, M. tirolensis, and M. vagans were morphologically very similar and are here allocated to what we term the M. engiadina species group. The group also includes M. sp. E, based on its position in the PCA plot of traditional morphometric characters and in the trees generated from cox1, ITS2, and AFLPs. Since M. alpicola and M. vagans were consistently lumped in all data sets, we declare M. alpicola a junior synonym of M. vagans, following the priority rule of the ICZN (1999). In principle, the failure to delimit M. vagans as separate species could also result from the low sample size (n=3), since this is a known issue of Bayesian clustering methods such as STRUCTURE (Kalinowski 2011). However, BAPS delimited other small-sample clusters in our data set that are clearly biologically relevant (e.g., sexual and parthenogenetic populations of M. cryptoglacialis; M. sp. G and M. hrabei). We are thus confident that M. vagans and M. alpicola constitute just one species and advocate the use of BAPS alongside STRUCTURE in future delimitation studies, independent of the taxonomic rank under investigation.

Reanalysis of the traditional morphometrics data set tentatively indicated two overlapping clusters, with M. tirolensis and M. ticinensis on one side, M. engiadina, M. rubrofusca, and M. vagans on the other side, and M. sp. E in between (Supplementary Fig. S2a, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). This is congruent with the phylogenetic pattern observed in the cox1 tree (Fig. 3). In the ITS2 tree, however, some M. engiadina individuals clustered with M. alpicola and M. rubrofusca, others with sexual M. ticinensis. This mito-nuclear discordance is unlikely to be the result of incomplete lineage sorting because M. rubrofusca/M. vagans and sexual M. ticinensis are distantly related. Rather, we discuss three different hybridization scenarios (Fig. 7a–c).

F

Deciphering admixture in the Machilis engiadina species complex. Three hypotheses explaining incongruence among sequence markers in M. engiadina, a) without taking into account signals of admixture inferred from AFLP markers, b) under the assumption of four independent hybridization events, and c) under the assumption of fixation of extant lineages from a hybrid swarm. d–e) Principal coordinates plots based on 466 AFLP markers from 158 individuals. Gray dashed lines indicate the triangle of putative parental species—a pattern described by Ma and Amos (2012). f) STRUCTURE plots for K=2 to K=6. g) Triangle plot from the STRUCTURE results at K=3. h) Distribution map including occurrence data from this study and the literature.

Scenario 1: Hybridization between the M. rubrofusca/M. vagans lineage and sexual M. ticinensis, followed by the persistence of a hybrid lineage (i.e., M. engiadina), could represent a reasonable explanation. We found statistical support for this scenario using a subset of the final AFLP matrix. The software NewHybrids identified the samples of M. engiadina as either F1 or F2 hybrids when M. rubrofusca and sexual M. ticinensis were predefined as parental species (for details, see Supplementary Table S4a, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). However, when M. tirolensis, parthenogenetic M. ticinensis, and M. sp. E were included as unknowns, they were identified as paternal lineages alongside sexual M. ticinensis with equally high posterior probabilities (Supplementary Table S5b, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). With the data at hand, we cannot evaluate whether this is due to the close relation of M. ticinensis, M. tirolensis, and M. sp. E or a more complex evolutionary history, possibly involving repeated hybridizations among more than two species. In the absence of a formal test for multiple hybridizations, we used the results of descriptive statistics and admixture analyses to discuss two additional scenarios.

Scenario 2: A principal coordinates analysis of the reduced AFLP matrix revealed four major clusters (Fig. 7d). By comparing the two-dimensional pattern (PCo1 vs. PCo2) with the results of a simulation study of Ma and Amos (2012), we built the following hypotheses: First, M. rubrofusca, M. tirolensis, and parthenogenetic M. ticinensis formed a triangle of parental species. Second, M. engiadina was a hybrid (50:50; H1 in Fig. 7b) of M. rubrofusca and M. tirolensis since it lies halfway on the line connecting these two species. Third, M. vagans was a backcrossed hybrid (~80:20; H3 in Fig. 7b) of M. rubrofusca and parthenogenetic M. ticinensis since it lies on the line connecting these two species. Fourth, sexual M. ticinensis was a backcrossed hybrid (H2 in Fig. 7b) between M. tirolensis, parthenogenetic M. ticinensis, and a third species (M. rubrofusca, M. vagans, or M. engiadina) since it almost completely lies within the triangle. Machilis engiadina is an especially good candidate for the third parent species since it shares some AFLP loci with sexual M. ticinensis (those represented by PCo3; Fig. 7e) but none with M. tirolensis. Fifth, M. sp. E was a backcrossed hybrid (H4 in Fig. 7b) of M. tirolensis and another species from within this species group (possibly the ancestor of M. ticinensis, to which its cox1 sequence is most closely related). These hypotheses are congruent with the results of an admixture analysis in STRUCTURE (Fig. 7f). All lineages except M. sp. E were recognized by the clustering algorithm (K=7), which would not be expected from young hybrids but could be due to separate evolutionary histories of the hybrids following hybridization. A triangle plot (Fig. 7d) captures the overall picture when three parental populations are assumed, basically mirroring the pattern of the first two principal coordinates (Fig. 7a). We thus hypothesize four independent hybridization events and consequently four instances of hybrid speciation (H1–H4 in Fig. 7b). This scenario provides a reasonable explanation for the mito-nuclear incongruence in M. engiadina and for the admixture proportions seen in the STRUCTURE plots. Hybrid speciation is increasingly accepted as a promoter of diversity in animals (Herder et al. 2006; Dasmahapatra et al. 2012; Papadopulos et al. 2013; Andújar et al. 2014) but difficult to demonstrate. Besides genetic evidence, claiming hybrid speciation requires evidence for reproductive isolation of the putative hybrid lineage from the parental species and, ideally, evidence that isolation was a consequence of the hybridization event (Schumer et al. 2014). Here, reproductive isolation of the putative hybrid lineages has been warranted by parthenogenesis in two species (M. engiadina and M. vagans). Moreover, switches to parthenogenesis are known to occur frequently after hybridization in animals (Suomaleinen et al. 1987; Paland et al. 2005; Gómez-Zurita et al. 2006). In the third putative hybrid species (sexual M. ticinensis), reproductive isolation may be only temporarily warranted by parapatry with respect to its parental species (M. tirolensis, parthenogenetic M. ticinensis, and possibly M. engiadina; Fig. 7h). In times of sympatry, however, reproductive isolation may have been absent. Possibly, incomplete reproductive isolation of sexual M. ticinensis from its parental species could even explain why this lineage shows admixture with a third species.

Scenario 3 (Fig. 7c) differs from Scenario 2 in that only one ancestral hybridization event is postulated (e.g., between the ancestors of M. rubrofusca and M. tirolensis). When two distantly related genomes fuse, hybrid swarms consisting of lineages with different admixture proportions can emerge (Cockayne 1926; Stebbins 1959; Seehausen 2004). Given available ecological-niche space, hybrid swarms can act as cradle of novel, highly adapted genotypes and are therefore discussed as potential precursors of adaptive radiations (Seehausen et al. 2003; Hudson et al. 2010; Genner and Turner 2012). The extant representatives of the M. engiadina species group may represent lineages that emerged from a hybrid swarm and became fixed via parthenogenetic reproduction, while most other lineages went extinct. The fact that the two putative ancestral lineages were distantly related and the cline-like pattern at K=2 in the admixture plot may be interpreted as in favor of this hybrid swarm scenario.

We cannot here evaluate the likelihood of Scenarios 1–3. However, this does not impede our inference that M. engiadina originated via homoploid hybrid speciation, since M. engiadina and sexual and parthenogenetic M. ticinensis share the same ploidy level (Gassner et al. 2014). Homoploid hybrid speciation is thought to occur rarely compared with polyploid hybrid speciation (but see Schwarz et al. 2005; Gompert et al. 2006; Mavárez et al. 2006; Nice et al. 2013), because a clear reproductive barrier against both parental species is missing in the absence of polyploidy (Mallet 2007). Consequently, homoploid hybrid speciation becomes more likely when hybrids find alternative ecological niches or become reproductively isolated via asexual reproduction (Rieseberg et al. 2003; Abbott et al. 2013). Both conditions are met in M. engiadina and M. vagans: First, their distributional ranges extend beyond those of their putative parental species (Fig. 7h) indicating a shift in their potential ecological niche space, and second, no males have been found so far.

Overall, the M. engiadina species group may be an emerging example of how parthenogenesis can stabilize hybrid genotypes and underpins the idea that hybrid and ecological speciation can mutually facilitate biological diversification (Seehausen 2004; Kearney 2005). In terms of species delimitation, hybridization provides a reasonable explanation for the discordance found among secondary species hypotheses for the M. engiadina species group. Moreover, our hybridization scenarios support the persistence of parthenogenetic hybrid lineages as independently evolving units and thus as separate species. For now, we set aside the formal description of parthenogenetic M. ticinensis (henceforth M. sp. F) as a new species, mainly because sexual M. ticinensis was represented only by one population. Until additional samples of sexual M. ticinensis become available, we cannot exclude the possibility that we sampled an aberrant population of M. ticinensis and that parthenogenetic M. ticinensis (M. sp. F) is in fact congruent with M. ticinensis s. str. We also postpone the formal description of M. sp. E until further material becomes available.

Mito-nuclear discordance in the M. helleri species group.— Three nominal species, M. lehnhoferi, M. helleri, and M. hrabei, belong to what we term the M. helleri species group. We found incongruence between mitochondrial and nuclear trees (see Fig. 3) in several populations in this group. We first tested for signals of admixture by separately comparing the AFLP profiles of the three possible species pairs. NewHybrids unambiguously assigned all putative hybrid individuals to pure parental species (Supplementary Table S5, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr), thus rejecting the hybrid hypothesis. However, in the comparison of M. helleri and M. hrabei, one individual of M. helleri (from population “Leopoldsteiner See”) had a higher posterior probability for having a backcrossed genotype than for being pure M. helleri. Even though we do not consider recent hybridization as the reason for the extensive mito-nuclear discordances in this species group, the one putatively backcrossed individual may indicate recent or ongoing hybridization between M. helleri and M. hrabei.

Alternative explanations for the mito-nuclear discordance in this species group are 1) incomplete lineage sorting, 2) introgression of mitochondrial or ribosomal DNA, and 3) phylogenetic error. Incomplete lineage sorting is a stochastic process that affects populations randomly over the entire geographical distribution (Funk and Omland 2003; Toews and Brelsford 2012). By contrast, all populations of M. lehnhoferi and M. hrabei that harbored a M. helleri-like mtDNA haplotype were located at the margins of the species’ geographical distribution (Supplementary Fig. S5, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). In accordance with established patterns of mitochondrial introgression (Toews and Brelsford 2012), we therefore conclude that mitochondrial haplotypes of M. helleri occasionally introgressed into both M. hrabei and M. lehnhoferi. We consider phylogenetic error to be an unlikely explanation for incongruence here due to the large gap observed between intra- and interspecific genetic distances among cox1 haplotypes.

The M. hrabei-like ITS2 sequences found in two populations of M. lehnhoferi could be explained by incomplete lineage sorting, introgression of ribosomal DNA, or phylogenetic error in the ITS2 tree. It is almost impossible to disentangle these effects in phylogenetic reconstructions of young sister species. In the absence of any signal of admixture in the AFLP data set and any geographical pattern in line with a historical or recent contact zone (Supplementary Fig. S5, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr), all three scenarios seem equally plausible.

In terms of species delimitation, incomplete lineage sorting, mitochondrial and nuclear introgression, and phylogenetic error provide reasonable explanations for the patterns of incongruence in our data, and thus, the species status of M. helleri, M. hrabei, and M. lehnhoferi is confirmed.

Additionally, our data indicate the possibility of ongoing speciation between western (including M. pulchra) and eastern populations of M. helleri. Even though western mtDNA haplotypes are still present in some populations of the eastern lineage, Gassner et al. (2014) found a consistent difference in chromosome numbers (western: n=52; eastern: n=50). Thus, both lineages potentially constitute independently evolving metapopulations as defined under the unified species concept. The western lineage corresponds to M. helleri as it includes topotypic samples and includes M. pulchra consistently across all data sets. We therefore declare M. pulchra a junior synonym under M. helleri, following the priority rule of the ICZN (1999). The eastern lineage likely corresponds to M. aciliata (see Table 1) based on morphological similarity of our samples to the original description of M. aciliata (Janetschek 1955). We were not able to obtain topotypic samples of M. aciliata and therefore refer to the eastern lineage of M. helleri as M. sp. G until further material has become available.

Morphological oversplitting.—According to their original descriptions, M. inermis, M. ladensis, M. robusta, and M. vicina differ from each other in pigmentation patterns on legs, head, and maxillary palps. However, these descriptions were based on few individuals from single or few populations (Wygodzinsky 1941a). By sampling multiple populations between type localities, we found that the pigmentation patterns did not fit discrete categories. We therefore suspected synonymy of the four nominal species already when defining primary species hypotheses. This hypothesis has been corroborated by the results of traditional morphometrics and all molecular analyses. We give M. inermis priority over the other three names and declare M. ladensis, M. robusta, and M. vicina as junior synonyms under M. inermis. We also point out that four additional species descriptions (i.e., M. anderlani, M. nigrifrons, M. vallicola, and M. oblitterata) potentially fall in the morphological variation now present in M. inermis. Unfortunately, no samples were found at the type localities of these nominal species (see Table 1). Future studies should be alert to potential synonyms of M. inermis among Machilis species characterized by a primary ovipositor also outside the Eastern Alps, that is, M. zangherii Janetschek, 1949, M. macedonica Janetschek, 1957, M. lindbergi Wygodzinsky, 1959, and M. ingens Bitsch, 1963. Should synonymy with any of these apply, the Alpine endemism of M. inermis will need to be reconsidered.

Cryptic species uncovered.—Both sequence markers and AFLPs supported the presence of at least two species within M. glacialis. No type material was available (Table 1), but allocating the name glacialis to a biological entity is feasible because M. glacialis was exclusively described from its type locality. Moreover, samples from the type locality fitted the original description given by Verhoeff (1910), and no additional Machilis species was found at this location. Consequently, we assign the name glacialis Verhoeff, 1910 to the final species hypothesis including topotypic samples and—temporarily—term the other, morphologically cryptic final species hypothesis M. sp. A. Reanalysis of M. sp. A using traditional morphometrics (Supplementary Table S3, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr) and qualitative morphology (see formal description of M. cryptoglacialis sp. n. in Appendix 1) revealed consistent differences from M. glacialis, therefore uncovering it from morphological crypsis. Machilis sp. A included two mitochondrial subclades divided by 9.6% mean net-between-group distance that were recognized as separate secondary species hypotheses by the Bayesian Poisson tree processes model. These subclades corresponded, on the one hand, to individuals from five Central and, on the other hand, one Southern Alpine population. While in the former only females were found (n=26), four males were sampled in the latter (n=8), possibly indicating a pattern of geographical parthenogenesis (Vandel 1928). Clonal reproduction is further supported by the homogeneity of cox1 sequences among Central Alpine populations compared with the observed sequence diversity in the sexual population. Even though southern and central populations each had private ITS2 sequences, reciprocal monophyly was not supported, and one male from the southern population was recognized as a separate secondary species hypotheses by the Bayesian Poisson tree processes model. Similarly, reciprocal monophyly was not supported in the Bayesian phylogeny of AFLPs, while Bayesian clustering (BAPS) recovered Central and Southern Alpine lineages as separate secondary species hypotheses. In the absence of additional samples from the area between central and southern populations, we decline to describe two separate species. Instead, we erect one new species, M. cryptoglacialis sp. n. (see printed Appendix 1), and highlight the possibility of ongoing speciation among Central and Southern Alpine populations.

Previously unrecognized species beyond cryptic diversity.— Our sampling strategy included areas rarely sampled for jumping bristletails, especially in the Southern and Southeastern Alps. We thereby sampled several previously unknown Machilis species. Among these, just M. sp. B was sampled repeatedly and is here described as M. albida sp. n. (see printed Appendix 1). In contrast, M. sp. C, M. sp. D, and M. sp. E were sampled just once each, and we thus set aside any potential formal description until more samples become available to characterize morphological and genetic variation. Concluding, we highlight the underestimated diversity in the genus Machilis and the potential importance of the Southern Alpine area as a hotspot of Alpine endemic biodiversity.

Geographical parthenogenesis.—In several species, Central Alpine populations were putatively all-female while conspecific populations at the eastern or southern margin of the Alps were bisexual, and this pattern often reflected genetic or morphological differentiation or both. In M. ticinensis, Northern Alpine parthenogenetic populations (M. sp. F) may have experienced a different evolutionary history of hybridization compared with Southern Alpine sexual M. ticinensis. In M. cryptoglacialis sp. n., Central Alpine parthenogenetic populations differed strongly in mtDNA and traditional morphometric characters from Southern Alpine sexual individuals. In M. mesolcinensis, Central Alpine parthenogenetic populations differed strongly in traditional morphometrics from Southern Alpine sexual individuals. Within M. helleri, individuals previously determined as M. pulchra syn. n. were restricted to the Central Alps and are putatively parthenogenetic, while sexual M. helleri primarily occurred along the Eastern Alpine foothills. These patterns are in line with the concept of geographical parthenogenesis (Vandel 1928). It is also well-known that peripheral Alpine areas served as refuges (so-called massifs de refuge) during Pleistocene glaciations (Holdhaus 1954; Pauls et al. 2006; Margraf et al. 2007) and that recolonization of previously devastated areas has often been associated with increased frequency of asexuality in several animals and plants (Kearney et al. 2006; Hörandl 2009). Consequently, we hypothesize that several Alpine Machilis species survived Pleistocene glaciations in southern or eastern massifs de refuge, where they retained sexual reproduction. At the same time, parthenogenetic populations may have survived on so-called nunataks (see Wachter et al. 2012) or may have recolonized Central Alpine areas following deglaciation. Either way, these mechanisms may have triggered the genetic and morphological divergence between Central Alpine and peripheral populations, thus promoting diversification and speciation in the genus Machilis.

Reanalysis of Traditional Morphometric Data

Using the final species hypotheses, we identified traditional morphometric character combinations that yielded the highest possible classification success for individual specimens in four species groups. Optimal combinations included three to eight characters and produced group-specific classification rates between 77% and 97% (for details, see Supplementary Table S6, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). Classification rates specific for final species hypotheses ranged between 53% and 100% (original), and 29% and 100% (cross-validated). Six type specimens of four species (M. engiadina, M. fuscistylis, M. helleri, and M. rubrofusca) were correctly classified based on the discriminant function (Supplementary Table S7, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). For three of these species (M. fuscistylis, M. rubrofusca, and M. engiadina) which are sympatric with at least some of the morphologically most similar species (Fig. 1; Supplementary Fig. S2, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr), the classification results are especially relevant, in that the names to use for the collected material are herewith confirmed (the fourth species, M. helleri, is not sympatric with the most similar species, and therefore there is no doubt about the name). In contrast, the syntype of M. montana was wrongly classified as M. glacialis (Supplementary Table S7, available on Dryad at http://dx.doi.org/10.5061/dryad.tf7qr). We argue that the morphometric characters used here are just not able to reliably distinguish M. montana from M. glacialis. Anyway, this has no bearings on the species status of samples identified as M. montana in this article—the two species are clearly distinct in all molecular data (Figs. 3–6)—nor on the names to use for the collected specimens since the two species are not sympatric and geographically well separated by the Rhine–Splügen line. Classification rates were generally lower in species groups 1 and 2, where hybridization and introgression were frequently found. Overall, the optimal character combinations improve the quantitative assessment of species identities compared with subjective descriptions of qualitative morphological characters like pigmentation patterns. However, they also highlight the morphological similarity among Alpine jumping bristletails in the genus Machilis.

Conclusion

In this study, we revised species limits in a taxonomically challenging arthropod group. Of 25 primary species hypotheses investigated, none was recovered by traditional morphometrics, and only three were supported by full congruence among molecular data sets. We found high levels of mito-nuclear discordance for multiple species and identified evolutionary explanations for the observed patterns. Hybridization is inferred to be common in Alpine Machilis species. In the M. engiadina species group, past hybridization events combined with switches to parthenogenesis may have generated mosaic genotypes and thus promoted lineage diversification via hybrid speciation. In the M. helleri species group, widespread patterns of mitochondrial introgression and at least one backcrossed individual provide evidence for ongoing hybridization.

We describe two new species, of which one is uncovered from morphological crypsis, and highlight the potential presence of four additional species to be described, three of them from the Southern Alps. Patterns of geographical parthenogenesis mirrored instances of genetic and morphological divergence among Central Alpine and marginal Alpine populations in several species, indicating that they have been affected by glacial cycles during the Pleistocene. Ultimately, repeated demographic shifts imposed by climatic oscillations together with these species’ predisposition for hybridization and parthenogenesis may provide a framework for explaining the outstanding amount of Alpine endemism in the genus Machilis. Based on the presence of up to six previously unrecognized species in our data set, we assume that the diversity of Alpine Machilis species is still underestimated, especially in the Southern Alpine region. We stress the importance of describing this diversity for its conservation.

Supplementary Material

Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.tf7qr.

Acknowledgments

The authors thank Barbara Thaler-Knoflach, Gregor A. Wachter, Lukas J. Rinnhofer, Hannes Rauch, Michael Url, Dominik Wagner, Andrea Lassacher, Patrizia Arthofer-Peer, and Gernot Kunz for help with field sampling. They also thank Clemens Folterbauer for support in the wetlab, Herbert Zettel and Erhard Christian for advice on applying the ICZN Code to specific taxonomic problems, Vito Zingerle for support with fundraising, and the Autonomous Province of South Tyrol (Amt für Landschaftsökologie) and the administration of the Swiss National Park for granting sampling permits. Gerhard Tarmann and Daniel Burckhardt helped locating type specimens and granted loans. Luis F. Mendes supported T.D. with taxonomic advice throughout this project. Frank E. Anderson, Karl M. Kjer, Felix A. H. Sperling, and an anonymous reviewer offered important comments on an earlier version of the article.

Funding

This project was funded via stipends of the University of Innsbruck [to T.D. and M.G.] and by the Autonomous Province of South Tyrol [project-ID: 1/40.3; 27 January 2014].