Phylogeographic and genetic insights into Sinonychia martensi: an endemic cave-dwelling harvestman in Beijing

Abstract

Caves are one of the most exciting environments on earth, often considered an evolutionary laboratory due to the suite of convergent adaptive traits (troglomorphisms) of organisms inhabiting them. Sinonychia martensi Zhang & Derkarabetian, 2021, is the first and only Travunioidea species recorded in China and is endemic to Beijing, being known from multiple caves. However, nothing is known regarding its phylogeographic or evolutionary history. In this study, we assessed the species boundaries of S. martensi from nine caves using morphological and molecular methods to elucidate its phylogenetic position and genealogical relationships. We also investigated the genetic diversity, population genetic structure and demographic history of S. martensi to clarify the population-level relationships and make inferences about historical phylogeography. The results indicate that the species from different caves all belonged to S. martensi but represent different populations. These populations exhibit strong population structure and low genetic diversity. Cave populations may share a common ancestor and multiple independent invasions to different caves. The diversification within S. martensi was likely driven by climate change and subtropical evergreen broadleaf forests associated with the middle Miocene. This study highlights the need for further conservation efforts and exploration in Beijing caves.

Introduction

Caves are some of the most adverse environments on Earth. The restricted access to food, the extreme conditions of constant darkness and humidity make these habitats very challenging for living organisms [40]. Cave-dwelling species have historically attracted the attention of evolutionary biologists due to their bizarre ‘regressive’ characters [81], convergent evolution, small population size and low reproductive capacity relative to species found outside caves. Additionally, their species numbers tend to face the risk of extinction due to climate change and human factors [48]. Despite these challenges, many taxa, including arthropods and vertebrates, have colonized these subterranean environments [15]. Within the arachnid order Opiliones, the Travunioidea includes many cave-obligate taxa, with species from 14 genera showing some degree of troglomorphism, possessing homoplastic morphological features that evolve in response to cave life [18]. Convergent evolution in troglomorphic morphology has been recorded among different lineages in the order Opiliones, such as the paranonychid harvestman from the western United States, taxa in the Laniatores family Phalangodidae and the European Ischyropsalis C.L. Koch, 1839, among others [19, 34, 73]. A strong positive correlation between the degree of troglomorphism and the time spent since cave invasion was found in the genus Sclerobunus Banks, 1898, raising the interesting possibility of predicting taxon age from a simple troglomorphism index in this group [19]. However, cave ages cannot be relied upon to place upper limits on cave species ages, as studies have shown that troglomorphic species are sometimes older than the caves they inhabit [38, 80].

The convergent evolution in troglomorphic morphology significantly impacts trait identification, especially when the degree of intraspecific differentiation is great and the morphological differences of closely related species are small. There are many shortcomings in determining species by traits alone. With the development of molecular data, population genetics and ecology, various species delimitation methods have been derived. Previously, some cave-dwelling harvestmen have been regrouped and redefined using a variety of species-delimitation methods combined with morphological data [33]. Most species delimitation methods based on a single locus like COI are prone to errors [33]. Based on COI, the morphological identification of Minisge sagai Cruz-López, Monjaraz-Ruedas & Francke, 2019 living at shallower depths (20–200 m) was consistent with the species delimitation results, whereas the species delimitation results were not consistent with the morphological results for Minisge kanoni Cruz-López, Monjaraz-Ruedas & Francke, 2019 distributed at middle depths (400–600 m) without much gene flow between caves [14]. Species delimitation may overestimate the number of species in cases of unusually high genetic variation [26].

Geographic barriers [36, 83], ecological differences [42], historical processes [8, 35], and the dispersal potential of species [28, 61], are the main evolutionary forces affecting the population genetic structure and phylogeographic distribution of extant species. Understanding these factors in studies of cave-dwelling species can enhance our understanding of their distribution dynamics, which is critical for developing effective conservation strategies. To date, phylogeographic analyses have been conducted on many Opiliones [5, 16, 17, 58, 59, 79].

Sinonychia martensi [85], recorded in Beijing’s Tangren Cave, is the first species of the family Cladonychiidae Hadži, 1935 and superfamily Travunioidea Absolon & Kratochvíl, 1932 described from China by Zhang & Derkarabetian [85]. Subsequently, the experimental team found it in other caves near Tangren cave (Fangshan and Mentougou Districts), but it has not yet been found outside of caves or in the southern area. Based on morphological study, Zhang & Derkarabetian [85] proposed that this species was closely related to Speleonychia [6], a highly troglomorphic species found in lava tubes in North America and not to any other taxa from South Korea or Japan. It was therefore established as a new genus: Sinonychia [85], although it has not been included in any comprehensive molecular phylogenetic analyses.

