Research Article
Print
Research Article
Comparative phylogeography of two troglobitic Coleoptera (Leiodidae, Leptodirini) species from Romania based on mitochondrial DNA
expand article infoRuxandra Năstase-Bucur‡§, Giuliana Allegrucci|, Valerio Ketmaier, Ionuţ Cornel Mirea§#, Oana Teodora Moldovan‡§
‡ Emil Racovitza Institute of Speleology, Cluj-Napoca, Romania
§ Romanian Institute of Science and Technology, Cluj-Napoca, Romania
| University of Rome Tor Vergata, Rome, Italy
¶ Unaffiliated, Berlin, Germany
# Emil Racovita Institute of Speleology, Bucharest, Romania

Corresponding author: Ruxandra Năstase-Bucur ( ruxandra.nastase.bucur@gmail.com )

Corresponding author: Giuliana Allegrucci ( allegrucci@uniroma2.it )

Academic editor: Arnaud Faille
© 2022 Ruxandra Năstase-Bucur, Giuliana Allegrucci, Valerio Ketmaier, Ionuţ Cornel Mirea, Oana Teodora Moldovan.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation: Năstase-Bucur R, Allegrucci G, Ketmaier V, Mirea IC, Moldovan OT (2022) Comparative phylogeography of two troglobitic Coleoptera (Leiodidae, Leptodirini) species from Romania based on mitochondrial DNA. Subterranean Biology 42: 61-78. https://doi.org/10.3897/subtbiol.42.73524
ZooBank: urn:lsid:zoobank.org:pub:731EE7BC-58F1-45D7-BBB8-84321F876D97
Open Access

Abstract

About 50 species of cave-obligate Leptodirini (Leiodidae) beetles have been described so far in Romania, most of them populating caves in the Apuseni Mountains (north-western Romania) and the Southern Carpathians. In this contribution, we present the first molecular phylogeographic study of the two troglobiotic Pholeuon species from the Apuseni Mountains. The two species are Pholeuon (s.str.) leptodirum and Pholeuon (Parapholeuon) gracile, endemic to Bihorului Mountains and Pădurea Craiului Mountains, respectively. To examine the genetic divergence within and between the two species we sequenced 571 bp of the mitochondrial COI gene in a total of 145 specimens, 56 specimens of the first species (collected in five caves) and 89 of the second species (collected in eight caves) across their geographic ranges. We found very low genetic variation, four haplotypes in P. leptodirum and seven haplotypes in P. gracile, and a maximum of 0.7% and 0.9% intraspecific divergence, respectively. However, a significant genetic divergence of 6.55% was found between species. The results are consistent with previous definitions of the two species based on morphological characters, while caution should be taken in considering attributions to different subspecies. Our research contributes to the phylogeographic information of troglobitic beetles, providing a solid basis for future comparison with other terrestrial or aquatic cave adapted species.

Keywords

Carpathians, cave beetles, cytochrome oxidase gene I, Pholeuon, population genetics

Introduction

A detailed knowledge of the biology and ecology of the species and their genetic structure at the population level plays a crucial role in minimizing the effect of loss of biodiversity. Important target of global conservation efforts are endemic species. These, with their specific climatic and environmental requirements and generally limited dispersal capacity, are particularly vulnerable to extinction (Myers et al. 2000; Lamoreux et al 2006). Obligate subterranean species, cave-adapted so called troglobionts, are particularly vulnerable to the pollution produced on the surface, which percolates soil and layers of limestone, contaminating the subterranean habitats (Wood and Perkins 2002; Manenti et al. 2021). Subterranean environments are inhabited by a specialized fauna living in relatively stable conditions that creates a unique biological laboratory where evolutionary and ecological processes can be studied in situ (Mammola 2019). Subterranean terrestrial habitats are characterized by absence of light, saturated air humidity, constant temperatures and often scarce food resources (Howarth et al. 2008). The ecological balance is, therefore, fragile and any disturbance can potentially cause alteration of the natural conditions, fragmentation of the habitat and populations, with the possible extinction of the troglobitic species.

Romania hosts many unique karst landscapes and caves, and several types of endemism are found in its subterranean fauna. The occurrence and distribution of the Romanian cave fauna can be to a larger extent explained by paleogeography and ecology of the group (Moldovan 2008). Among them, the Leptodirini (Coleoptera, Leiodidae) species include strictly cave-adapted beetles and represent a conspicuous group of species endemic to one or a few caves in a karst area. Indeed, eight Leptodirini genera with 49 endemic species inhabit Romanian caves (Moldovan et al. 2020). Their limited distribution in a fragmented karst landscape make Romanian Leptodirini a model group for studying speciation and historical biogeographic processes both at fine and large scales. Molecular studies on Leptodirini species are still very limited in number, and the most comprehensive study was undertaken by Ribera et al. (2010); this comprised 57 Leptodirini species from all the major lineages distributed in the Western Mediterranean area, including two Romanian species. The molecular clock approximation suggested that the main Western Mediterranean lineages originated in Early-Mid Oligocene and that the ancestral species were already present in the geographical areas in which they are found today (Ribera et al. 2010).

