Regional climates shape the biogeographic history of a broadly distributed freshwater crab species complex

The evolutionary importance of paleoclimate regimes has been noted in biogeographic studies. However, little is known about how paleoclimate differences shaped the biogeographic pattern and diversification history of the freshwater fauna in important zoogeographical boundary regions. Here, we aim to investigate how past regional climatic differences have shaped the biogeographic history of the inland aquatic fauna in China using an endemic freshwater crab species complex found on both sides of the Qinling Mountains–Huaihe River Line (QHL), a critical ecological boundary in eastern China, as a model system.

This study focusses on the Qinling-Huaihe Line (QHL) which is a reference line used by geographers to distinguish between northern and southern China at a latitude of about 33° North (where Qinling refers to the Qinling Mountains, and Huaihe refers to the Huaihe River). The QHL is also an important boundary between the Palaearctic and the Oriental biogeographic zones and includes the Qinling-Dabie orogenic belt (QDB), the Yangtze River Plain, and the southern Anhui mountains (Chen, 2004;Cheng, 1987;He, 1994;Tan, 2011;Zhu et al., 2003). This region was significantly influenced by the uplift of the Qinghai-Tibet Plateau, while had a profound impact on the development of the Great River Basin (e.g. the Yangtze and Yellow Rivers, and the Huaihe River in the QHL; Dong et al., 2011Dong et al., , 2015Meng & Zhang, 2000;Shi et al., 2013;Zhang et al., 2015).
The climate in this region is heavily influenced by differences in the intensity of the East Asian monsoon that has created different ecological conditions in either side of the QHL. During the Pleistocene glacial and interglacial periods, the northern part of the QHL tended to be temperate with a cool dry climate, while the southern part of the QHL tended to be hot and wet with either a subtropical or tropical climate, and a temperature difference of about 6-7°C Jin et al., 2010). The uplifted mountain range served to block the flow of the warm and humid southern air mass, resulting in a significantly different levels of precipitation in the south (which became wetter) than in the north (which became drier; Liu, 2005).
The resulting ecological and geological complexity of the QHL today provides an ideal region for testing how ecological changes, Pleistocene climatic oscillations, and topographic variations have affected the biogeographic history of the species found there.
Recent research in the QHL has focused on lineage divergence in non-aquatic species on both sides of the QHL, and its role as a major geographical barrier (e.g. Liu et al., 2017;Song et al., 2016;Yan et al., 2010). In contrast, aquatic species might be expected to have responded differently to past climate fluctuations (i.e. show more genetic structuring) than their non-aquatic counterparts because aquatic habitats have experienced a different biogeographic history and are much less stable than terrestrial systems (Lam et al., 2018;Ribera & Valladares, 2011;Ribera & Vogler, 2000;Wang et al., 2013;Wang et al., 2012). The QHL and the surrounding regions are rich in freshwater systems, and the abundance and complexity of these habitats has driven the geographical divergence of a number of aquatic species. For example, this boundary has played a critical role in the geographical division of the fish fauna in the Palaearctic and Oriental zoogeographic regions (Shaanxi Institute of Zoology et al., 1987). However, little is actually known about the role that past climatic changes on both sides of QHL have played in promoting phylogeographic subdivisions and divergence. The present study focuses on the role of the QHL in shaping the historical diversification of aquatic invertebrates which have been largely overlooked to date.
Sinopotamon Bott, 1967, is endemic to China and is the most speciose and widely distributed genus of freshwater crabs found anywhere in the world (see Chu et al., 2018;Dai, 1999;Ng et al., 2008;mid-Pleistocene, and (b) isolation of different populations of S. yangtsekiense s.l. in a number of separate refuges during the LGM.

K E Y W O R D S
eastern China, freshwater crab, genital divergence, paleoclimate shift, phylogeographic structure, zoogeographical boundary Fang et al., 2013;Shih et al., 2016). For this study we focused on the widely distributed species Sinopotamon yangtsekiense sensu lato, which is found on both sides of the QHL, and exhibits a great deal of variation across its range, and may actually prove to be a species complex (Dai, 1999;Dai & Chen, 1981). For example, morphological variations of the G1 of S. yangtsekiense led Dai (1999) to recognize three subspecies, S. y. yangtsekiense, S. y. tongbaiense, and S. y. shanxianense, based on geographic differences between populations. However, this is not supported by the molecular evidence because these three subspecies do not form a monophyletic group, and the three taxa form two separate lineages (Zheng et al., 2006). In the present study, therefore, we refer to this taxon as S. yangtsekiense s.l. until the taxonomy of this species complex can be formally resolved. Interestingly, our morphological studies indicate that the differences in the morphology of the G1 of these three taxa may be related to past regional climate changes. For example, differences in the size and position of the female vulvae may be an adaptation to slow the evaporation rate of water in response to a drier climate, which has prompted corresponding co-evolutionary changes in the shape of the male genitalia necessary to ensure continued successful mating. The significant variation of the male G1 morphology in populations of S. yangtsekiense found to the north and south of the QHL, and the availability of a comprehensive genetic dataset for a wide range of specimens from this region, make this species complex an ideal research subject to study historical biogeography in these freshwater habitats.
Here, we reconstruct the biogeographical history of the S.