In this study, we investigated the phylogenetic position of S. martensi in Travunioidea using phylogenomic data, inferred the relationships among populations and estimated the genetic diversity and population structure in the group. Additionally, we used a phylogeographic framework to infer the possible historical events that promoted the divergences of the major lineages and shaped the population genetic diversity in the genus. Our study revealed populations exhibit strong population structure and low genetic diversity. We hypothesize that before this, the common ancestor of S. martensi was forced to invade separate caves and survive to this day due to drastic climate changes.

Materials and methods

Study area and specimen sampling

We sampled 139 individuals from nine caves, most of which were located in Beijing’s Fangshan District, with the exception of one cave, Qiubo Cave, situated in the Mentougou District (Fig. 1). The rock types in these two districts are mainly sedimentary and metamorphic rock, with sedimentary rock in particular occupying a large proportion. Specimens were manually collected during multiple trips to the caves, where they were found under stones or on walls (Fig. 2a, b, c). For data analysis, S. martensi specimens from different caves were designated using English letter combinations as follows: Bianfu Cave (BFD), Siyu Cave (SYD), Xi Cave (XD), Xionggu Cave (XGD), Qiubo Cave (QBD), Tangren Cave (TRD), Jinhua Cave (JHD), Beixin Cave (BXD), and Nishui Cave (NSD). Samples used in the experiment were collected with the assistance of laboratory members from 2018 to 2022. Detailed information on the number of samples collected from each cave, along with their geographical locations and environmental characteristics, can be found in Supplement 1. All specimens were preserved in 95% ethanol and stored at -20 ℃ upon their return to the laboratory.

Fig. 1

Geographic distribution of nine cave populations and view of Xi Cave entrance. The letters indicated the initials of the cave name. QBD: Qiubo Cave; TRD: Tangren Cave; JHD: Jinhua Cave; XGD: Xionggu Cave; XD: Xi Cave; SYD: Siyu Cave; BFD: Bianfu Cave; BXD: Beixin Cave; NSD: Nishui Cave

Fig. 2

Living specimens (a-b) and habitat of S. martensi, in Xionggu cave (c) and Qiubo cave (d) (©A. Ruoyi Xiao, B. Shanfeng Zhang)

Taxonomy