Romanian Leptodirini group is considered to derive from the Dinaric ancestors, before the separation of the Carpathian Mountains from the Dinarides (Jeannel 1931; Decu and Negrea 1969; Moldovan and Rajka 2007). The Romanian Carpathians are divided into three geographical units: The Apuseni Mountains (North-West), the Eastern Carpathians, and the Southern Carpathians. The Apuseni Mountains are characterized by an amazing variety of relief micro forms, both above and underground, and a considerable number of caves (3,960 caves, with a cave density of 3.5 caves/km2) (Cocean 2000). The unique features of the fragmented karst in this region, interrelated with the underground hydrological network, could present a great opportunity for the dispersal of the species. Whereas the complex geological settings, from local to regional context, potentially function as geographical barriers for the dispersal of subterranean fauna (Moldovan 2008).

Romanian Leptodirini were studied, so far, for various aspects: ecological, taxonomic, either by classical taxonomy (the revision of genus Drimeotus; Moldovan 2000) and by morphometrics (Racoviță 1996; Racoviță 1998–1999; Racoviță 2010, 2011), historical biogeography (Moldovan and Rajka 2007), cytogenetic (Buzilă and Marec 2000), molecular phylogeography (Bucur et al. 2003). Studies on molecular phylogeography of Leptodirini were carried out on the phyletic series of Drimeotus with three genera, Drimeotus, Pholeuon, and Protopholeuon, all three endemics to Apuseni Mountains (Bucur et al. 2003). Drimeotus and Pholeuon include three and two subgenera, respectively, while Protopholeuon is a monospecific genus. Species belonging to Pholeuon and Protopholeuon are considered more troglomorphic (sensu Christiansen 2012) (longer appendages, slender body) than Drimeotus (Moldovan et al. 2007). The study of Ribera et al. (2010), including representatives from Romanian Carpathians, established the monophyly of the phyletic series Drimeotus. On the contrary, the subgenera belonging to Drimeotus and Pholeuon genera were not monophyletic (Bucur et al. 2003), offering multiple independent colonization events by surface ancestors as a possible explanation for their actual distribution.

In the present study we carried out a comparative phylogeographic analysis of the two cave adapted Leptodirini subgenera from the Apuseni Mountains. In particular, we considered two mountain ranges, Pădurea Craiului and Bihorului, with the species belonging to their corresponding endemic subgenera: Pholeuon (Parapholeuon) and Pholeuon (s. str). Until present, three species have been described from each subgenus, Parapholeuon with P. gracile, P. moczaryi, and P. angustiventre and Pholeoun s. str. with P. angusticolle, P. knirschi and P. leptodirum. For this study, one species from each subgenus, P. gracile and P. leptodirum, both including several populations/subspecies, were considered. The analyses included all three described subspecies of P. gracile while only five of eleven described subspecies of P. leptodirum were included (Racoviță 2011). The description of subspecies is rather controversial as it is solely based on morphometric characters, and their taxonomical status reaches beyond the scope of this paper.

The main aims of the present study were: i) to test the congruence between the subspecies identified on the basis of morphometric measurements and the putative molecular species delimited using DNA barcoding (mitochondrial cytochrome oxidase subunit I, COI); and ii) to investigate the degree of genetic differentiation both within and between the two cave-adapted species, and to analyse it according to the geographic distribution of populations.

Materials and methods

All analyzed populations are endemic in two massifs of the Apuseni Mountains, Bihorului and Pădurea Craiului, and have limited distribution in one or few caves only (Fig. 1 and Table 1).

Figure 1. 

Sampling sites of P. (Parapholeuon) gracile in Pădurea Craiului Mountains (blue) and P. (s. str.) leptodirum in Bihorului Mountains (red). Codes refer to the caves’ name, as indicated in Table 1.

Table 1.
Download as
CSV
XLSX

Leptodirini subspecies (as proposed by Racoviță 2011) included in the study with the codes used throughout the paper, the hydro-karstic basins, altitude and slope where the caves are located. Sample size (N), haplotypes and number of specimens sharing each haplotype are also indicated.

Cave Name Code River valley Basin Subspecies Altitude (m a.s.l.) Geographic slope N Haplotype (no. specimens)
Pădurea Craiului Mountains – P. (Parapholeuon) gracile
Cubleș CUB Vida Holod P. g. gracile 440 right 11 H1 (9), H2 (2)
Vizu II VIZ Vida Holod P. g. bokorianum 350 left 12 H3 (5), H6 (6), H5 (1)
Tocoș TOC Runcșor Roșia P. g. chappuisi 585 right 11 H1 (8), H2 (2), H7 (1)
Întorsuri INT Runcșor Roșia P. g. chappuisi 575 right 11 H1 (9), H2 (2)
Ciur Ponor CPO Albioara Roșia P. g. chappuisi 510 left 10 H3 (10)
Doboș DOB Albioara Roșia P. g. chappuisi 465 left 11 H3 (11)
Vălău VAL Albioara Roșia P. g. chappuisi 355 right 12 H3 (9), H4 (3)
Gruieț GRU Șteazelor Roșia P. g. chappuisi 300 left 11 H1 (9), H2 (2)
Bihorului Mountains – P. (s.str.) leptodirum
Coliboaia COL Sighiștel Sighiștel P. l. jeanneli 560 right 11 H8 (11)
Măgura MAG Sighiștel Sighiștel P. l. hazayi 550 right 12 H8 (12)
Corbasca COR Sighiștel Sighiștel P. l. moldovani 500 left 11 H8 (6), H9 (5)
Fânațe FAN Bulzului Crișul Băița P. l. leptodirum 560 right 10 H8 (8), H11 (2)
Secătura SEC Bulzului Crișul Băița P. l. problematicus 1080 right 12 H10 (12)

