Subterranean Deuteraphorura Absolon, 1901, (Hexapoda, Collembola) of the Western Carpathians — Troglomorphy at the northern distributional limit in Europe

An integrative approach employing molecular, morphological and geographical data were applied to species delimitation among Deuteraphorura congeners occupying caves of the Western Carpathian Mts. A new species of Deuteraphorura from the Western Carpathians is described. D. muranensis sp. nov. belongs among species with 4 pso at the hind margin of the head and possesses highly troglomorphic features. It is conspicuous with its distinctly elongated claws and long, hair-like body chaetae. The status of the new species was confirmed by DNA barcoding based on the mitochondrial COI marker. Populations of D. kratochvili (Nosek, 1963), the most widespread species, were studied in detail. Both ABGD and PTP analyses brought results congruent with geography, i.e. the molecular and geographic distance of the populations were positively correlated. However, some molecular separation based on pairwise distance and the number of substitutions was indicated within two of the studied populations. Despite the indistinct morphological differences, the tested populations were well isolated both geographically and genetically, which indicates that each studied population may represent a cryptic species. The troglomorphy of cave Collembola at the northernmost border of the distribution of cave-adapted species in the Europe is discussed. It is clear that the level of troglomorphy is closely associated with conditions of the microhabitat occupied by the individual subterranean species. The results of our study enhance the importance of the Western Carpathians regarding the diversity pattern of obligate cave species in Europe.


Introduction
In Europe, a high diversity of cave-adapted species has been observed in so-called "hot-spot" areas of the southern European mountains, with a decreasing trend towards northern regions, where only a few troglobionts have been documented [1,2,3]. The distribution of obligate cave species in the Western Carpathian Mts, i.e. beyond troglobiont-rich areas, points out the importance of these karst regions as the northernmost distributional range borders of terrestrial troglobionts [4]. Caves represent extreme environments with specific conditions that generate selection pressure on cave inhabitants and lead to the convergence of morphological traits [5]. Variation in these features, even in closely related cave species, is a function of microhabitat partitioning, distinct especially in aquatic environments [6]. Changes in morphology associated with cave life may affect several structures related to darkness, food scarcity, air temperature and the presence of wet clayey sediment. Such changes are manifested in body size, body slenderness, eyes, pigment, appendage length, foot complex and wing development as progressive (the development of systems not seen in surface-dwelling organisms) or regressive (loss of systems occurring in surface-dwelling organisms) adaptations [7,8,9]. Along with general troglomorphisms widespread across cave-adapted taxa (larger body size, reduced eyes and pigment, elongated appendages and claws), specific modifications appear in cave Collembola: shorter, thinner and pointed tenent hairs on tibiotarsi compared to their edaphic counterparts, reduction and basal shift of the inner teeth on the claw, and hypertrophy and/or multiplication of antennal sensilla [8,10].
Subfamily Onychiurinae as typically edaphic Collembola, are unpigmented and eyeless forms in which other troglomorphisms rarely develop. The exceptions are the three troglobiotic onychiurids, namely Ongulonychiurus colpus (Thibaud & Massoud, 1986) discovered in the deep abysses of Picos d´Europa (Spain); a new Ongulonychiurus from a deep cave in the Biokovo Mts (Croatia) [11]; and Troglaphorura gladiator Vargovitsh, 2019 [12], recently described from the Snezhnaya Abyss in the Caucasus Mts (Georgia). All species are characteristic with highly troglomorphic traits.
The genus Deuteraphorura (Poduromorpha, Onychiuridae) includes 83 species [13], 49 of which have been described from caves. Along others, the number and distribution of pseudocelli and parapseudocelli (hardly visible soft cuticular structures) over the body, play an important role in morphological identification of species of the family Onychiuridae. However, most onychiurids are well-characterized by abovementioned features, high variability of morphological characteristics and frequently also their left-right asymmetry (e.g. body chaetotaxy, number of pseudocelli and parapseudocelli on the body segments) is often in some species of the genera Protaphorura and Deuteraphorura.
Integrative taxonomy is a synthetic approach combining evidence from independent data sources (e.g. morphology, DNA, ecology, geography) in order to reveal biological features relevant for species delimitation [14,15]. This integrative approach has been successfully applied to collembolan taxa, where the colour pattern was previously the only helpful characteristic used in their identification, e.g. [16,17,18,19]. Many studies have reported congruent results when a combination of morphological and molecular approaches was used in species delimitation, e.g. [20,21,22,23,24]. In contrast, low morphological but clear genetic differences [17,25,26,27,28], or clear morphological but low genetic differences [29], have also been documented in Collembola. In this regard, the subfamily Onychiurinae seems to represent a serious challenge due to the uniform colouration, similar ecology and habitat preferences of its representatives and the often unstable and/or asymmetric chaetotaxy and pseudocellar pattern over the body, especially in genera Protaphorura and Deuteraphorura (e.g. [30]).
We assumed that the divergence of cave Deuteraphorura may be higher than previously expected, including older phyletic relicts in the territory of the Western Carpathians.
This contribution aimed at 1) evaluating the congruence between morphological traits, DNA barcode and geographical data for species delimitation of cave populations of the genus Deuteraphorura occupying the Western Carpathian caves; 2) clarifying the species status of Deuteraphorura kratochvili (Nosek, 1963), the most widespread troglobiotic species of the genus in the study territory; and 3) describing the new remarkably troglomorphic species of the genus Deuteraphorura discovered in the Muránska planina Plateau karst area. The level of troglomorphy in subterranean representatives of the genus Deuteraphorura is discussed regarding their distribution over the northernmost boundary of the range of terrestrial troglobionts in Europe.