External morphological identification of all specimens was performed using a Leica M205A stereomicroscope [85]. Male genitalia were extracted by first removing the genital operculum and carefully excising the surrounding tissue to expose the glans. The sharp dissection needle was then used to cut open the sheath of the penis, allowing for easier examination. A total of 11 male specimens were examined under the stereomicroscope. To date, no males have been found in Jinhua Cave, Xionggu Cave, or Qiubo Cave. Two male genitalia from each of the remaining six caves were dissected and photographed (including one specimen from Tangren Cave). Photographs of male genitalia were captured using an Olympus microscope equipped with a KUY NICE CCD camera. Individual images were then compiled into a composite image using Helicon Focus (http://www.heliconsoft.com/helicon/heliconfocus.html) and further edited with Adobe Photoshop CS3.

DNA extraction, amplification and sequencing

We created two molecular data sets: one based on Ultraconserved Elements (UCEs) and the other using traditional multi-locus Sanger sequencing. For all newly sequenced samples, genomic DNA was extracted using the QIAGEN DNeasy Blood & Tissue Kit, and DNA quantity was assessed with a Qubit™ Fluorometer.

For two specimens of S. martensi we performed low coverage genome sequencing with short reads. For these samples, genomic DNA was sent to Novogene Co. Ltd. for library preparation using the Truseq Nano DNA HT Sample Preparation Kit (Illumina USA), followed by sequencing on the Illumina NovaSeq platform with 150 bp paired-end reads and an insert size of approximately 350 bp.

Genome assembly was performed using the PLWS pipeline [88]. Quality trimming of reads was carried out using bbduk.sh (BBTools, [7]), with reads shorter than 15 bp, containing more than 5 Ns, or with poly-A or poly-T tails of at least 10 bp being trimmed. Genome contigs were assembled using multiple k-mer strategies (k-mer as) in Minia v3.2.1 [10]. Redundans v0.13c [63] was employed to identify and remove contigs representing high heterozygosity. Contig scaffolding and gap filling were performed with BESST v2.2.8 [70] and GapCloser v1.12 in the SOAPdenovo2 suite [50] respectively. Sequences shorter than 500 bp were excluded from subsequent analyses.

Following the PHYLUCE workflow [24], UCEs were bioinformatically extracted from the two assembled genomes using the Arachnida probe set [25, 76]. During the alignment of the Arachnida probes to the two Sinonychia genomes, coverage and identity were both set to 75, with 500 bp on either side extracted. The extracted UCE loci were then aligned back to the probe set with both min coverage and min-identity set to 65 to remove duplicated UCEs and determine the final orthologs for phylogenetic analyses.

Of all the samples collected, 72 were used for multi-locus Sanger sequencing experiments. To minimize operational errors, ZHC043 and ZHC045 were subjected to repeated experiments. We selected one mitochondrial DNA (mtDNA): cytochrome c oxidase subunit I (COI) and three nuclear DNA (nuDNA): histone H3 (H3), internal transcribed spacer subunit II (ITS2), and 28S rDNA for the molecular phylogenetics of samples from nine caves. The COI gene was used for subsequent population structure analysis. DNA was extracted from the right leg using the Tiangen kit following standard protocol. Polymerase chain reaction (PCR) was performed in 25 μl reactions consisting of 6.9 μl ddH₂O, 12.5 μl mix, 0.8 μl of each primer and 4 μl DNA template.

PCR products were verified by 1% agarose gel electrophoresis. For COI amplification, we used the primers LCO1490 and HCO2198 with the following PCR conditions: 94 ℃ for 5 min, 35 cycles at 94 ℃ for 30 s, 45 ℃ for 40 s and 72 ℃ for 1 min, followed by a final step of extension step at 72 ℃ for 7 min and a hold at 4 ℃. The H3 PCRs were performed with the primers H3aF and H3aR under the conditions: 94 ℃ for 2 min, 35 cycles at 94 ℃ for 30 s, 52 ℃ for 30 s and 72 ℃ for 1 min, with a final extension at 72 ℃ for 3 min and a hold at 4 ℃. For ITS2, the primers were 5.8S2 and 28S2, with PCR conditions of 95 ℃ for 2 min, 40 cycles at 95 ℃ for 30 s, 48 ℃ for 30 s and 72 ℃ for 1 min, followed by a final extension at 72 ℃ for 5 min and a hold at 4 ℃. The 28S primers were ZX1 and ZR2, with PCR conditions of 94 ℃ for 3 min, 40 cycles at 94 ℃ for 30 s, 62 ℃ for 45 s and 72 ℃ for 1 min 15 s, followed by a final extension at 72 ℃ for 5 min and a hold at 4 ℃ (Table 1).

Table 1 The primers of the target genes

Sequences chromatograms of the target fragments were edited and visualized using the Mesquite 3.61 package [51] Chromaseq [52], Phred [21, 32]) and Phrap [31]. We performed manual quality checks on the chromatograms of each sequence and organized the target fragments into separate folders. All sequences were submitted in GenBank, with accession numbers provided in Supplement 2. We used MAFFT v7.450 [41] to align sequences with the “G-INS-I” strategy and checked for stop codons in COI and H3 by translating them into amino acid sequences using Geneious Prime 2022.1.1 [43]. Before phylogenetic analysis, all four amplified target genes were concatenated with PhyloSuite v1.2.2 [87].

Phylogenetic analyses

We included the two short-read assemblies in a dataset with all previously sequenced UCE samples of Travunioidea (and outgroups) from the UCE analyses of Derkarabetian et al. [18]. The assembled and aligned UCE contigs were processed using default settings in phyluce 1.7.1 [24] with conservative gblocks settings (–b1 0.5 –b2 0.85 –b3 4 –b4 8). Only loci with at least 50% of samples represented in a locus were included in the final analyses after manual inspection in Geneious Prime 2022.1.1 (www.geneious.com). A concatenated matrix was used for maximum likelihood analyses in IQ-TREE v2 [56], employing 1000 ultrafast bootstrap replicates [37] and merging partitions with identical models identified via ModelFinder [46].

For the Sanger data, all aligned COI sequences were analyzed for genetic distance in MEGA X [45]. The best substitution model for our dataset was selected using the Bayesian information criterion (BIC) in Partitionfinder within Phylosuite v1.2.2. Phylogenetic analyses for the concatenated matrix, including COI, H3, ITS2 and 28S rDNA sequences, were conducted using both Bayesian inference (BI) and maximum likelihood (ML). COI was partitioned by codon positions. The ML analysis was performed using IQ-TREE v1.6.8 [56] with 5000 ultrafast bootstraps within Phylosuite v1.2.2. BI trees were constructed with MrBayes v3.2.6 [69] within Phylosuite v1.2.2 with two independent runs and four Markov Chain Monte Carlo (MCMC) chains. Speleonychia sengeri [6] (Cladonychiidae Hadži, 1935) and Briggsus flavescens Briggs, 1971 (Cladonychiidae Hadži, 1935) were used as outgroups and their sequences were downloaded from GenBank. The analysis ran for 10⁶ generations, sampling every 1000 generations, with the initial 25% discarded as burn-in. Convergence was assessed in Tracer v1.7.2 [66], considering an effective sample size (ESS > 200) as satisfactory. The obtained dendrograms were visualized and edited in FIGTREE v1.4.3 [65].