Sampling

Specimens of P. gracile, belonging to the three described subspecies (P. gracile s. str., P. g. chappuisi, and P. g. bokorianum), were collected from eight caves located in four different valleys of Pădurea Craiului Mountains (Table 1, Fig. 1). Five populations of P. leptodirum (P. leptodirum s. str., P. l. jeannely, P. l. hazayi, P. l. moldovani, P. l. leptodirum and P. l. problematicus) were collected in two different valleys in Bihorului Mountains (Table 1, Fig. 1).

In each cave, between 10 and 12 individuals were collected for the analysis. The total number of specimens was 145 (89 individuals for the first species and 56 for the second one). Specimens were preserved in 95% ethanol until DNA extraction was processed.

DNA extraction, PCR and sequencing

Total genomic DNA was extracted from the entire specimens using the DNeasy Blood and Tissue Kit (Qiagen), following the producer’s protocol. A fragment of 571 base pairs (bp) of the mitochondrial Cytochrome Oxidase I (COI) gene was amplified using LCO1490 and HCO2198 primers (Folmer et al. 1994).

Double stranded amplifications were performed in a 50 µl reaction volume containing buffer, 5 µl dNTP’s 10 mM, 0.5 µl primer 10 mM (each primer), 0.4 µl TAQ polymerase (5U/µl) and 38.6 µl purified water. Each PCR cycle (of a total of 30 cycles) consisted of a denaturation step at 94 °C for 1 min, annealing at 50 °C for 1 min and extension at 72 °C for 7 min. PCR products were purified following the manufacturer’s protocol for the PCR-Nucleospin Gel and PCR Clean-Up (Macherey-Nagel). Both strands were sequenced on an automated sequencer.

Genetic analyses

Sequences were aligned and edited with BioEdit (v. 7.2), the number of transitions and transversions were analysed with DNAsp (v 5.10.1) (Librado and Rozas 2009). We used MEGA7 (v.7.0) (Kumar et al. 2016) to analyse interspecific haplotype diversity. Genetic structure between populations of the two species or within populations of the same species was analysed with F-statistics using Arlequin (Excoffier et al. 2010) by calculating the following parameters: haplotype diversity (h), the absolute haplotype frequencies, and the nucleotide diversity (πn). The minimum spanning network was built using PopArt software (Leigh and Bryant 2015).

Mean pairwise intra- and interspecific distances were determined using MEGA (Tamura et al. 2013). Analysis was conducted using uncorrelated p-distance. The analysis involved 145 nucleotide sequences. All codon positions were included.

