Diversity of Pholcus Spiders (Araneae: Pholcidae) in China’s Lüliang Mountains: An Integrated Morphological and Molecular Approach

Simple Summary Pholcus is the most diverse spider genus in Pholcidae, and is widely distributed in the Palaearctic, Indo-Malayan, Afrotropical, and Australasian Regions. Previously, the Pholcus spiders have not been recorded from the Lüliang Mountains of North China. We undertook an expedition there for the first time. Phylogenetic analyses of DNA sequence data from four gene fragments (COI, H3, wnt, 28S) suggested that Pholcus from the Lüliang Mountains were grouped into nine well-supported clades. We adopted an integrative approach, including morphology and four methods of molecular species delimitation (ABGD, GMYC, bPTP, and BPP), to investigate species boundaries. Such analyses identified the nine clades as nine separate species, of which eight are new to science. All of them belong to the P. phungiformes species group. Abstract Spiders of the genus Pholcus were collected for the first time during an expedition to the Lüliang Mountains in Shanxi Province, North China. Phylogenetic analyses of DNA sequence data from COI, H3, wnt, and 28S genes allowed us to group them into nine well-supported clades. We used morphology and four methods of molecular species delimitation, namely Automatic Barcode Gap Discovery (ABGD), the Generalized Mixed Yule Coalescent (GMYC), Bayesian Poisson Tree Processes (bPTP), and Bayesian Phylogenetics and Phylogeography (BPP), to investigate species boundaries. These integrative taxonomic analyses identified the nine clades as nine distinct species, comprising Pholcus luya Peng & Zhang, 2013 and eight other species new to science: Pholcus jiaocheng sp. nov., Pholcus linfen sp. nov., Pholcus lishi sp. nov., Pholcus luliang sp. nov., Pholcus wenshui sp. nov., Pholcus xiangfen sp. nov., Pholcus xuanzhong sp. nov., and Pholcus zhongyang sp. nov. The species occur in geographic proximity and show many morphological similarities. All of them belong to the P. phungiformes species group. The records from the Lüliang Mountains represent the westernmost distribution limit of this species group.


Introduction
The concept of "biodiversity hotspots" was developed by conservation biologists to highlight areas with exceptional concentrations of endemic species simultaneously experiencing exceptional loss of habitat [1,2]. Although such hotspots are determined mainly by high species richness of vascular plants and vertebrate species, it is assumed that they also harbor correspondingly high diversity of invertebrates [1]. Two such hotspots fall within China's territories: first, the mountains of Southwest China, and second, portions of what is termed "Indo-Burma" which includes parts of Yunnan, Guangxi, Guangdong, and Hainan [2,3]. Even though North China and Northeast China are not considered to be biodiversity hotspots, descriptions follow Huber [21] and Yao et al. [18,28]. This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN. The LSID for this publication is: urn:lsid:zoobank.org:pub:4CAA98A6-2CF3-44C8-94A5-4B363E23B603. The distribution map was generated with ArcGIS 10.2 (ESRI Incorporated Company, Redlands, USA).

Morphological Observation
Specimens were examined and measured with a Leica M205 C stereomicroscope. Left male pedipalps were photographed (unless otherwise indicated in figure legends). Epigynes were photographed before dissection. Vulvae were treated in a 10% warm solution of potassium hydroxide (KOH) to dissolve soft tissues before illustration. Images were captured with a Canon EOS 750D wide zoom digital camera (24.2 megapixels) mounted on the stereomicroscope mentioned above and assembled using Helicon Focus 3.10.3 image stacking software [27]. All measurements are given in millimeters (mm). Leg measurements are shown as: total length (femur, patella, tibia, metatarsus, tarsus). Leg segments were measured on their dorsal side. The specimens studied were preserved in 75% ethanol and deposited in the College of Life Science, Shenyang Normal University (SYNU) in Liaoning, China. Terminology and taxonomic descriptions follow Huber [21] and Yao et al. [18,28]. This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN. The LSID for this publication is: urn:lsid:zoobank.org:pub:4CAA98A6-2CF3-44C8-94A5-4B363E23B603. The distribution map was generated with ArcGIS 10.2 (ESRI Incorporated Company, Redlands, CA, USA).