Molecular species delimitation

We performed the molecular species delimitation for S. martensi from nine caves using three methods: Assemble Species by Automatic Partitioning (ASAP) [64],Generalized Mixed Yule Coalescent (GMYC) [60], and Bayesian Phylogenetics and Phylogeography (BPP) [82], based on the four genes obtained.

ASAP, an advancement in Automatic Barcode Gap Discovery (ABGD), is based on pairwise genetic distances and does not require any prior information about the number of species, phylogenetic trees, or predefined genetic distances. We performed ASAP analysis using three models, Jukes-Cantor (JC69), Kimura (K80) TS/TV and P-Distance, respectively.

For GMYC analysis, we used the single-threshold version of the likelihood approach. The best nucleotide substitution model for each gene fragment was first determined using Modelfinder within Phylosuite v1.2.2, and then an ultrametric tree was constructed as the tree model. The aligned sequences and ultrametric tree were obtained using BEAUti v1.10.4 in BEAST v1.10.4 [77], and three output files were created. The.xml file was run in the BEAST v1.10.4, generated three output files and the file (.log.txt) was checked the convergence in Tracer v1.7.2. The phylogenetic tree was generated with TreeAnnotator v1.10.4 in BEAST v1.10.4. The phylogenetic tree was analyzed using the R package “splits” to obtain species delimitation results, which were compared with those from other methods.

In account for the possibility of different gene trees and species trees due to incomplete lineage sorting of species, we used BPP, a multigene-based species delimitation method. Four genes were included in the analysis and they were COI, H3, ITS2 and 28S. All individuals were divided into seven hypothetical species based on geographic distribution and phylogeny. Unguided Bayesian species delimitation (A11) was performed using BPP v 4.3 to explore changes in species tree topology and species delimitation models. To improve accuracy, we used three priors in this study [86]:

• θ ~ G (1, 10), τ ~ G (1, 10);

• θ ~ G (2, 2000), τ ~ G (2, 2000);

• θ ~ G (1, 10), τ ~ G (2, 2000).

Molecular dating

Divergence times were estimated using the COI and ITS2 dataset with BEAST v1.10.4. We applied mean substitution rates previously estimated for Gonyleptidae harvestman: COI: 0.0055 substitutions per million years (Mya) (range:0.003–0.008) and ITS2: 0.0004 substitutions per million years (Mya) (range: 0.0002–0.0006) [9]. BEAST analyses were conducted with the GTR + G model for COI and the JC model for ITS2. The tree and clock models were linked (site model unlinked). A strict clock with a fixed value and a Yule model process were used as tree priors. All parameters were set and saved as a.xml file. The.xml file was run in the BEAST v1.10.4, generated three output files and the file (.log.txt) was checked for convergence in Tracer v1.7.2. The analysis ran 7 × 10 ⁷ generations, with the initial 10% discarded as burn-in. Parameters’ effective sample size (ESS) values were viewed in Tracer v1.7.2 to confirm that ESS > 200. Maximum clade credibility (MCC) species and gene trees were annotated using TreeAnnotator v1.10.4 and edited in Figtree v1.4.3.

Population structure, diversity and demography

COI sequences of 139 samples were used for population structure analysis. We constructed haplotype networks for COI using the TCS methods [11] in the software PopArt v1.7 [47]. DNAsp v6.12.03 [49] was used to calculate the number of haplotypes, nucleotide diversity (π), haplotype diversity (Hd), nucleotide differences (K), and the number of segregating sites (S) for each cave.

The demographic history of S. martensi was studied using mismatch distribution in Arlequin v3.5 [22] to determine the distribution of frequencies of pairwise differences. Neutrality indices, Tajima's D and Fu's Fs are widely used in molecular analysis to detect whether populations are under recent expansion. Tajima's D test [78] is derived from mutation (segregating sites) frequencies, while Fu's Fs test [29] is derived from haplotype distribution [67]. All parameters were evaluated based on 1000 bootstrap replicates.

Analysis of molecular variance (AMOVA; [23]) implemented in Arlequin v3.5 was employed to estimate the genetic variation among and within populations. F-statistic values (F_ST; differentiation index) between populations were calculated based on COI datasets using Arlequin v3.5. The conventional population F_ST comparisons were measured with 1,000 permutations at a significance level of 0.05.

Results

Taxonomic status of S. martensi

Sinonychia martensi was assigned to the Cladonychiidae family based on its distinctive morphological characteristics. Additionally, its intestinal midgut structure and the male genitalic morphology bear a strong resemblance to those of Briggsus and Speleonychia [85]. However, the convergent nature of cave life may have contributed to the morphological similarities between geographically distant taxa concerning the typical troglomorphic traits like eye loss, depigmentation, and appendage elongation. To confirm the phylogenetic position of S. martensi within Travunioidea, and to further explore the evolutionary history of the genus, we conducted a molecular phylogenetic study based on ultraconserved elements (UCEs). UCE analyses revealed Sinonychia as the sister lineage to Speleonychia + Briggsus + Isolachus, all from the Pacific Northwest of North America (Fig. 3), largely confirming the hypotheses of Zhang and Derkarabetian [85] based on genital morphology.