Using PAST software (Hammer et al. 2001) a Mantel test (Mantel 1967) with 5,000 simulations was carried out to test for an isolation-by-distance (IBD) signature (a positive correlation between geographic and genetic distances (Wright 1943; Slatkin 1993).

The hierarchical distribution of genetic variation was characterized using analysis of molecular variance (AMOVA). This method apportions genetic variation within and among groups, estimating Φ-statistics (Weir and Cockerham 1984; Excoffier et al. 1992; Weir 1996) that are analogous to Wright’s hierarchical fixation indices (FST) under the island model of gene flow (Wright 1951). Three-level AMOVA was conducted in ARLEQUIN 3.5.1.2 (Excoffier et al. 1992, 2005) using an FST-like estimator (Fixation Index). For each of the two considered species, AMOVA was run three times considering different groups of populations on the basis of the geological and geographic characteristics. In particular, samples were partitioned by the river basin (Table 1), populations within geographic regions and inside each population. The tests included permutation of inferred haplotypes among groups (FCT); individual haplotypes among populations but within group (FSC); inferred haplotypes among populations (FST).

In order to test for the monophyly of the two Leptodirini species we carried out phylogenetic analysis within a Bayesian framework. J model test (Dariba et al. 2012) was used to perform a hierarchical likelihood ratio test and calculate approximate Akaike Information Criterion (AIC) values of the nucleotide substitution models.

Phylogenetic analysis was performed using Bayesian inferences as implemented by the software MrBayes 3.2.7 (Ronquist et al. 2012). Two simultaneous searches, comprising four Markov chains (MCMC) each and starting from a randomly chosen tree were run for 1,000,000 generations and sampled every 100 generations.

Convergence on a common phylogenetic topology by separate Bayesian searches was checked using Tracer 1.7 (Rambaut et al. 2018). The effective sample size (ESS) of all parameters showed values above 1,000 (values much higher than the threshold of statistical significance, ESS > 200) in both simultaneous searches, indicating that MCMC had converged. Out of 20,000 trees, the first 1,000 were discarded as burn-in, and posterior probabilities (PP) were calculated from post-burn-in trees. The tree was rooted with a species of Leptodirini from Sardinia, Ovobathysciola sp.; we were provided with one specimen.

Results

Sequences were obtained from a total of 145 individuals, 89 for P. gracile and 56 for P. leptodirum (Table 1).

The alignment consisted of 571 bp and defined 46 variable sites, of which 44 were parsimony informative. The nucleotide diversity among all sequences was πn = 0.033. The sequences are deposited in Genbank with the Accession Numbers OL457148–OL457159 (for H1–H11 and Ovobathysciola, respectively).

Intraspecific variability

We identified seven haplotypes for P. gracile (numbered H1–H7; Table 1), separated from each other by a maximum of eight mutations, whose geographic distribution is represented in Fig. 2. The haplotypes were separated in two haplogroups according to the genetic divergence between them. One group consisted of specimens of Doboş (DOB), Valău (VAL), Ciur-Ponor (CPO) from Albioara Valley/Sohodol basin, and Vizu II (VIZ) from Vida Valley/Holod basin. The second group was represented by Cubleş (CUB) from Vida Valley/Holod basin, Tocoș (TOC)/Sohodol basin, Ȋntorsuri (INT) from Runcșor Valley/Runcșor basin, and Gruieţ (GRU) from Șteazelor Valley/Roșia basin.

Figure 2. 

Geographic distribution of haplotypes in P. gracile and P. leptodirum. Individuals of DOB and CPO caves are fixed for the haplotype H3 in P. gracile; individuals of SEC cave are fixed for the unique haplotype H10 in P. leptodirum (abbreviations for the names of caves as in Table 1).

The most widespread haplotypes were H1 and H3, each showing a frequency of 40% in all the analyzed specimens. In particular, individuals from Ciur-Ponor and Doboş caves were fixed for H3 haplotype that was also identified in two other caves, Valău (75%) and Vizu II (42%). Haplotype H1 was shared by the individuals from Cubleş (82%), Gruieţ (82%), Ȋntorsuri (82%), and Tocoş (73%) caves. Haplotypes H5 and H7 were identified in a single specimen, each from Vizu II and Tocoș caves, respectively.

Haplotypes H6 and H4, present with a frequency of 50% and 25%, appeared to be exclusive of Vizu II and Vălău caves, respectively. Haplotype H2 was spread in four caves, Cubleş, Gruieţ, Ȋntorsuri, and Tocoş, with a frequency of 18% in each case. The haplotype diversity for P. gracile was Hd = 0.684 and nucleotide diversity was πn = 0.003.

The haplotype network for P. gracile is illustrated in Fig. 3a. The identified haplotypes are divided in two clusters: one comprising haplotypes H1, H2 and H7, and the other haplotypes H3–H6. The two clusters are differentiated by three mutations. Within each cluster, the haplotypes are separated by only one mutation.

Figure 3. 

Minimum spanning networks a for P. gracile with haplotypes H1–H7 and b for P. leptodirum with haplotypes H8–H11. Multiple mutational steps between haplotypes H1–H3 and H8–H9 could be either un-sampled haplotypes or extinct ones. Haplotype numbers are as in Table 1.

FST values for P. gracile are presented in Table 2, in which 19 of 28 inter-populational comparisons provided significance at 0.05 level of probability. All negative values indicated a lack of genetic differentiation between the respective populations. FST values over 0.80 indicated a certain degree of genetic differentiation and reduced gene flow between populations included in the two clusters. Mean p-genetic uncorrelated distances between populations belonging to P. gracile ranged between 0.1 and 0.9%.

Table 2.
Download as
CSV
XLSX

P. (Parapholeuon) gracile – computing conventional FST from haplotype frequencies (values in bold indicate significance at the 0.05 level after Bonferroni correction).

CUB VIZ TOC INT CPO DOB VAL
VIZ 0.52155
TOC -0.07556 0.45111
INT -0.10000 0.52155 -0.07556
CPO 0.82896 0.43893 0.75370 0.82896
DOB 0.83636 0.45372 0.76364 -0.10000 0.00000
VAL 0.63047 0.25069 0.56005 0.63047 0.15691 0.63047
GRU -0.10000 0.52155 -0.07556 -0.07556 0.82896 0.83636 0.63047

Analysis of molecular variance (AMOVA), suggested some degree of genetic structure within each population (FST = 0.856, P = 0). Genetic variation among different geographic groups and among populations within each clade was 49.4% and 36.2%, with FCT = 0.493 (P > 0.05) and FSC = 0.715 (P = 0), respectively. Mantel test did not sug­gest a clear isolation by distance across the sampled region (R2 = -0.002, P > 0.05).

In P. leptodirum only four haplotypes have been identified (numbered H8–H11; Table 1), separated from each other by a maximum of five mutations, whose geographic distribution is shown in Fig. 2. The most frequent haplotype was H8, spread in all the analyzed cave populations, except for Secătura Cave. In this case, all sampled individuals were fixed for haplotype H10. In particular, samples from Coliboaia and Măgura caves were fixed for haplotype H8 that was also identified in Fânaţe (80%) and in Corbasca (55%) caves. Moreover, Fânaţe and Corbasca caves showed also haplotype H11 and H9 with a frequency of 20% and 45%, respectively. The haplotype diversity for P. leptodirum is Hd = 0.517 and nucleotide diversity was πn = 0.001. The haplotype network for this species is illustrated in Fig. 3b.

For the genetic structure of P. leptodirum, the indicators of genetic population structure (FST) are presented in Table 3; of 10 inter-populational comparisons, 7 provided significance at the 0.05 level of probability. The highest FST values (>0.70) were for comparisons between the population from Secătura Cave with the other four populations, indicating a degree of genetic differentiation between populations and a limited level of gene flow.

Table 3.
Download as
CSV
XLSX

P. (s. str.) leptodirum – computing conventional FST from haplotype frequencies (values in bold indicate significance at the 0.05 level after Bonferroni correction).

COL MAG COR FAN SEC
MAG 0.00000
COR 0.40000 0.41463
FAN 0.12438 0.13669 0.19767
SEC 1.00000 1.00000 0.73742 0.83762

Mean genetic distance between populations belonging to P. leptodirum ranged from 0.1 to 0.7%.

Analysis of molecular variance (AMOVA) suggested, also in this case, some degree of genetic structure within the populations (FST = 0.623, P = 0). Genetic variation among different geographic groups and among populations within each group was 11.9% and 50.4%, with FCT = 0.119 (P > 0.05) and FSC = 0.572 (P = 0), respectively. Mantel test did not sug­gest a clear isolation by distance across the sampled region (R2 = -0.002, P > 0.05).

Interspecific variability

As expected, genetic variation between P. gracile and P. leptodirum was much greater than intraspecific variation, with a mean genetic distance of 6.55%. Analysis of molecular variance (AMOVA) carried out considering the taxonomic assignment of each population and not the valley, suggested that the two species are well differentiated showing a genetic variation of 95.6%, with FCT = 0.955 (P = 0.002).

The phylogenetic analysis, carried out considering TrN + G model, as suggested by J model test, strongly supports the monophyly of the two species and their clear genetic separation. P. leptodirum showed a certain degree of homogenization among its populations, although each was described as a different subspecies. On the other hand, P. gracile showed a higher genetic divergence between the analyzed populations, forming two distinct clades. However, in this case the subspecies were not genetically supported, since Cubleș (P. g. gracile) and Vizu (P. g. bokorianum) did not form different clades, but are linked to the other populations, representing P. g. chappuisi (Fig. 4).

Figure 4. 

Bayesian tree constructed from 145 individuals of P. gracile and P. leptodirum from the Apuseni Mountains, belonging to 13 populations (caves). The genetic separation of the two species is clear. As outgroup an Ovobathysciola sp. from Sardinia was used.

Discussion

The genetic variability at COI DNA barcode detected in populations of both analyzed species of Pholeuon is quite low. We found a maximum of eight mutations between haplotypes of P. gracile and a maximum of five mutations separating the haplotypes of P. leptodirum. This result agrees with other studies concerning the analysis of COI DNA barcode in cave dwelling species. Intraspecific genetic variation in seven species of Bathysciola (Coleoptera, Leiodidae, Leptodirini) from Central-Southern Italian Apennines and Pre-Apennines ranged from 0 to 1.5% (Latella et al. 2017), while in Dolichopoda cave crickets (Orthoptera, Rhaphidophoridae) intraspecific genetic variation ranged from 0 to 1% (Allegrucci et al. 2005, 2014, 2021). The investigation of two species of the Tetracion troglobitic millipede (Diplopoda, Callipodida, Abacionidae) revealed a maximum of 1.4% intraspecific genetic divergence (Loria et al. 2011). Also, populations’ COI genetic divergence levels in the troglobitic Darlingtonea kentuckensis (Coleoptera, Carabidae, Trechinae) were found at 1.3% (Boyd et al. 2020).

However, despite the low variability, the studied cave populations showed a significant level of genetic structure. AMOVA analysis evidenced significant partitioning of variation within and among populations in both studied species. This result is not surprising for troglobionts because it reflects the possible barriers between the different caves and/or groups of caves to which populations are confined. On the other hand, genetic variation is not significantly partitioned among geographic groups in both species (FCT shows P > 0.05 in both cases) and Mantel test did not show a phylogeographic pattern. These results could be explained by the evolution of caves both in Pădurea Craiului Mountains inhabited by P. gracile and in Bihorului Mountains where P. leptodirum is found.

The formation of the caves is strongly related to hydrological network development and the tectonics, both at regional and local scales. The area in Pădurea Craiului Mountains where the caves are located is a highly tectonic region (Orașeanu 2020), with several stages of cave systems evolution. Caves can evolve quite differently from one another regarding water input, rock type, geotectonic features, and local hydrological system although being part of the same basin (Rusu 1981).

The hydrographic network of Pădurea Craiului Mountains is not well organized due to very intense processes of karstic caption (Orășeanu 2020) that promotes the formation of geological and hydrographic barriers, at the same time preventing gene flow between cave populations. For example, Vizu and Cubleş caves, on one hand, and Doboş/Ciur Ponor and Valău caves, on the other hand, are located in the same valley, but on the opposite slopes. Vizu Cave population showed a haplotypic composition (H3, H5, H6) completely different from Cubleș Cave (H1 and H2). The two caves are located on different slopes of the same valley, less than 1 km apart from each other. Cubleş and Vizu caves could belong to different stages of evolution and development of the Vida River (Orășeanu 2020). Both of them were carved by tributary streams (Cubleş by Blajul and Vizu by Viduţa), controlled by the local water table (e.g., different incision rates). Because the relative altitude varies for the two caves, with Cubleş at ~50 m and Vizu at ~2 m above the present waterflow the river Viduţa might act as a hydrographical barrier even during the incision of the valley and, therefore, could promote genetic isolation and differentiation of populations. We cannot assume that the paleogeographic changes were the only factors, but we can hypothesise that these could be the first step towards the isolation of populations. Ciur Ponor and Doboș caves are located on the same side of Albioara valley, specimens of the two populations share the same haplotype, H3. In Vălău cave population, located on the opposite slope, haplotype H3, but also the unique haplotype H4 have been identified. Since, in both cases, caves differ in aplotype composition, the rivers could act as hydrographic barriers, preventing the dispersal of populations. On the contrary, Tocoș, Întorsuri and Gruieț populations share haplotypes, as no hydrographic barrier separates these caves (Table 1, Fig. 2).

Bihorului Mountains, with caves hosting P. leptodirum, are mostly comprised by limestones, dolomites, conglomerates, and eruptive rock (Seghedi 2004). Still, the development of the karst network in Bihorului Mountains is characterized by intense fragmentation of the carbonate rocks with the development of large-scale karst systems and a petrographic mosaic shaping the relief. Populations of Măgura/Coliboaia caves, located on the same side of the valley, share the same haplotype (H8). At the same time, specimens from the Corbasca cave, located on the opposite slope, have the haplotype H8 and the unique haplotype H9. The populations from Fânațe and Secătura caves are situated on the same side of Bulzului valley. The caves formed in the same type of limestone are geographically close (~1 km). The first cave is inhabited by P. l. leptodirum, while the second one by P. l. problematicus. Their haplotypic composition is completely different, with Fânațe Cave population showing the most common haplotype H8 and the unique haplotype H11 in low frequency, while Secătura Cave population is fixed for the unique haplotype H10, suggesting a complete lack of gene flow. Although the two caves have certain similarities and are located at different altitudes (Secătura at 740 m.a.s.l. and Fânațe at 580 m.a.s.l.), they could have been populated at different stages by the hypothetical ancestral populations that were already genetically differentiated.

All these local features could establish geologic and hydrogeologic barriers, even for caves that are geographically close. So, even the smallest change at a certain time, in the local evolution of a cave could, be a limiting factor for the dispersal of cave populations (Sánchez-Fernández et al. 2018). Different slopes of the same river valley could represent a geographic barrier strong enough to hamper gene flow in species with poor dispersal capabilities. This could potentially be even reinforced when unsuitable habitats (i.e. impermeable strata, dry conditions, small voids etc) must be crossed. This assumption needs to be validated with additional analyses conducted on a more extensive sampling design, that would include multiple basins with multiple caves on the opposing valley slopes.

In conclusion, in both of the analyzed species, P. gracile and P. leptodirum, the genetic divergence of the COI DNA is too low to discriminate between different subspecies, although in some cases a certain degree of intraspecific genetic structure has been found (for example, between the proposed subspecies P. g. bokorianum and P. l. problematicus), suggesting that a reassessment of their status is needed. This result was expected because the genetic divergence at DNA barcoding is not informative about the species status when recently diverged species are compared and complete lineage sorting has not yet been achieved (DeQueiroz 2005; Lencioni et al. 2021). Episodes of gene flow between the different proposed subspecies represent a possible explanation preventing complete lineage sorting and most of the populations belonging to the two considered species here are not completely isolated. Moreover, the identification of species based on one single genetic marker can be incongruent with species identification using morphological characters (Moritz and Cicero 2004; Matz and Nielson 2005; Allegrucci et al. 2014; Lencioni et al. 2021).

In the Bayesian analysis of the relationships between the analyzed taxa, the two Pholeuon species are monophyletic and well differentiated from each other, with P. gracile showing a higher differentiation than P. leptodirum. As far as interspecific variation is concerned, the mean genetic distance was between 6.4 and 7.7%, with a mean value of 6.55%. Following Hebert et al. (2003) COI divergence ranges from below 1% to 16–32% be­tween species of beetles within the same genus with an average sequence divergence of 11.2. Our value falls within the intermediate range, suggesting that the two species are well differentiated and possibly for a long time. It is not rare to find high genetic differentiation between species of troglobionts as a result of allopatric speciation and the formation of hydro-geographic barriers due to changes in the surface landscape along the different climatic periods. Leptodirini appears to be an ancient cave group, as demonstrated from previous papers (Caccone and Sbordoni 2001; Ribera et al. 2010; Latella et al. 2017) with different species separated for long periods of time and accumulating genetic divergence. High genetic differentiation has been found also in a group of nine species of Bathysciola, where the interspecific genetic distances ranged from 3.1% of up to 15.1% (Latella et al. 2017). In this case, the two main lineages of Bathysciola were considered and divergence times suggested a Miocene deep cladogenesis. The two clades were shaped by the geological events during the Pliocene and the climatic changes of the Pleistocene (Latella et al. 2017). On the other hand, low levels of mitochondrial diversity have been found in species of the Pholeuon genus belonging to Drimeotus phyletic lineage although they showed deep cladogenesis with species belonging to the Drimeotus genus, included in the same phyletic lineage and revealing a possible split in the late Miocene (Bucur et al. 2003).

In conclusion, based on the observed genetic structure between the different populations of the two Pholeuon species further studies including more populations and species are needed to understand the genetic variation patterns of the group and provide valuable information for the life histories and conservation of Leptodirini in Romania. The data presented in this contribution – albeit preliminary in terms of sampling and limited in terms of genetic markers– confirm the importance of subterranean environment as a reservoir of biodiversity at a microgeographical scale. Such biodiversity should be hence managed accordingly.

Acknowledgements

The authors are grateful to Kenesz Marius and Sitar Cristian for the help with the sampling. This work was supported by a grant of Ministry of Research and Innovation, CNCS–UEFISCDI, project number PN-III-P4-ID-PCCF-2016-0016 (DARKFOOD), within PNCDI III. The authors are thankful to the reviewers for their comments and suggestions towards improving our manuscript.

References

  • Allegrucci G, Todisco V, Sbordoni V (2005) Molecular phylogeography of Dolichopoda cave crickets (Orthoptera, Rhaphidophoridae): a scenario suggested by mitochondrial DNA. Molecular phylogenetics and Evolution 37: 153–164. https://doi.org/10.1016/j.ympev.2005.04.022
  • Allegrucci G, Rampini M, Di Russo C, Lana E, Cocchi S, Sbordoni V (2014) Phylogeography and systematics of the westernmost Italian Dolichopoda species (Orthoptera, Rhaphidophoridae). Zookeys 437: 1–23. https://doi.org/10.3897/zookeys.437.7917
  • Allegrucci G, Rampini M, Chimenti C, Alexiou S, Di Russo C (2021) Dolichopoda cave crickets from Peloponnese (Orthoptera, Rhaphidophoridae): molecular and morphological investigations reveal four new species for Greece, The European Zoological Journal 88(1): 505–524. https://doi.org/10.1080/24750263.2021.1902005
  • Boyd OF, Philips TK, Johnson JR, Nixon JJ (2020) Geographically structured genetic diversity in the cave beetle Darlingtonea kentuckensis Valentine, 1952 (Coleoptera, Carabidae, Trechini, Trechina). Subterranean Biology 34: 1–23. https://doi.org/10.3897/subtbiol.34.46348
  • Bucur R, Kosuch J, Seitz A (2003) Molecular phylogenetic relationships of Romanian cave Leptodirinae (Coleoptera: Cholevidae). Atti del museo civico di storia naturale di Trieste 50: 231–265.
  • Buzilă R, Marec F (2000) Conservative chromosome number in the Cholevidae family of cave beetles. Mémoires de Biospéologie XXVII: 21–23.
  • Cocean P (2000) Munții Apuseni. Procese și forme carstice. Editura Academiei, București, Romania, 253 pp.
  • Darriba D, Taboada GL, Doallo R, Posada D (2012) jModelTest 2: more models, new heuristics and parallel computing. Nature Methods 9: e772. https://doi.org/10.1038/nmeth.2109
  • Decu V, Negrea Şt (1969) Aperçu zoogéographique sur la faune cavernicole terrestre de Roumanie. Acta Zoologica Cracoviensia XIV(20): 471–545.
  • Excoffier L, Smouse PE, Quattro JM (1992) Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics 131(2): 479–491. https://doi.org/10.1093/genetics/131.2.479
  • Excoffier L, Laval G, Schneider S (2005) Arlequin ver. 3.0: An integrated sotfware package for population genetics data analysis. Evolutionary Bioinformatics Online 1: 5–11. https://doi.org/10.1177/117693430500100003
  • Hebert PD, Ratnasingham S, de Waard JR (2003) Barcoding animal life: cytochrome c oxidase subunit 1 divergences among closely related species. Proceedings of the Royal Society of London – Series B: Biological Sciences 270(Suppl 1): S96–S99. https://doi.org/10.1098/rsbl.2003.0025
  • Howarth FG, Moldovan OT (2018) The ecological classification of cave animals and their adaptations. In: Moldovan OT, Kováč Ľ, Halse S (Eds) Cave Ecology. Springer, Amsterdam, 41–69. https://doi.org/10.1007/978-3-319-98852-8_4
  • Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R (1994) DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Molecular Marine Biology and Biotechnologies 3(5): 294–299.
  • Jeannel R (1931) Origine et évolution de la faune cavernicole du Bihar et des Carpathes du Banat. Archivio Zoologico Italiano, Atti XI Congr. Internaz. Zool. Padova 16: 47–60.
  • Kumar S, Stecher G, Tamura K (2016) MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. Molecular Biology and Evolution 33(7): 1870–1874. https://doi.org/10.1093/molbev/msw054
  • Lamoreux J, Morrison J, Ricketts T, Olson DM, Dinerstein E, McKnight MW, Shugart HH (2006) Global tests of biodiversity concordance and the importance of endemism. Nature 440: 212–214. https://doi.org/10.1038/nature04291
  • Latella L, Sbordoni V, Allegrucci G (2017) Three new species of Bathysciola Jeannel, 1910 (Leiodidae, Cholevinae, Leptodirini) from caves in Central Italy, comparing morphological taxonomy with molecular phylogeny. Insect Systematics & Evolution 49: 409–442. https://doi.org/10.1163/1876312X-00002175
  • Lencioni V, Rodriguez-Prieto A, Allegrucci G (2021) Congruence between molecular and morphological systematics of Alpine non-biting midges (Chironomidae, Diamesinae). Zoologica Scripta 50(4): 455–472. https://doi.org/10.1111/zsc.12480
  • Loria SF, Zigler KS, Lewis JJ (2011) Molecular phylogeography of the troglobiotic millipede Tetracion Hoffman, 1956 (Diplopoda, Callipodida, Abacionidae). International Journal of Myriapodology 5: 35–48. https://doi.org/10.3897/ijm.5.1891
  • Manenti R, Piazza B, Zhao Y, Schioppa E, Lunghi E (2021) Conservation Studies on Groundwaters’ Pollution: Challenges and Perspectives for Stygofauna Communities. Sustainability 13(13): e7030. https://doi.org/10.3390/su13137030
  • Mantel N (1967) The detection of disease clustering and a generalized regression approach. Cancer Research 27(2): 209–220.
  • Matz MV, Nielson R (2005) A likelihood ratio test for species membership based on DNA sequence data. Philosophical Transactions of the Royal Society B – Biological Sciences 360(1462): 1969–1974. https://doi.org/10.1098/rstb.2005.1728
  • Mammola S (2019) Finding answers in the dark: caves as models in ecology fifty years after Poulson and White. Ecography 42: 1331–1351. https://doi.org/10.1111/ecog.03905
  • Moldovan OT (2000) Révision de Drimeotus s.s. Miller, 1856 (Coleoptera, Cholevidae, Leptodirinae) de Transylvanie (Roumanie) avec description de deux nouvelles espèces et clé de détermination des taxa. Zoosystema 22(1): 139–152.
  • Moldovan OT (2008) Why so many Leptodirini (Coleoptera, Leiodidae) in Romania? In: Makarov SE, Dimitrijević RN (Eds) Advances in Arachnology and Developmental Biology – Papers dedicated to Prof. Dr. Božidar Čurčić (Eds) 12: 473–483.
  • Moldovan OT, Rajka G (2007) Historical biogeography of subterranean beetles – “Plato’s cave” or scientific evidence? Acta Carsologica 36(1): 77–86. https://doi.org/10.3986/ac.v36i1.210
  • Moldovan OT, Racoviță G, Dunay G (2007) Reconsidering Pholeuon C. Hampe (Coleoptera: Leiodidae: Cholevinae), with the description of a new subgenus. Zootaxa 1449: 31–43. https://doi.org/10.11646/zootaxa.1449.1.2
  • Moldovan OT, Iepure S, Brad T, Kenesz M, Mirea IC, Năstase-Bucur R (2020) Database of Romanian cave invertebrates with a Red List of cave species and a list of hotspot/coldspot caves. Biodiversity Data Journal 8: e53571. https://doi.org/10.3897/BDJ.8.e53571
  • Myers N, Mittermeier RA, Mittermeier CG, Da Fonseca GAB, Kent J (2000) Biodiversity hotspots for conservation priorities. Nature 403: 853–858. https://doi.org/10.1038/35002501
  • Orășeanu I (2020) Hidrogeologia carstului din Munții Apuseni. Editura Belvedere, Oradea, România, 350 pp.
  • Ribera I, Fresneda J, Bucur R, Izquierdo A, Vogler AP, Salgado JM, Cieslak A (2010) Ancient origin of a western Mediterranean radiation of subterranean beetles. BMC Evolutionary Biology 10: e29. https://doi.org/10.1186/1471-2148-10-29
  • Racoviță Gh (1996) Révision systématique des leptodirinae souterraines des Monts Apuseni. II. Le sous-genre Parapholeuon Ganglb. Du bassin de Crișul Repede (Monts Pădurea Craiului). Trav. Inst. Spéol. “Emile Racovitza” t. XXXV: 69–105.
  • Racoviță Gh (1998–1999) Révision systématique des Bathysciinae souterrains des Monts Apuseni. III. Le sous-genre Parapholeuon Ganglb. Du basin de Crișul Negru (Monts Pădurea Craiului). Travaux de l'Institut de Spéologie “Emile Racovitza” t. XXXVII–XXXVIIIȘ: 175–216.
  • Racoviță Gh (2010) Révision systématique des Leptodirinae souterrains des Monts Apuseni. VII. Le sous-genre Pholeuon (s. str.) du bassin de Crișul Negru (Monts du Bihor). Travaux de l'Institut de Spéologie “Emile Racovitza” t. XLIX: 3–27.
  • Racoviță Gh (2011) Révision systématique des Leptodirinae souterrains des Monts Apuseni. VIII. Apercu synthétique sur le genre Pholeuon. Travaux de l'Institut de Spéologie “Emile Racovitza” t. L: 37–59.
  • Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA (2018) Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Systematic Biology 67(5): 901–904. https://doi.org/10.1093/sysbio/syy032
  • Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP (2012) MrBayes 3.2: Efficient Bayesian Phylogenetic Inference and Model Choice Across a Large Model Space. Systematic Biology 61(3): 539–542. https://doi.org/10.1093/sysbio/sys029
  • Rusu T (1981) Les drainages souterrains de Monts Pădurea Craiului. Systematic Biology “E. Racovitza” XX: 187–205.
  • Sánchez-Fernández D, Rizzo V, Bourdeau C, Cieslak A, Comas J, Faille A, Fresneda J, Lleopart E, Millán A, Montes A, Pallares S, Ribera I (2018) The deep subterranean environment as a model system in ecological, biogeographical and evolutionary research. Subterranean Biology 25: 1–7. https://doi.org/10.3897/subtbiol.25.23530
  • Seghedi I (2004) Geological evolution of the Apuseni Mountains with the emphasis on the neogene magmatism. In: Cook NJ, Ciobanu CL (Eds) Au-Ag-telluride Deposits of the Golden Quadrilateral, Apuseni Mts., Romania.
  • Tamura K, Stecher G, Peterson D, Filipski A, Kumar S (2013) MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Molecular Biology and Evolution 30(12): 2725–2729. https://doi.org/10.1093/molbev/mst197
  • Weir BS (1996) Genetic data analysis II: Methods for discrete population genetic data. Sinauer Associates, Inc., Sunderland.
login to comment