Phylogenetic Analyses
Genomic DNA extraction and amplification were performed as in Yao et al. [29]. We targeted four DNA fragments for sequencing: the mitochondrial gene fragment encoding COI and three nuclear gene fragments encoding H3, wnt and 28S. Two species Pholcus paralinzhou  and P. taishan Song & Zhu, 1999 were selected as outgroups. Primers are listed in the Supporting Information (Table S1). DNA sequences were checked and edited with BioEdit 7.2.2 [30]. Bayesian inference (BI) and maximum likelihood (ML) methods were used to reconstruct phylogenetic trees using both COI and a combined dataset. BI analysis was performed with MrBayes 3.2.4 [31]. The GTR + I + G model was used for the concatenated data. Two simultaneous runs of four Monte Carlo Markov chains (MCMCs) with default heating parameters were performed for 4 million generations. Trees were sampled every 1000 generations with the first 25% of sampled trees discarded as burn-in. The results were checked using Tracer 1.6 [32] to ensure stationarity. ML analysis was conducted using RAxML 8.2.9 [33] under a GTRCAT model for all partitions, with 500 rapid bootstrap replicates followed by a thorough maximum likelihood tree search. The phylogenetic results using COI and the concatenated data were used to perform the analyses of molecular species delimitation below.

Molecular Species Delimitation
We applied four methods for molecular species delimitation. The Automatic Barcode Gap Discovery (ABGD) online version examines species delimitation with recursive partitioning using a range of prior intraspecific divergence and relative gap widths, estimating the threshold between intra-and interspecific genetic variation to generate species-level groupings. The ABGD analyses were conducted using both Jukes-Cantor and Kimura 2-P distance matrices with options: Pmin = 0.001, Pmax = 0.1, Steps = 10, X = 1.0, Nb bins = 20 [34].
The Bayesian implementation of the Poisson Tree Processes (bPTP) model tests species boundaries based on phylogenetic trees of individual genes. The bPTP method uses nucleotide substitution information and implements a model assuming phylogenetic tree branch lengths are generated by two classes of Poisson processes (intra-and interspecific branching events). This analysis was conducted on a web server (http://species.h-its.org/ ptp/, accessed on 17 January 2023) using individual gene trees. The MCMC was run for 100,000 generations, with a thinning of 100 and burn-in of 0.2 [35].
The Generalized Mixed Yule Coalescent (GMYC) model delimits species from an ultrametric tree of individual genes without prior definitions of species. The GMYC method identifies a time point on the best tree (the tree with the highest likelihood) where the branching rate shifts from speciation to the population coalescent process. This analysis was performed under the singlethreshold model using the R 4.2.2 package SPLITS (Species Limits by Threshold Statistics) [36].
Bayesian Phylogenetics and Phylogeography (BPP) is a multilocus coalescent species delimitation analysis, which requires data from multiple genes and pre-defined candidate species [37,38]. The BPP method accommodates species phylogeny as well as lineage sorting due to ancestral polymorphism and estimates the posterior distribution for different species delimitation models. We used BPP to test apparently conflicting results between the analyses mentioned above. Like Leaché and Fujita [39], we conducted four different sets of analyses with different values of α and β: (i) Gθ(1, 10) and Gτ(1, 10), assuming large ancestral population sizes and deep divergences between species; (ii) Gθ (2,2000) and Gτ (2,2000), assuming small ancestral populations and shallow divergences; (iii) Gθ(1, 10) and Gτ (2,2000), assuming large ancestral populations and shallow divergences; and (iv) Gθ (2,2000) and Gτ(1, 10), assuming small ancestral populations and deep divergences. The analyses were performed using the following settings: species delimitation = 1, algorithm = 0, finetune = 5. The reversible-jump MCMC analyses were run for 100,000 generations and sampled every two generations, with 25,000 samples being discarded as burn-in. Each set of α and β was run at least twice to confirm consistency.

Morphological Variations
A total of nine populations and 42 specimens were studied. Most morphological variation occurs in the shape of the uncus and the presence or absence of its distal apophysis, the presence or absence of the distal apophyses on the procursus, the shape of the vulval pore plate ( Figure 2B), and the presence or absence of the frontal apophysis of the male clypeus. Based on the above characters, the nine populations were identified as nine distinct species. More detailed diagnoses, descriptions, and illustrations are provided under the Taxonomic accounts (see below).

Phylogenetic Relationships
All loci were successfully sequenced for all individuals. A total of 120 sequences of four gene fragments from 40 ingroup members, and four sequences of two gene fragments from two outgroup members were generated. We obtained a concatenated alignment of 2607 bp (COI, 1189 bp; H3, 300 bp; wnt, 292 bp; 28S, 826 bp). The sequences are deposited in the GenBank under accession Nos. OQ706157-706196, OQ719631-719670, OQ719671-719710, OQ719758-719797 (Supporting Information, Table S2). Separate analyses of individual gene COI and concatenated data found compatible topologies (Supporting Information, Figure S1). The BI and ML analyses of the concatenated data supported the same topology. Figure 2A presents the tree from the BI analysis. The Pholcus phungiformes species group was clearly divided into nine well-supported major clades and these clades consisted of samples from a single population, therefore, we defined the nine clades as nine candidate species.

Species Delimitation
The ABGD analysis identified nine provisional species using both Jukes-Cantor and Kimura 2-P distance matrices, and the result was fairly consistent with morphology ( Figure 2A). The ML tree of COI was used to conduct the GMYC and bPTP analyses due to the same topology of the ML and BI trees. The GMYC analysis identified 10 species, with Pholcus linfen sp. nov. (as identified from morphology and ABGD) divided into two species (Figure 2A: number 5, brown and light brown boxes; Supporting Information, Figure S2B). The results from the bPTP were largely consistent with those from morphology and ABGD, except that Pholcus wenshui sp. nov. and Pholcus jiaocheng sp. nov. were recognized as a single species, while Pholcus luliang sp. nov. and Pholcus zhongyang sp. nov. were recognized as another species (Figure 2A: numbers 1 and 2, 3 and 4; Supporting Information, Figure S2A). The BPP analysis requires pre-defined species and phylogenetic relationships of these species. Based on the results of morphology, ABGD, and GMYC, as well as the ML tree of the concatenated data, we used BPP to validate the nine species. The BPP analyses found very high probabilities of speciation events for all of the nodes tested using all four prior combinations. In particular, four prior combinations produced speciation probabilities of one for most of the nodes (Figure 2A: BPP i-iv). These results were also consistent across multiple runs.   Considering all evidence, we conclude that there are nine species from the Lüliang Mountains. The phylogenetic tree derived from the concatenated data clearly divided the samples into nine deeply divergent clades ( Figure 2A). Moreover, the ABGD analysis supports speciation events among the nine species and the result is fairly consistent with morphology.
Although the GMYC analysis divided the species P. linfen sp. nov. (Figure 2A: number 5) into two species, this delimitation result is unreasonable because all four samples (W156-W159) are from the same population. In addition, a single speciation event for P. linfen sp. nov. is well supported by ABGD, bPTP, and BPP analyses. The morphological characters of specimens from the same population are also consistent, e.g., the prolateral membranous process of the procursus with a curved sclerotized apophysis, the procursus with a slightly sclerotized ventro-distal apophysis and without spine-shaped distal apophyses, and the nearly trapezoidal vulval pore plates.
The bPTP collapsed two species pairs each to a single species: P. wenshui sp. nov. and P. jiaocheng sp. nov. (Figure 2A: numbers 1 and 2); P. luliang sp. nov. and P. zhongyang sp. nov. (Figure 2A: numbers 3 and 4). Nevertheless, the other three molecular delimitation results clearly support their status as separate species. Furthermore, P. wenshui sp. nov. can be distinguished morphologically from P. jiaocheng sp. nov. by the prolateral membranous process of the procursus with a strongly sclerotized edge, the procursus with a spine-shaped distal apophysis, the uncus with a wide distal apophysis and a sawtoothed edge, and the epigynal plate strongly curved in the ventral view. P. luliang sp. nov. can be distinguished morphologically from P. zhongyang sp. nov. by the procursus with a spine-shaped distal apophysis and without sclerotized ventro-distal apophyses, the uncus with a slightly curved distal apophysis, and the nearly round vulval pore plates.
Finally, the BPP analyses unambiguously support the speciation events among those nine clades, with the posterior probabilities of most of the nodes are one in four prior combinations. For the two species pairs collapsed by the bPTP, the possibility of P. wenshui sp. nov. and P. jiaocheng sp. nov. recognized as a single species lies in BPP ii and iv, but their posterior probabilities are very low, 0.007 and 0.0007, respectively. Similarly, the possibility of P. luliang sp. nov. and P. zhongyang sp. nov. being delimited to one species was found only in BPP iv, and the posterior probability was only 0.001.
We employed morphology and four commonly used molecular methods for the P. phungiformes species group, and produced fairly consistent results, except for slight deviations arising from GMYC and bPTP analyses. On this basis, we can assert that combining morphological with molecular data allows for rapid and accurate assessment of species richness and therefore such an approach can be considered as an essential part of the conservationist's toolkit. It should however be noted that several molecular-based methods of species delimitation that have been proposed and often applied (e.g., Dincă et al.; Dumas et al.; Li et al. [40][41][42]) have yielded different conclusions. For this reason, the strengths and weakness of each of these methods may still need to be further explored and evaluated.

Comparison of Species Diversity within the P. phungiformes Species Group
Based on the samples collected from the nine populations in the Lüliang Mountains, each is a distinct species. This degree of diversity is comparable proportionally with that found in the Yanshan-Taihang Mountains and the Changbai Mountains. The Yanshan-Taihang Mountains are almost four times the size of the Lüliang Mountains, and they harbor 35 species, almost four times the diversity of the latter [17]. Similarly, both the total area and diversity of the Changbai Mountains [18] are approximately three times those of the Lüliang Mountains. The Taebaek Mountains of the Korean Peninsula are an exception. Their area is just over twice that of the Lüliang Mountains, but with 37 species recorded there [20], its diversity is more than four times that of the latter. Such disparity suggests there may well be species from the group that are yet to be discovered in the three Chinese mountain ranges. This belief is rooted in the fact that all four mountain ranges have similar landforms and ecological niches occupied by the spiders, e.g., rock walls in montane mixed forests. Furthermore, the Lüliang Mountains and the Yanshan-Taihang Mountains are located within the same latitudinal belt as the Taebaek Mountains (35 • to 40 • N).
The Lüliang Mountains probably represent the westernmost distribution limit of spiders of the P. phungiformes species group. None of this group could be found when we sampled specimens of Pholcus in the adjacent Shaanxi Province intensively and extensively in 2013, 2016, and 2019, covering all prefectures throughout the province. They were again absent in the collections from our 2022 expedition to the Qinling Mountains whose range extends from Shaanxi Province to the western part of Henan Province.

Geographic Proximity and Morphological Similarities within the P. phungiformes Species Group from the Lüliang Mountains
Of the nine species from the Lüliang Mountains, eight species are new to science. More significantly, they are in close geographic proximity. For instance, the sister species P. wenshui sp. nov. and P. jiaocheng sp. nov. (numbers 1 and 2 in Figure 1) were 49 km apart. Another two sister species, P. luliang sp. nov. and P. zhongyang sp. nov. (numbers 3 and 4), were found only 42 km apart. The species with the nearest distance in geography are P. linfen sp. nov. and P. xiangfen sp. nov. (numbers 5 and 6). They were found only 39 km apart. Additionally, P. lishi sp. nov. (number 7) was found approximately 49 km away from the nearest species P. zhongyang sp. nov. (number 4), and P. xuanzhong sp. nov. (number 8) was found 48 km away from the nearest P. wenshui sp. nov. (number 1).
All the new species are similar to each other in their morphology. For instance, in males, the procursus is highly complex distally but includes the same three structures: a prolateral membranous process, a dorsal membranous lamella, and dorsal spines (e.g., arrows 1, 3 in Figure 9C, arrows in Figure 9D). Six species, excluding P. linfen sp. nov. and P. xuanzhong sp. nov., possess a spine-shaped or sclerotized distal apophysis on their procursus (e.g., arrow 2 in Figure 9C). Five species, namely P. wenshui sp. nov., P. jiaocheng sp. nov., P. luliang sp. nov., P. zhongyang sp. nov., and P. lishi sp. nov., possess a curved distal apophysis on the uncus (e.g., arrow 1 in Figure 10C); the uncus of the three other species is nearly semicircular or elliptic (e.g., Figure 6C). Five species, namely P. wenshui sp. nov., P. jiaocheng sp. nov., P. luliang sp. nov., P. zhongyang sp. nov., and P. linfen sp. nov., have a small frontal apophysis on their clypeus. In females, the epigynal plate of all the new species except for P. xuanzhong sp. nov. is posteriorly curved.
Such geographic proximity, as well as the morphological similarities, suggests that this group might have undergone a recent radiation. Furthermore, all of these new species are found only on rock walls from the Lüliang Mountains. This group may have become specialized in living on rock walls, although we have not detected any particular adaptive or physiological traits. It is uncertain whether such specialization could have hampered their dispersal and gene interchange among populations between different rock walls in the mountainous region, thereby paving the way for geographic isolation and species radiation. These questions are challenging subjects for further investigation into this species group.
Natural history: The species was found on rock walls. Distribution: China (Shanxi, type locality; Figure 1).