Fig. 3

Phylogenetic position S. martensi in Travunioidea. Maximum likelihood tree based on the UCE dataset and the asterisk indicates ultrafast bootstrap = 100

Male genital morphology

The genital structure of this species is relatively simplified, with only two pairs of thick, recurved setae on the glans (Fig. 4). The apical pair of setae are smaller and nearly perpendicular to the axis of the stylus, while the subapical pair are larger and more strongly inclined towards the base of the penis [85]. Additionally, we observed a puzzling phenomenon in the sex ratio of collected specimens. The proportion of males to females varied significantly among the nine caves investigated as part of this study: in Bianfu Cave, the ratio can reach 7:3, but in the Xionggu Cave, we have not found a single male among the approximately 30 specimens collected. Based on the results of the morphological analyses, we were unable to distinguish differences between harvestmen from different caves and concluded that the harvestmen from nine caves were the same species.

Fig. 4

Dorsal view of penis for S. martensi in six caves. Except for Tangren cave, other caves have two individuals respectively

Phylogenetic analyses

We obtained sequences for COI, H3, ITS2, 28S measuring 620 bp, 330 bp, 429 bp, 1123 bp, respectively. The COI analysis of variable and conserved sites showed that populations of S. martensi were relatively conserved, with 108 out of 620 total sites being parsimony-informative. Inter-cave genetic distances (K2P model) were high, ranging from 0.2% to 12% based on COI datasets (Table 2; Supplement 3).

Table 2 Estimated COI divergence from MEGAX analyses based on the Kimura 2-Parameter model

The morphological results were consistent with intrafamilial phylogeny based on UCEs. We conducted phylogenetic analyses with BI (Fig. 5) and ML (Fig. 6) based on the concatenated dataset of four gene sequences, and both showed the monophyly of S. martensi with strong support (Bayesian Posterior Probability, BPP = 1). Both analyses recovered similar topologies, except for the sample ZHC095 JHD. Within the S. martensi group, we observed a major divergence separating species into three clades in both trees: the first clade formed by BFD, the second by SYD, and the third comprising the remaining caves. Most caves corresponded to monophyletic lineages, with the exceptions of JHD vs. TRD and BXD vs. NSD, which were genetically mixed. Consultations with local professional cavers revealed a water connection between Jinhua Cave and Tangren Cave, possibly explaining the genetic relationship between samples from these caves. An interesting line of inquiry would be to experimentally test whether flood waters carry floating debris capable of transporting individuals to new localities. Similarly, Nishui Cave and Beixin Cave, which are geographically close (650 m away from the map), are closely related phylogenetically and form a large clade. Geographic proximity generally coincides with closer phylogenetic relationships within the NSD-BXD group, suggesting some degree of connection due to their geographical location. Interestingly, samples ZHC043 and ZHC045 from Tangren Cave were placed unexpectedly within the large clade including Nishui Cave and Beixin Cave. To test for possible laboratory errors, we re-extracted DNA from samples ZHC043 and ZHC045 and redid PCR (naming these samples ZHC043A and ZHC045A). The results remained the same, raising intriguing questions.

Fig. 5

Bayesian tree based on the concatenated dataset. Numbers at nodes represent posterior probabilities

Fig. 6

Maximum likelihood tree based on the concatenated dataset. Numbers at nodes represent ultrafast bootstrap values

Species delimitation

Three molecular species delimitation methods analyses (ASAP, GMYC and BPP) were used to test the species hypothesis on samples from nine caves based on single gene fragments (COI, H3, ITS2 and 28S) and combined gene fragments [44, 54]. The results showed differences between the three methods and highlighted that different gene greatly influenced the outcomes.