Material and methods
Individuals of the genus Deuteraphorura were collected by hand in caves of the Western Carpathians, Slovakia (see futher text), on cave sediment, speleothems, bat guano, rotten wood and the surface of the water pools and immediately fixed in 96% ethyl alcohol. For morphological study, the specimens were separately mounted on permanent slides in Swann medium (Liquido de Swann) modified after Rusek [31] and studied through a phase-contrast Leica DM 2500 microscope equipped with DIL optics (differential interference contrast), a measuring eyepiece (micrometric ocular) and a drawing arm. The images were taken with an Olympus 5.1 megapixel camera and edited using Adobe Photoshop CS6. Chaetotaxy of the tibiotarsus is presented after Deharveng [32], and chaetotaxy of the labium after Fjellberg [33]. Claw length was measured along the internal side as the distance from the distal margin of the pretarsus to the top of claw. Claw width (lateral view) was measured as the width of the pretarsus in conjunction with tibiotarsus.
Five populations of Deuteraphorura were included in the molecular analyses; four populations of D. kratochvili from the following caves: Demänovská Ice Cave (Low Tatras), Stará Brzotínska Cave (Slovak Karst), Drienovská Cave (Slovak Karst) and Ochtinská aragonitová Cave (Revúcka vrchovina Mts), and one population of D. muranensis sp.nov. from Jelenia priepasť Abyss (Muránska planina Plateau) (Fig 1). Deuteraphorura kratochvili populations were examined for the purpose of the species redescription [30]. In this study, populations were defined as individuals collected in separate caves. The selected caves are located in four orographic units in Slovakia, with a minimal air distance of 12 km between the closest caves (Stará Brzotínska Cave-Ochtinská aragonitová Cave) and a maximum air distance of 110 km between the most distant caves (Drienovská Cave-Demänovská Ice Cave).
The sampling of Collembola material was carried out with the permission of the Ministry of the Environment of the Slovak Republic (No. 2661/2017-6.3).