| Population sampling, DNA extraction, and amplification
The taxonomy at the species level used in this study follows Dai (1999) and Ng et al. (2008). A total of 482 specimens of S. yangtsekiense s.l.
(including the three subspecies assigned by Dai, 1999) were collected from 34 localities throughout its distributional range that includes the QHL, the Yangtze Plain, and the Southern Anhui mountains. All samples were collected between 2004 and 2015, stored in 95% ethanol and grouped according to morphological traits and locality ( Figure 1; Table S1). The outgroup taxon, S. depressum, was collected from Kaihua County in Zhejiang Province, China. We sequenced mtDNA markers from 273 specimens, and microsatellite loci from 482 specimens.
The total genomic DNA was extracted from crab muscle or gill tissue using the DNeasy Blood & Tissue Extraction kit (Qiagen ® ).
Haplotype networks were plotted first by the median-joining method (Bandelt et al., 1999) and then by the maximum parsimony Steiner method (Polzin & Daneshmand, 2003) using Network v.5.0.0.1 (www.fluxu s-engin eering.com). The significance of each pairwise comparison was tested with 1000 random replicates, and a standard Bonferroni correction was applied for multiple tests. The relationship between genetic distance and geographical distance was tested using isolation-by-distance (IBD) analysis for all sampling sites using Mantel tests in IBD v.1.53 (Jensen et al., 2005). To verify whether there was significant gene flow between populations, we used the software package Migrate v.4.4.3 (Beerli, 2006(Beerli, , 2009 to estimate the effective number of migrants (M = m/μ, where m is the migration rate and μ is the mutation rate) and effective migrants per generation (Nem).
We generated a relative period for the lineage divergence by COX1 and NAD5. Molecular dating was performed in BEAST v.1.10.4 (Uchard et al., 2018). We conducted two runs of 10,000,000 generations (sample freq = 1,000). After checking for convergence in Tracer v.1.6 (Rambaut & Drummond, 2013), 25% of MCMC were discarded as burn-in. Divergence time was estimated with an uncorrelated lognormal model and a speciation Yule tree priority with a mutation rate substitution rate of (2.0 ± 0.6 stdev) × 10 −9 substitutions per neutral site per year, which is commonly used for the COX1 gene and NAD5 in crustaceans (Jesse et al., 2009). We set the common ancestor of

| Microsatellites data analyses
The software Micro-Checker v.2.2.3 (Van Oosterhout et al., 2004) was used to detect possible null alleles and potential errors in genotypic. Genetix v.4.05 (Belkhir, 2004) and GenePop v.3.4 (Raymond & Rousset, 1995) were used to evaluate the genetic diversity indices, including number of alleles (N), effective number of alleles (Ne), observed heterozygosity (Ho), and expected heterozygosity (He) and to test for Hardy-Weinberg equilibrium (HWE). The significance of the deviation of genotype frequencies from the HWE within populations is expressed by the deviation of F IS from zero based on estimates of multilocus genotypes. These statistics were assessed using FSTAT v.2.9.3 (Goudet, 2001).
Microsatellite data were unambiguously assigned to particular mtDNA lineages by neighbor-joining (NJ) analysis in MEGA v.7.0.26 (Kumar et al., 2016). Genetix v.4.0.5 was used to analyze the genetic variation patterns of eight microsatellite loci by Factorial Correspondence Analysis (FCA) and the genetic diversity values of individual multilocus were calculated, and displayed in threedimensional spatial patterns. The NJ tree of multilocus data was constructed by individual shared allele distances (Du & Lin, 2006). The population genetic structure without an outgroup was evaluated under the Bayesian framework via STRUCTURE v.2.3.4 (Pritchard et al., 2000) with the mixed model and the correlated allele frequency model. The tested K ranged from 1 to 10 independent runs of 1,000,000 MCMC iterations each, the first 100,000 iterations were discarded, and the average operation was performed 10 times.

| Morphological and statistical analysis
The first (G1) and second (G2) gonopods of male brachyurans are typically complex and exhibit considerable morphological variation among taxa (Cumberlidge, 1999;Dai, 1999;Davie et al., 2015;Guinot et al., 2013). This is especially true in freshwater crabs where the unique shape of the G1 is used as a high-weight species-level character in taxonomic identification, along with characters of the G2, and the sixth pleonal segment (A6; Chen et al., 2007;Dai, 1999;Ng, 2008). In this study, we focused on the morphology of the G1, G2, and A6 of males because of difficulties quantifying the characters of the female genitalia.
The terminology used follows Cumberlidge (1999) and the classification follows Dai (1999) and Ng et al. (2008) Table S3 in Supporting Information).
Morphological measurements would be expected to be correlated with the present climate conditions (Begliomini et al., 2017;Sun et al., 2017). We obtained, therefore, 64 climate factors for the sampling sites from 1981 to 2010 from the China Meteorological Data Sharing Service System (http://data.cma.gov.cn/site/index. html; see Table S4 and Appendix S1 for details of climatic factors).
We conducted a one-sample Shapiro-Wilk test and a Kolmogorov-Smirnov test to determine whether data conformed to a normal distribution and the homogeneity of variances, and one-way ANOVAs to compare each trait and factor among populations and groups. In an effort to test whether these morphological characters support distinct differences among populations, we standardized morphological traits and conducted a principal component analysis (PCA) to indicate differences among populations; the linear discriminant function analysis (DFA) was conducted using the lda function in the MASS package v7.3 (Venables & Ripley, 2002); and the clustering algorithm based on Euclidean distance was conducted via R package cluster v.2.1.0 and facteoextra v.1.0.5. The Pearson correlation coefficient was calculated to determine the significant relationships between the morphological variables and climatic factors. All statistical analyses were performed in Program R v. 3.6.1 (R Core Team, 2016).

| Demographic history and ecological niche modeling
To investigate the demographic history of S. yangtsekiense s.l., we used Arlequin to calculate Tajima's D, and Fu's Fs, and significance was evaluated through 10,000 merge simulations. The distribution of all nucleotide mismatches between DNA sequences was accessed using the mismatch distribution analysis, whereby a unimodal distribution indicated that the population had recently experienced rapid growth, and a multimodal distribution indicated a balanced population size (Rogers & Harpending, 1992). The bootstrap method was used to statistically test the validity of the extended model and to calculate the sum of squares (SSD) between the observed and expected distributions (Excoffier & Lischer, 2010). Finally, the demographic history of the population size over time was calculated using the Extended Bayesian Skyline Plot (EBSP; Ho & Shapiro, 2011) implemented in BEAST by assessing the time variation of effective population size. We performed two independent runs for 100,000,000 generations, sampling one tree every 10,000 generations. We applied a unique model of evolution (GTR) and an identical strict substitution rate for all partitions (mean substitution rate value, see above).
Ecological niche modeling (ENM) was performed for the Last Glacial Maximum (LGM) and presented for all populations.

| Phylogenetic analyses and population genetic structure
The network relationship ( Figure 2a

| Morphological variation
Current morphological changes to determine the likely factors that contribute to these variations (the G1 characteristics of lineages are shown in Figure 1). The ANOVA results determined that the five morphological traits and CVs were statistically significant (p ≤ 0.01) among the lineages and populations. Correlation analysis confirmed that morphological traits do not vary with the sample size (p = 0.097~0.945) and provides strong support for the results of the morphological analysis. Two significant PCs found the total variance of the morphological measures to be 87.94% and 12.02%, respectively (see Figure S4 in Supporting Information). Discriminant function analysis showed that individuals were split into two distinct clades with strong support (membership probability >0.95). Overall, TA B L E 1 Analysis of molecular variance (AMOVA) for the mitochondrial data (mt) and microsatellite loci (SSR) in Sinopotamon yangtsekiense s.l. PCA, DFA, and cluster analysis indicated that the partitioning was basically consistent with genetic analysis ( Figure S4). However, the finding that population 8 clustered with Lineage B is doubtful because it is likely that it is an artifact arising from using data from a single individual rather than from a larger sample. All the morphological data exhibited a normal distribution and were, therefore, used for subsequent correlation analyses.

Sum of squares
The influence of precipitation and temperature on morphologi-

| Historical demographic inferences and ecological niche modeling
The mtDNA data for three different datasets (  analyses and the Extended Bayesian skyline plots ( Figure 5). EBSP suggested that demographic expansion occurred approximately 15 kya, which corresponded to the end of the LGM occurring in continental China (Anhuf et al., 2006). Bio9, and Bio18. The predicted distribution pattern for the presentday populations closely matched the actual occurrences of this species ( Figure 6). In addition, the distribution of lineage A during the LGM indicated that the populations then had a significantly smaller range compared to present day levels. For example, during the LGM, the distributional range of the S. yangtsekiense complex was limited mostly to regions of low elevation in the eastern part of its range.

| DISCUSS ION
Our results show that the most basal split within the S. yangtsekiense

| Interspecific relationship and splitting
The QHL represents a natural zoogeographic boundary in China dividing the moist and warm region southeast of the line from the dry and cooler region northwest of the line. This division has increased following a critical paleoclimate transition event-the intensification of the East Asian monsoon (0.85 Ma; Guo et al., 1993;Zhao et al., 2020). Interestingly, the divergence time between the two lineages of the S. yangtsekiense complex was estimated to be 0.78 Ma (95% HPD: 0.48-1.06 Ma), which is largely consistent with sharp climatic changes. This, taken together with our historical demography analyses, infers that the climatic differences on both sides of the QHL may have acted as a phylogeographic barrier, splitting and genetically isolating the existing population, in an event that subsequently led to the evolution of two independent lineages. Interestingly, this dividing line broadly coincides with the Palearctic / Oriental boundary along the QHL (Cheng, 1987;Zhang, 1999).
The unexpectedly close relationship of populations within lineage A that are nevertheless separated by a distributional gap between TBM and YR ( Figure 1) may be attributed to these areas serving as common refugia during the glacial period. This gap between population in TBM and YR might be due to the competitive displacement of S. yangtsekiense by S. honanense which is also found there (Dai & Chen, 1981). It is possible that S. honanense may have caused the removal of the formerly established Nanyang basin habitat of S. yangtsekiense through direct or indirect competitive interactions (Altshuler, 2006;Cheng et al., 2009;Douglas et al., 1994;Kfir, 1997). This phenomenon of competitive exclusion is also found between other closely related species (Altshuler, 2006;Cheng et al., 2009;Douglas et al., 1994;Kfir, 1997) especially where the ranges of two species overlap in hybrid zones (Rieseberg et al., 2000).
The historical river capture events in eastern China could explain the close relationship within lineage B (Li et al., 2005;Ren et al., 1959), where the shared haplotypes and strong gene flow indicate that there were no obvious gene barriers. Historically, the Huaihe and Yangtze Rivers were connected and flowed southward into the Guhong River (Clark et al., 2004;Yang, 2006). Although the exact timing of the individual river connections is controversial, it seems to have occurred at some point in the late Pleistocene and produced significant changes to the drainage network in eastern China (which is consistent with this study). In addition, during the Pleistocene in eastern China numerous freshwater lakes formed, such as Lakes Luoma, Hongze, Gaoyou, and Taihu (Cao, 1989;Zhan, 2008;Zhu et al., 2003). The intermittent connectivity of hydrological features and summer floods produced homogenization that allowed more migration opportunities for populations of S. yangtsekiense in the eastern region of China (Yan et al., 2013;Yang et al., 2009;Yu et al., 2014).

TA B L E 2 Neutrality tests and mismatch analysis (Spatial expansion and Demographic expansion)
Long distance dispersal of freshwater crabs in the past has resulted in the reproductive isolation of populations through reduced gene flow leading eventually to the fixation of novel adaptations and allopatric speciation. In addition, the isolation of small populations of crabs may have led to the fixation of characters that may not be adaptive due to genetic drift and the founder effect (Chung et al., 2014;Wang et al., 2017;Yong et al., 2015). It is possible that the dispersal of the S. yangtsekiense complex has led to significant population fragmentation in parts of its range that were associated with numerous founder events. If this were the case, then the allelic diversity and the number of heterozygotes in the population would both significantly decrease as the number of founder events separating populations increased (Clegg et al., 2002;Cun & Wang, 2010).
We inferred that Qinling-Dabie orogenic belt acted as a dispersal corridor for the S. yangtsekiense complex, and that strong founder effects had occurred in lineage B, and these colonization routes  (Table S1).

| The morphological divergence driven by climatic change
The morphological analysis infers that the S. yangtsekiense complex has been affected by changes in temperature and humidity / precipitation (Figure 4). High humidity and high rainfall are associated with an elongated G1TA, while lower humidity and less rainfall are associated with a curved G1TA and an accompanying shortened G2SS.
Freshwater crabs are semi-terrestrial and are often found out of water where they dig burrows in river banks and forage on nearby land (Dai, 1999;Dai & Chen, 1981). When the habitat of these crabs is subjected to low levels of humidity and precipitation and even prolonged drought, reducing the evaporation rate from the female openings by altering their size and position, is a selective pressure that could lead to morphological variation not only in females, but also in the variation of the G1 shape of the males that must mate with them.
Compared with other non-sexual traits, the variation of genital traits caused by climate changes can lead to reproductive isolation and can promote lineage divergence. We inferred, therefore, that climate change on both sides of the QHL is likely the intrinsic cause of genital divergence in these crabs, which in turn has reduced gene flow and has led to phylogeographic divergence. The genital morphology of males and females has co-evolved, at least in marine crabs (Guinot et al., 2013), where variation in genital traits has prevented hybridization and enhanced reproductive isolation (Brennan & Prum, 2015;Eberhard, 2011). Notably, we inferred that genital coevolution in the S. yangtsekiense complex is consistent with Yao et al. (2020), and appeared to have been subjected to rapid divergence by relation to historical and climatic drive. We proposed that the changes in the genital of the S. yangtsekiense complex is geographically characterized as the northwest is more curved than the southeast. We speculate that the gradual drying of the environment in the northwest of the QHL has resulted in a shift of the size and position of the female genital vulvae towards the midline of the sternum, presumably to reduce the evaporation of water during egg production.
These changes in the size of the opening and the position of the female genital vulvae were accompanied by the co-evolution of the shape of the male genitalia (the G1 TA) which became more angled (bent) to facilitate alignment with the female openings necessary for successful mating.

| Historical demography and glacial refugia
ENM showed that the distributional range of S. yangtsekiense s.l. was significantly constricted during the last glacial maximum (0.78 Ma) and then gradually recovered over time to occupy its present distribution ( Figure 6). The neutrality tests and mismatched distributions indicated that lineages A and B have experienced a recent expansion ( Table 2). The EBSP analysis also showed that S. yangtsekiense s.l. has undergone a rapid demographic expansion after the last glacial maximum (~15 kya; Figure 5), similar to that described for other freshwater species found in this region (Horreo et al., 2018;Li et al., 2014;Lopes-Lima et al., 2018;Salvi et al., 2014;Ling et al., 2012). Tajima's D, however, is not significantly negative, indicating that the population has undergone recent expansion after which it was subdivided, undergoing either an extensive migration and expansion and/or contraction (De Queiroz, 2007;Sharma et al., 2014;Wang & Ge, 2010). Local genetic diversity changes with increasing geographic latitude might be, in part, due to historical contractions or subdivisions of the original population (Yang et al., 2016). Similarly, the relatively low genetic diversity found in several distinct populations (1, 4, 10, 23, and 24) suggests that these may have experienced a bottleneck effect in the past.
Freshwater crabs grow and reproduce best in warm and humid environments, and this provides two explanations for why Pleistocene climate fluctuations have been major factors affecting crab populations in the past. For example, S. yangtsekiense s.l. grows fastest in warm temperatures (between 20 and 28°C), and slowest in cooler temperatures (between 10 and 20°C; Chen et al., 1994). The climate during the last glacial maximum 18-21 kya contrasts sharply with the climate of the past 17,000 years and with the most recent glacial period (Cullen, 1981). According to a series reports, the temperature during the last glacial period in eastern China was between 4 and 12°C (Shi et al., 1997;Xu et al., 2010;Zhen et al., 2014), which is cold enough to reduce the growth rate and increase the mortality rate. In addition, cold temperatures reduce the fertility rate in freshwater crabs because females lay fewer eggs than normal (only 30-100), and females have even been observed to eat their own eggs when subjected to cold conditions and food shortages (Chen et al., 1994;Dai, 1999).
According to the EBSP and mismatch distribution analyses, populations of S. yangtsekiense s.l. could have undergone a recent population expansion across the QDB from west to east ( Figure 5), followed by bottleneck events associated with population contractions (Rogers & Harpending, 1992). Our study provides corroborative evidence for cumulative founder effect events during the eastward expansion.
The QDB mountain range provided numerous glacial refuges for freshwater species during the past ice ages (Hewitt, 2000;Zhang et al., 2008;Zhao et al., 2020). The results of our genetic diversity analysis confirmed the presence of at least three glacial refugia (YR, TBM, and HMR) within the QDB during the LGM. The population expansion of S. yangtsekiense s.l. from refugia in the QDB during the late Pleistocene is evidenced by intensified genetic divergence during this time. This scenario was similar to that reported for the recent colonization of the QBD by S. acutum (Fang et al., 2015). The climatic impact of the glacial period on the QHL in China was different from its impact on Western Europe and North America. In the latter two regions, significant areas were covered by ice sheets during the glacial period, forcing many species to migrate to warmer areas to the south. When temperatures warmed up during the interglacial periods, species recolonized the northern areas (Hewitt, 2004). In contrast, during the glacial period, the QHL in China was not significantly glaciated, but the change in climate altered the ecosystems there dramatically, leaving only a few sheltered pockets (refugia) within the QDB and the Huangshan region (HMR) that remained environmentally suitable for some species where more tolerable conditions persisted. Interestingly, we found that these refugia (e.g. TBM, YR) either coincided with, or were adjacent to, the existing distributional range of populations of freshwater crabs. The high annual precipitation rates in the highlands of the QDB (such as the Qinling and Tongbai mountains) during the glacial period established local humid conditions where populations of S. yangtsekiense s.l. could survive (despite their inland location; Ju et al., 2007;Knutti et al., 2004). A similar pattern of species distributions coinciding with refugia in the QBD was also found in studies of amphibians (Zhang et al., 2008) and reptiles (Ding et al., 2011). Dai and Chen (1981)  The high divergence of the lineages A and B in the S. yangtsekiense complex, and the congruence of genetic and morphological data, suggest that this complex probably contains a cryptic species. It is likely that the existence of cryptic diversity in the S. yangtsekiense complex has been underestimated, given the rapid diversification within Sinopotamon (since 0.8 Ma; Ji et al., 2016). We have not revised the taxonomic status of the S. yangtsekiense complex here because this is beyond the scope of the present work whose focus is on the biogeography of this taxon in relation to historical climate changes. Clearly, further studies are needed to clarify the unstable taxonomy of the S. yangtsekiense complex.

ACK N OWLED G EM ENTS
We thank Prof. Kaiya Zhou (Nanjing Normal University) for invaluable comments on this study. We are also grateful to Prof. Hongxin

DATA AVA I L A B I L I T Y S TAT E M E N T
The COX1