The ASAP results based on the three models (JC69, K80 and P-Distance) were highly similar. Therefore, we combined and selected the best of them. Based on mtDNA COI, ASAP tended to classify OTUs by cave: specimens from BFD, SYD, and XD were each divided into one OTU, specimens of TRD and JHD, connected by hydrological systems, were divided into one OTU, specimens of BXD and NSD, geographically close, were divided into one OTU, and specimens of QBD and XGD, phylogenetically sister to each other, were divided into one OTU (Fig. 7). Based on ITS2, ASAP produced a relatively conservative division with only two OTUs, aligning more closely with the morphological classification results. The 28S-based ASAP analysis divided all samples into six OTUs, primarily based on caves, except for ZHC031_SYD in Siyu Cave, which was a single OTU. GMYC analysis, based on the principle of tree topology, classified samples from nine caves into 8 OTUs (COI), 30 OTUs (H3), 8 OTUs (ITS2), and 10 OTUs (28S). Similar to ASAP, GMYC tended to classify specimens within each cave as an OTU based on COI (Fig. 7). Delimitation results were more complex for H3, ITS2, and 28S, but both ITS2 and 28S classified specimens from BFD and SYD each as a single OTU. The BPP delimitation method, based on multi-locus sequence data, divided all samples into five OTUs. Similar to ASAP, BPP also grouped the same caves as a unit. Jinhua Cave and Tangren Cave, Qiubo Cave and Xionggu Cave were classified as one OTU respectively, Nishui Cave, Beixin Cave, and Xi Cave were classified as one OTU.

Fig. 7

Species delimitation based on molecular data. In the same analytical method, the same color represents the same OTU group. The phylogenetic tree on the left is the ML tree based on the concatenated dataset. The purple triangles at nodes represent ultrafast bootstrap values

Molecular dating

Our molecular dating based on COI data indicated that S. martensi diverged from its closet Cladonychiidae relatives, S. sengeri and B. flavescens, approximately 33.38 million years ago (Mya) with a 95% highest posterior density (HPD) interval of 26.21–40.54Mya. Intraspecific divergences within S. martensi were dated to around 15.51Mya (95% HPD = 12.18–18.83Mya) during the middle Miocene (Fig. 8). Among the nine caves in this study, the earliest divergence was observed between populations in Bianfu Cave and Siyu Cave, which occurred approximately 8.94 Mya (95% HPD = 6.28–11.60Mya). Subsequently, around 8.60 Mya, another major divergence event led to the formation of a large branch including populations from TRD, JHD, QBD, NSD, BXD, XD, and XGD.

Fig. 8

Calibrated COI Bayesian tree obtained in the BEAST analysis. Numbers above the lineages and in the brackets represent estimated divergence dates, the 95% highest posterior density (HPD) of each node respectively. Posterior probabilities less than 95 are shown with red node bars. Inset map of Fangshan District and Mentougou District with black triangles representing sampled caves

Intraspecific analyses

Intraspecific analyses were conducted for the nine populations of S. martensi, revealing varying levels of genetic variation. A total of 22 haplotypes were identified for the COI gene. The haplotype network shows a clear division between different caves where S. martensi is present. Only haplotype 11 is shared by three caves, while other haplotypes are generally exclusive to a single cave (Fig. 9). The haplotype network diagram indicates that the haplotype branch clustering is consistent with the phylogenetic tree, and the 139 samples exhibit a strong population structure. The mean haplotype diversity (Hd) across all populations was 0.3631 ± 0.122. The mean nucleotide diversity (π) was 0.00287 and the average number of nucleotide differences (K) was 1.79 (Table 3).

Fig. 9

Haplotype network constructed using TCS software based on COI sequences. Twenty-two haplotypes were found in S. martensi. Each haplotype is represented by a circle, and the size of each circle is proportional to the number of individuals with that haplotype. Each line connecting two haplotypes represents a single mutation, and dots indicate hypothetical missing haplotypes

Table 3 Diversity indices and neutrality test obtained for S. martensi

The overall mean F_ST is 0.927 (p < 0.01), significantly larger than 0.25. Lower F_ST values were observed between geographically adjacent populations (Beixin Cave and Nishui Cave F_ST = 0.233; Table 4) and connected populations (Tangren Cave and Jinhua Cave F_ST = 0.075; Table 4). F_ST values close to 1 or equal to 1 between most caves indicate substantial genetic differentiation or completely separate populations (e.g., Xi Cave vs Siyu Cave, Siyu Cave vs Xionggu Cave, Xi Cave vs Qiubo Cave). Pairwise F_ST derived from the AMOVA suggested relatively high and significant levels of genetic differentiation among populations but not among individuals of the same populations (Table 5).

Table 4 F_ST values among pairs of populations

Table 5 Results of analysis of molecular variance (AMOVA) in S. martensi based on COI data

To uncover the demographic history of the nine populations, neutrality tests were conducted using Tajima’s D and Fu’s Fs statistics based on the COI datasets. The Tajima’s D values for other cave samples were not significant (p > 0.05). Only the values of Tajima's D and Fu’s Fs in Bianfu Cave were significantly negative, supporting an excessive number of rare haplotypes resulting from a recent population expansion (Fig. 9; Table 3). The significantly negative neutrality test values and unimodal mismatch distribution pattern of Bianfu Cave suggest that the population fits the sudden expansion model (Fig. 10; Table 3). The mismatch distribution of all samples may indicate that the nine populations in China fit a neutral evolution model (Fig. 10).