Molecular data analysis
To prevent contamination, all DNA laboratory work was conducted under sterile conditions with the use of barrier tips. Total DNA was extracted using the Qiagen DNeasy Blood & Tissue Kit according to the modified manufacturer's protocol (see [34]). Polymerase chain reaction (PCR) [35] was carried out using a 12.5 μL of reaction volume consisting of 1 μL of template DNA (not quantified), 10 × PCR Buffer (TopBio s.r.o., CZ), 12.5 mM of dNTP mix, 5 μM of each primer and 0.125 units of Taq polymerase (TopBio s.r.o., CZ) on a GenePro (Bioer Co., Ltd, China) thermal cycler. A fragment (670 bp) of the COI gene was amplified using the universal primers LCO1490 (5'-ggt caacaa atcataaagatattg g-3') and HCO2198 (5'-taa act gggtga ccaaaaaat ca-3'; [36]). The thermal cycling conditions were as follows: 94˚C for 1 min followed by 35 cycles (94˚C for 30 sec, 47˚C for 35 sec and 72˚C for 50 sec), followed by 1 min 30 sec at 72˚C. After verification via agarose gel electrophoresis, the reaction products were purified using Exo I/FastAP (Thermo Fisher Scientific Inc.). The sequencing of the purified products was performed using LCO1490 at the SEQme s.r.o. in Dobris, Czech Republic, by the Sanger method. In cases when the primer failed to produce high quality chromatograms, reverse primer sequencing was employed. Sequences were manually edited and trimmed of unreadable short stretches (ca 30 bp at the 5' and 3' ends) with Bioedit v.7 [37]. Since none of them contained stop codons or indels in ORF, all were considered to be true mitochondrial and not nuclear copies. All the sequences were verified as being consistent with Onychiuridae congeners using the GenBank BLASTn search (the Mega Blast algorithm with the default setting). Sequences were aligned with the MEGA 7.0.26 [38] software by the Muscle (Codons) algorithm using the Invertebrate Mitochondrial GeneCode and default parameters. Standard DNA barcoding distance analysis was conducted using the Tamura-3 parameter method [39] recommended by model selection. A neighbour-joining tree [40] with the Tamura-3 parameter method [39] was constructed, and the robustness of the tree nodes was assessed by bootstrap analysis with 1000 replications. Deuteraphorura inermis (Tullberg, 1869) (GenBank code KY231117.1) was used as the outgroup.
Both barcoding gap-and evolutionary models were applied for the COI marker. Automatic Barcode Gap Discovery (ABGD) clustered the sequences into candidate species based on pairwise distances by detecting barcoding gaps for the COI marker [41]. Prior intraspecific divergence varied from 0.001 (Pmin, a single nucleotide difference) to 0.14 (Pmax, threshold value) applied by Porco [42]. The relative gap width was set to 1, with 20 recursive steps, 40 bids for a graphic histogram of distances, a K2P model [43] for distance calculation and other parameters as the default.
The Poisson tree processes (PTP) model based on significant differences in the number of substitutions was performed using on-line software [44]. Identical sequences were removed to avoid incorrect likelihood calculations, and an unrooted NJ tree was constructed using the IQ-TREE web server with the default parameters (http://iqtree.cibiv.univie.ac.at/). Correlation between the genetic and geographical distances of the populations was evaluated by means of the Mantel test (999 permutations) using the GenAlEx 6.5 program. Haplotype diversity (h) was calculated using DnaSP 5 [45]; subsequently, a haplotype network for the Deuteraphorura populations was constructed using the program Network 5.0.0.3 (website fluxus-engineering.com) (results not shown). The presence/absence of shared haplotypes was used to verify gene flow among the populations.
All the new sequences are available in GenBank (accession numbers for D. muranensis: MN727813 -MN727820 and for D. kratochvili: MN727821 -MN727844).

Nomenclatural acts
The electronic edition of this article conforms to the requirements of the amended International Code of Zoological Nomenclature, and hence the new names contained herein are available under that Code from the electronic edition of this article. This published work, and the nomenclatural acts it contains, have been registered in ZooBank, the online registration system for the ICZN. The ZooBank LSIDs (Life Science Identifiers) can be resolved and the associated information viewed through any standard web browser by appending the LSID to the prefix "http://zoobank.org/". The LSID for this publication is: urn:lsid:zoobank.org:act: AD211C26-DCCA-4042-A677-EBB1CBC528D6. The electronic edition of this work was published in a journal with an ISSN, and has been archived and is available from the following digital repositories: PubMed Central, LOCKSS.

Results
Molecular analysis of the mitochondrial cytochrome oxidase I (COI) fragment was carried out in five Deuteraphorura populations, the morphology and chaetotaxy of which were examined for the purpose of D. kratochvili redescription [30].
Genetic diversity was calculated between and within species and populations, respectively (Table 1). In D. kratochvili, intraspecific divergence was 11%, while it approximated zero within particular populations, except for Demänovská Ice Cave, with 3.5%. A significant positive correlation between the geographical and genetic distances was confirmed by the Mantel test (R 2 = 0.53; p = 0.001). In D.muranensis from Jelenia priepasť Abyss, intrapopulation divergence reached 4.8%.
Among 32 D. kratochvili sequences 13 haplotypes were detected, and none of them was shared among the populations. Each population was represented by 2-4 haplotypes, with the exception of Drienovská Cave, with only one haplotype.
ABGD detected prior intraspecific divergences (P) between 0.001 and 0.049, while 10 and 7 MOTUs were recognized in the initial and recursive partitions, respectively. Barcoding gaps at K2P distance were observed between 0.01-0.055, 0.065-0.075 and 0.095-0.11 (Fig 2). Initial and recursive partitions were congruent at about P = 0.03, while the number of MOTUs was 7. ABGD species delimitation corresponded to particular populations; however, populations from the Jelenia priepasť Abyss and Demänovská Ice Cave were further divided into two partitions each.
PTP delimitation based on haplotype sequences identified 7 species/MOTUs (support 0.88-1.0), thus fitting the results of the AGBD analysis. Species delimitation approaches were summarized on an NJ-tree (Fig 3).
Biology. Deuteraphorura muranensis sp. nov was collected in two caves of the Muránska planina Plateau, central Slovakia, almost exclusively on the water surface of sinter pools, but also on bat guano and organic baits in the deeper parts of the caves. The new species possesses distinctly elongated claws, one of the troglomorphic adaptations to cave life. These caves are type localities of another troglobiotic collembolan, an undescribed species of the genus Pseudosinella, likewise characteristic with conspicuously long claws.
Etymology. The species name is derived from the name of the closest village Muráň, with an old castle, and the karst region Muránska planina Plateau, where the caves inhabited by this species are situated.

Remarks on the new species
D. muranensis sp. nov. represents a well-defined species morphologically, genetically and geographically. It is the only representative of the genus Deuteraphorura with 4 pso on the hind Subterranean Deuteraphorura -Troglomorphy at the northern distributional limit in Europe PLOS ONE | https://doi.org/10.1371/journal.pone.0226966 January 15, 2020 margin of the head and the absence of pso on Th. I. Several species with 4 pso on the hind head margin have been described from caves, however, with pso on Th. I tergum. By the dorsal abdominal pso formula, the new species is most similar to congeners occupying caves in Spain ( Table 2). Parapseudocellar and pseudoporal formulas have been used as determining characteristics in several Onychiuridae taxa, e.g. [22,46,47]; however, these characteristics are often only barely visible, e.g. [48], as in the case of D. muranensis sp. nov. In addition, the new species shows a high level of left/right asymmetry in dorsal chaetotaxy and the distribution of the ventral abdominal pso (Fig 4B).
The claws of the new species (evaluated as average length/width ratio 3.3) are distinctly longer compared to other cave-dwelling congeners from the Western Carpathians-D. schoenviszkyi (Loksa, 1967) and D. kratochvili-with average length/width ratio 2.6 [30]. Subterranean Deuteraphorura -Troglomorphy at the northern distributional limit in Europe

Species delimitation and cryptic diversity in cave Deuteraphorura in the Western Carpathians
D. kratochvili represents the most widespread cave species, with a distribution range covering the Western Carpathian Mts, where the existence of two morphological forms, D. kratochvili Subterranean Deuteraphorura -Troglomorphy at the northern distributional limit in Europe and D. cf. kratochvili, with distribution limited to the adjacent karst areas was assumed [49,30]. In our former study [30], morphological variability was found in all Deuteraphorura populations, but the studied characteristics were not population-specific; thus, we considered them all as belonging to the same species, D. kratochvili. However, there is another obligate cave congener in the Western Carpathian caves that has not been successfully collected for more than ten years. Deuteraphorura schoenviszkyi (Loksa, 1967) is endemic to the Slovak and Aggtelek Karst and represents a morphologically well-defined species with distinctly elongated claws [49]. Unfortunately, no molecular analysis has yet been done, since occurrence is very rare in the caves of the area.
An integrative approach has been repeatedly highlighted as the best possible future of taxonomy, employing morphological and molecular attributes as well as attributes of the extended phenotype, behaviour and ecology [15,19,22,23,24]. We applied morphological, molecular and geographical approaches to resolve the taxonomic status of the cave Deuteraphorura representatives occurring in the Western Carpathians. Four D. kratochvili populations and one population of the new species were used in this integrative study. Although D. kratochvili occupies most Western Carpathian caves, the sampling of a sufficient number of individuals for both molecular and morphology analyses was highly demanding in most of them due to difficult access to individual cave sites and the rather dispersed distribution of the individuals in caves in general.
In the present study, groups of specimens defined by ABGD and PTP analyses respectively, corresponded to particular caves (populations). However, an additional separation was indicated within two populations, namely Jelenia priepasť Abyss and Demänovská Ice Cave. Similarly, relatively large intraspecific distances (4% or greater) have been reported within several collembolan species, e.g. [16,17,24,26,27]. Interpopulation divergence in D. kratochvili was higher than between the populations of P. janosik Weiner, 1990, from several Western Carpathian caves [34], suggesting longer isolation and evolution of D. kratochvili, thus probably representing an "old troglobiont", i.e. a troglobiont that is a descendant of an old phyletic lineage [4]. Subterranean Deuteraphorura -Troglomorphy at the northern distributional limit in Europe Since reliable species boundaries are difficult to delimit by genetic distance alone, and morphological differences are indistinct in D. kratochvili, we used a geographical criterion, which proved to be an efficient delimitation method for onychiurids [22]. The results of both molecular analyses and the geographical distance between localities suggested the existence of genetically isolated forms of D. kratochvili with positively correlated geographical and genetic distances. A similar correlation has been documented in several cave collembolan populations [21,27,34] and coincides with the intensity of gene flow. In Collembola, gene flow was confirmed across large geographical distances, but in contrast, also a minimal sharing of haplotypes, and thus reproductive isolation may be demonstrated between adjacent or even sympatric species that are indistinct morphologically [17,50]. In geographically distant and physically isolated populations, some intra-specific genetic distance has been reported, even though only 3-4% [27,34]. The intraspecific distance in D. kratochvili (11%) was higher than minimal interspecific divergence in Chinese Protaphorura documented in two species separated by more than 1000 km [22]. Based on these results, it is clear that D. kratochvili represents a complex of several species. The absence of reliable morphological characteristics and/ or the high level of their variability were the reasons why traditional taxonomic methods failed to reveal these cryptic species.
An increasing trend in the number of collembolan species is obvious, since molecular methods have become an integral part of taxonomy. The diversity of Collembola is greatly underestimated owing to the prevalence of cryptic species [51]. Whereas 14% was adopted as the threshold for species delimitation [42], it has not been strictly defined. Ubiquist species with large distribution ranges, including several more or less geographically isolated populations, have a high intraspecific (interlineage) variability reaching more than 20%, and particular lineages may represent cryptic species [25,26]. Along with cryptic (morphologically uniform) species with surprisingly high interspecific distances, morphologically, geographically and/or ecologically well-defined species with low interspecific distance have also been observed [29]. In some species previously identified as cryptic taxa by molecular distance, additional morphological differences were found after subsequent more thorough, detailed morphological re-examination of individuals [19]. Finally, if only genetic distance is applied for species delimitation, some species could be neglected, as low divergence between them would be assigned as intraspecific variability.
Generally, without thorough study of morphology and the use of geographical distance data species delimitation may be biased. Due to the lack of any morphological differences and minor interspecific (JPA populations = D. muranensis vs. OAC populations = D. kratochvili), major intraspecific (D. kratochvili) molecular divergence as well as interspecific distance smaller than intraspecific, we consider D. kratochvili to be a complex of cryptic species in which long-term isolation in individual karst systems likely led to gene-flow cut-off.

Troglomorphy of cave Collembola at the northernmost distributional border in Europe
In Europe, troglobiont diversity increases towards the southern mountains (Pyrenees, Alps, Dinaric Mts), reaching a maximum in so-called "hot-spot" areas [3]. The present study contributes to the recent findings that the Western Carpathians are the northernmost limit of the occurrence of the troglomorphic taxa in Europe, inhabited by several obligate cave Collembola species [4,30,34,49,52,53]. The distribution pattern of "hot-spot" areas is basically associated with the productivity of surface ecosystems, spatial heterogeneity and historical conditions covering the effect of the past climatic and paleogeographic events [54]. The presence of such a hotspot reflects not only the local diversity of troglobionts but may also be associated with the high level of troglomorphy in these obligate subterranean forms. A heterogeneous community of highly troglomorphic collembolan taxa of the families Tomoceridae, Entomobryidae, Arrhopalitidae and Neelidae may occupy the caves of a "hot-spot" area [55,56,57,58,59], while the level of troglomorphy is only moderate in troglobionts at the borders of the distributional range [4,49,53]. Morphological adaptations seem to be the most distinct in Entomobryidae, while basically inconspicuous in Onychiurinae, a subfamily which is "uniform" in general morphology and ecology. Moreover, regressive morphological changes, such as loss of pigmentation and eyes characteristic in edaphic taxa, could not be assigned as troglomorphisms. One morphological manifestation of troglomorphy in Collembola is elongation of body appendages [7,60]. In Onychiuridae, elongation of morphological structures in troglobiotic representatives is most distinct in the foot complex. Generally, the ratios claw length/width, claw length/Tita length or claw length/body length have been used to estimate the level of troglomorphy, e.g. [12,30,34,59]. Despite the fact that the subfamily Onychiurinae is extremely edaphic with several eutroglophile forms, only three highly troglomorphic species are known: two species of Ongulonychiurus and Troglaphorura gladiator. Whereas both Ongulonychiurus are distributed within hot-spot areas (Picos d´Europa, Spain and Biokovo Mts, Croatia), Troglaphorura was discovered in the Western Caucasus, representing indeed an important "hotspot" area of the subterranean fauna with a great potential for further important discoveries of the highly specialized subterranean invertebrates.
Endemic distribution of particular species in a "hot-spot" area is not indeed a "guarantee" of a high level of troglomorphy in the local taxa, as shown in several Deuteraphorura species from caves of Navarra (Spain) [61]. Absolonia gigantea (Absolon, 1901), occupying caves of the Balkan Peninsula, displays the same level of troglomorphy as Protaphorura janosik Weiner, 1990 [34] and P. cykini [62] from the Western Carpathians in Slovakia and the Baikal territory in Siberia, respectively, achieving a larger body size compared to other congeners. Moreover, recent discoveries of new troglomorphic taxa in previously "biospeleologically neglected" localities [12] document that the area with a high degree of troglomorphy in Collembola is larger in Europe than was defined by Culver et al. [3]. On the other hand, discoveries of troglomorphic taxa and the diversity of troglobionts in this group also depend, among other things, on the number of caves in the individual region, the local rate of biospeleological exploration and the sampling effort.
The highly morphologically adapted subterranean species of the genus Deuteraphorura described in this paper was found in two caves of the Muránska planina Plateau, a small karst area in central part of the Western Carpathians (213.2 km 2 in area) with~500 caves. We are not able to describe the possible governing factors that led to the evolution and endemic distribution of this highly troglomorphic species in this area. In fact, Christiansen & Culver [63] noted the relation between the level of troglomorphy and distribution-range size in Collembola, i.e. that increasing troglomorphy decreases their ability to disperse. Similarly, D. schoenviszkyi has endemic distribution and displays a higher level of troglomorphy than D. kratochvili, which has greater distribution range [30,49]. In addition, increasing troglomorphy indicates early (pre-glacial) times of first cave colonization [63]. Of course, this remains to be tested using a phylogenetic approach based on a broader molecular dataset originating from multiple Deuteraphorura populations across the Western Carpathian caves.