Fig. 10

Mismatch distribution of S. martensi for entirety nine caves (a) and Bianfu cave (b)

Discussion

Phylogenetics of Sinonychia

Caves are a unique and important type of ecosystem on Earth, providing a stable and enclosed habitat for cave organisms. This isolation contributes to species endemism and uniqueness, drives convergent evolution of morphologies, and serves as a natural laboratory for studying ecology, evolution and biomedicine [12]. Cave animals have historically evolved within these environments, maintaining stable populations. In particular, S. martensi specimens from the nine caves in Beijing are almost identical in morphological structure, particularly the male genitalia (Fig. 4).

Populations of S. martensi show variable genetic distances ranging from 0–1.5% within populations (except ZHC043 and ZHC043A) and 0.2%–12% between different populations based on COI datasets (Supplement 3; Table 2). The COI-based genetic distance range between the nine caves is slightly higher than that of other surface-dwelling Opiliones species [1, 26, 30, 72]. This may be due to the isolated and closed cave environment, which hinders gene flow between populations, coupled with the low dispersal ability of Opiliones species. Over evolutionary time, this leads to the formation of caves as "island" populations [75]. Despite significant genetic divergence, no morphological differences were found, indicating niche conservatism.

Although genetic distances between species were large, phylogenetic results based on multilocus data suggested that S. martensi individuals from the same caves were more closely related to each other than to individuals from different caves (Figs. 5 and 6). The cave, as a relatively closed environment, not only protects the organisms inside but also isolates them from outside influences. We found that populations from Siyu Cave and Bianfu Cave are more distantly related to the other seven caves. However, there may be a genetic exchange between these two caves due to underground passages. Populations from Beixin Cave and Nishui Cave are the closest genetically, and inbreeding within them is expected. However, the genetic exchange between Tangren Cave and these closer caves is hard to understand (Fig. 1), even though repeated validation confirmed the accuracy of these results. Further exploration is needed to understand these patterns, considering factors such as the number of samples collected, the presence of underground passages including rivers, and the activities of other organisms.

Due to the extreme similarity in morphology, the use of molecular information to delimit species boundaries and to verify the validity of morphospecies has become an effective means of rapidly identifying and describing species, thereby accelerating the discovery and understanding of the immensely diverse biological world [74]. However, taxa restricted to cave habitats may present special challenges to species delimitation [33]. In this study, we performed ASAP, GMYC and BPP analyses based on four genes. The results were not as expected, as the different methods for defining species were not entirely consistent. ASAP results were relatively conservative, mostly categorizing caves as a whole, especially for COI. In contrast, only the COI-based GMYC delimitation resulted in cave-related results, with the rest having more categorical units, and even 30 OTUs in the H3-based delimitation. It may be because GMYC is so sensitive to the delimitation of conserved genes and branch length that it is difficult to distinguish them. BPP used the consistent lineage formed by different genes as the basis for species delimitation, providing more phylogenetic information compared to single-gene analysis. The BPP analysis categorized all samples into 5 OTUs. Multiple factors can affect the results of species delimitation, divergence among included samples, type and number of genes, number of haplotypes, geographic distance and others. These factors interact in complex ways to influence delimitation results [84].

In this study, species delimitation based on molecular data did not align with morphological differences, which is a pattern that can be present in cave environments. Even though previous researchers delimited different harvestmen inhabiting deeper passages and classified them into separate species based on COI, they did not identify differences in morphology and ultimately defined them as the same species, Minisge kanoni. [14]. We were unable to distinguish the harvestmen from the nine caves based on morphological results; the phylogeny and species delimitations were mostly based on the caves as the taxonomic unit, especially Siyu Cave and Bianfu Cave. We concluded that the harvestmen from the nine caves were all S. martensi, and that they showed high population structure. Similar examples have been described in previous studies, where BPP has been used to delimite populations and evolutionary lineages within populations [55, 57, 62].

Cave evolution in Sinonychia

Most caves are formed by chemical reactions between circulating groundwater and surrounding rocks, resulting in isolated and strongly zonal environments characterized by dim light, oligotrophic conditions, and small daily and annual variations in temperature and humidity [48]. Due to these challenging conditions, colonization within such an environment requires highly specialized adaptations [39]. Consequently, obligate cave species tend to have small population sizes and low reproductive capacity compared to species outside caves, making them vulnerable to extinction from climate change and human activities [20]. Globally, most caves (approximately 93%) are not protected areas [71]. During our collection, the internal environment of Qiubo Cave and Xi Cave was significantly disturbed by human activities (Fig. 2d). Nucleotide and haplotype diversity in Xi Cave and Qiubo Cave was 0, probably attributable to human activities, although further evidence and analysis are needed to confirm this hypothesis (Fig. 2d,Table 3). Human disturbances can greatly affect the genetic diversity of cave species, which generally have small population sizes [2].

Caves can be seen as separate “small worlds” with distinct environmental characteristics, impacting the morphology and genetic characteristics of the organisms within. S. martensi located in different caves belong to different populations, and the origin of these populations and the links between them are worth exploring. Molecular dating shows the diversification within S. martensi occurred in the Miocene at approximately 15.51 Mya (Fig. 8), a period characterized by a warm and moist subtropical monsoon climate, along with mountain building and river incision promoting the development of large caves in subtropical East Asia [89]. This period also coincides with the development of the East Asian subtropical evergreen broadleaf forests (EBLFs) and Neogene climate changes [48]. The new tectonic movement laid the foundation for the emergence of the present-day karst caves in Beijing. The long peneplain process of the Paleogene period resulted in the emergence of a highly meandering Cenozoic River that impacted the mountainous area of western Beijing. The uplift of the West Mountain in Beijing during the Neozoic period, coupled with the strong humid climate, caused severe erosion of the near-plain, and the filling material was eroded, forming extensive cave areas. During the Quaternary, the rapid uplift of the West Mountain, combined with constant river cutting, led to the formation of multi-layered caves in some areas. A study of the spider genus Nesticella suggested that middle Miocene climate change promoted the origin of a subterranean lifestyle in Asian middle latitudes [3]. Here, we emphasized the intrinsic connection between surface creatures and cave organisms, as referred to by Li et al. [48]: S. martensi originated from terrestrial species inhabiting East Asia subtropical EBLFs. During the Miocene, the distribution of surface organisms decreased with the rise of East Asia subtropical EBLFs and the establishment of the monsoon climate [71]. The harsh seasonality likely led to the extinction of surface populations, while local caves with relatively stable temperature and humidity served as ideal refugia for local species, in which they became isolated (the Climate Relict Hypothesis or Pleistocene Effects Model; [4]).

This scenario suggests these cave populations may share a common ancestor followed by multiple independent invasions to different caves or with underground dispersion. Biotic colonization of caves is not a random process but is subject to periods of acceleration and decrease, in association with climate and vegetational changes and the establishment of seasonal climate in subtropical East Asia [48]. The climate-vegetation-relict model proposed for the subtropical East Asian cave biota by Li et al. [48] could well explain the colonization of S. martensi.

The significant intraspecific structure might have resulted in inflation of the true existence of cryptic species because population structure tends to result in more similar gene trees across loci than expected under neutral coalescence. Genetic distance and haplotype data reveal a complex population structure of S. martensi, reflecting low genetic diversity, high genetic differentiation, and varied haplotype number and diversity. Most haplotypes are restricted to a single cave, with the exception of haplotype 11, which is shared among Tangren Cave, Beixin Cave, and Nishui Cave, consistent with phylogenetic results. This pattern matches the high genetic distance and F_ST values. Our results reveal restricted gene flow among S. martensi populations. We speculate that the lack of population expansion among the nine populations of S. martensi may be due to the difficulties in dispersal between cave systems, which isolates the species into independent island populations. This geographical isolation hinders gene flow and promotes genetic differentiation. Additionally, the extreme environmental conditions of caves (dim light, oligotrophic, small daily and annual variations in temperature and humidity) create a relatively closed environment, which further prevents population expansion.

At present, we can conclude that each cave population represents an independent evolutionary lineage, deeply diverged from each other based on our data. We attempt to reconstruct the process of population changes in caves. According to the divergence time tree, a major branch of SYD and BFD formed at 8.94 Mya, and a major branch of TRD, JHD, QBD, NSD, BXD, XD and XGD was formed at 8.60 Mya (95% HPD = 6.60– 10.60 Mya) before present (Fig. 8). During the subsequent evolutionary process, the population sizes of S. martensi from the nine caves remained almost stable until about 0.5 Mya, when they underwent a sequential process of decreasing and increasing (Fig. 10). Bianfu Cave was the first to undergo differentiation and the only one that may have experienced population expansion, with a downward trend in population size from 0.00012 Mya to the present (Fig. 9; Table 3).

Conclusions

Caves are excellent evolutionary laboratories for biological studies, but the highly conserved morphology of S. martensi and the large genetic differentiation, likely due to past geological events, are typical of many short-range endemic taxa, particularly those restricted to caves. Our population structure analysis is based only on COI, a mitochondrial gene. In the future, we hope to use more loci for deeper lineage and geographic research to explain the evolutionary history of S. martensi. Our research indicates that cave organisms have disadvantages in terms of range of movement and gene flow, making caves crucial refuges for such “vulnerable” organisms. Protecting cave environments is imperative; otherwise, when caves are no longer “refuges”, these ancient creatures may not be visible to future generations.