A comprehensive phylogeography of the widespread pond snail genus Radix revealed restricted colonization due to niche conservatism

Abstract To clarify the effect of niche conservatism on evolutionary history, we focused on freshwater snails, which have different ecological and phylogenetic properties from previously tested taxa. We conducted a phylogenetic analysis using 750 lymnaeid individuals from 357 sites of eleven Radix species. Then, we estimated the ancestral distribution using the geographic coordinates and colonization routes. In addition, a statistical test of the colonization distances in the latitudinal and longitudinal directions was performed. We also conducted ecological niche modeling for two widely distributed species using climatic data. Ancestral geographic reconstruction estimated the origin of the genus to be around the Indian subcontinental region and showed that latitudinal immigration distances were shorter than longitudinal immigration distances in the diversification process. Ecological niche models suggested that the current distribution was restricted by climate, with annual mean temperature and precipitation of the driest month as particularly strong factors. Niche conservatism to the climate can affect the diversification of freshwater snails.


| INTRODUC TI ON
Whether or not ecological traits have restricted current biodiversity patterns is a critical issue. In particular, niche conservatism (NC), namely the retention of niche-related ecological traits through evolutionary history (see Wiens et al., 2010), is an attractive theme that bridges ecology and evolutionary biology (Peterson et al., 1999;Wiens et al., 2010;Wiens & Donoghue, 2004;Wiens & Graham, 2005). In a biogeographical context, NC, especially as a process that leads to current biodiversity patterns, is often suggested to be important in limiting distribution transitions (Losos, 2008;Pearman et al., 2008;Wiens & Graham, 2005). In particular, NC could explain the latitudinal diversity gradient (LDG), which is one of the principal themes of historical biogeography (Fine, 2015;Hillebrand, 2004;Mannion et al., 2014;Mittelbach et al., 2007;Pontarp et al., 2019).
In the context of the LDG, NC is often recognized as the basis of the tropical conservatism hypothesis (TCH; Pianka, 1966;Wiens & Donoghue, 2004;Wiens & Graham, 2005). Organisms in tropical regions that occupied greater areas in the past (~34 Ma; Bowen et al., 2015;Cronin, 2010) have restricted evolutionary transitions over different climatic regions due to NC, leading to the LDG (Farrell et al., 1992;Peterson et al., 1999;Wiens & Donoghue, 2004). Recently, the TCH has been evaluated as a potential driver of the LDG in many phylogeographic studies (e.g., Buckley et al., 2010;Duchêne & Cardillo, 2015;Economo et al., 2018Economo et al., , 2019Kerkhoff et al., 2014;Owens et al., 2017). Moreover, NC can also drive diversity gradients from a species origin outside the tropics, and NC has universal importance for explaining biodiversity patterns .
In contrast to the clear evidence for the contribution of NC to the diversity patterns of many taxa, the patterns of some taxa are often explained by other mechanisms such as competition, carrying capacities, and the places of the origin of organisms (Pyron & Wiens, 2013;Ramos Pereira & Palmeirim, 2013;Rolland et al., 2014;Siqueira et al., 2016). More importantly, studies comparing a large number of taxonomic groups have concluded that the contribution of NC to the LDG was fairly limited (Boucher-Lalonde et al., 2015;Jansson et al., 2013). Although phylogenetic studies are useful in understanding the importance of NC Wiens & Donoghue, 2004;Wiens & Graham, 2005), most of these studies have been conducted with vertebrates, plants, and insects due to limited taxon sampling (Fine, 2015;Jablonski et al., 2017). However, as ecological differences such as the dispersal mode of each taxon can interactively alter the effect of NC on the evolutionary process (Ackerly, 2003;Eiserhardt et al., 2013;Kubota et al., 2017), the relative importance of NC in evolutionary history can differ depending on the taxon considered (Chiu et al., 2020). While further studies on various taxa with different ecological traits may be useful (Chiu et al., 2020), the importance and contribution of NC to the diversification process are poorly understood, except for some taxa (e.g., Amphibians: Kozak & Wiens, 2010;Reptiles: Pyron & Burbrink, 2009;Stephens & Wiens, 2009).
Here, we focused on the common freshwater snail genus Radix to address this issue. Radix belongs to the family Lymnaeidae (Gastropoda) and has a wide distribution range covering most of the Old World, from the tropics to the subarctic regions (Aksenova et al., 2018). In addition, Radix species has high morphological diversity including phenotypic plasticity in their shell (Pfenninger et al., 2006;Terry & Duda, 2021;Ward et al., 1997), and thus many morphological species have been described (Vinarski et al., 2020). However, many species are currently considered to be synonyms, and nine species and one undescribed species have been detected based on the molecular taxonomy (Aksenova et al., 2018). In general, freshwater snails are often dispersed by water currents, and can rarely be transported by wind and animals such as birds, amphibians, insects, and mammals (Bespalaya et al., 2020;Kappes & Haase, 2012;van Leeuwen et al., 2012;Rees, 1965;Walther et al., 2008). As a result, they have low active and high passive dispersal potential (Kappes & Haase, 2012), and often have a wide distribution range. Perhaps as a consequence, freshwater snails seem to have relatively high niche flexibility, including niche shifts within the species (Cordellier & Pfenninger, 2008Kisdi, 2002;Torres et al., 2018). These ecological characteristics may mean that colonization is unlikely to be strongly hindered by causes other than climatic factors (e.g., the significant lack of dispersal ability or the habitat dependence on interspecies relationships such as symbiosis) and may suggest that the distribution transitions of freshwater snails could be based on more purely probabilistic processes than that of well-studied vertebrates having low dispersal ability such as amphibians and reptiles. This is quite different from the taxa traditionally used to examine NC and diversification. The effects of NC on the historical diversification of freshwater snails have been largely unexamined; therefore, the verification of the effect of NC using freshwater snails could provide new insights into the contribution of NC to diversification.
In this study, we collected comprehensive Radix materials and sequences from their whole distribution area and conducted phylogenetic analyses to estimate the historical distribution transition.
We took into account the phylogeographical structure within each K E Y W O R D S dispersal, freshwater snail, latitudinal diversity gradient, Lymnaeidae, Mollusca, passive disperser, tropical conservatism hypothesis species, which many conventional NC studies using phylogenetic approaches have not fully addressed. In addition, we quantitatively evaluated the latitudinal and longitudinal (L/L) dispersal for each distribution transition event by reconstructing the historical distribution change under specific latitudes and longitudes. The distribution transitions of each lineage have previously been evaluated as either categorical data or indicator values (e.g., Duchêne & Cardillo, 2015;Kerkhoff et al., 2014). These evaluations do not adequately assess colonization within the same climate zone, even though the climate is continuous. In contrast, our approach may be able to accurately estimate the L/L dispersal distances for each colonization event within the same climate zone. In other words, more colonization events are expected to have a longer dispersal distance in the longitudinal direction than in the latitudinal direction if the adaptation to the new niches of climatic temperature restricts the distribution transition.
Furthermore, we also performed ecological niche modeling to identify the climatic characteristics restricting the colonizations. Finally, we integrated these approaches and validated the contribution of NC as a process, focusing on the transitions in the distribution area.

| Sampling and DNA methods
We collected data from 750 lymnaeid individuals from 357 sites across nine countries and regions-278 of the individuals were sampled by us and the remaining data were obtained from GenBank ( Figure 1 and Table S1). The L/L information of samples from GenBank was identified using published papers and registered information from GenBank (Table S1). When L/L information was not available, it was assigned based on the location name of the sampling site (e.g., the specific lake, village, or township) ( Figure 1). This could be done accurately for most samples, as detailed site names were usually provided. These materials covered the greater part of the Radix distribution area.

| Phylogenetic methods
We conducted a phylogenetic analysis using CO1 sequences (278 newly collected by us and 472 from GenBank). Eight species were selected as outgroups from the eight most closely related genera with reference to a previous phylogenetic study (Aksenova et al., 2018). These sequences were aligned with MUSCLE v3.8 (Edgar, 2004), and there was no gap in the alignment. CO1 phylogenetic trees were determined using both the Bayesian inference (BI) and maximum likelihood (ML) methods. Before both analyses, the same sequences were stacked using FaBox1.41 (Villesen, 2007), and as a result, 455 haplotypes were detected in the genus Radix (Table   S1). Next, we selected the appropriate partitioned models of sequence evolution using PartitionFinder2 (Lanfear et al., 2016; Table   S3). Based on these models, the BI analysis was performed using MrBayes5d version 3.1.2.2012.12.13 (Tanabe, 2012a), an extended software of MrBayes v3.1.2 (Ronquist & Huelsenbeck, 2003), with two simultaneous runs. We discarded non-convergence trees after examining convergence and effective sample size (ESS; larger than 200) using Tracer v. 1.6 (Rambaut et al., 2013) and the remaining samples were used to estimate phylogeny (for detailed settings, see Table S3). ML analysis was performed using IQ-TREE version 1.6.7 (Nguyen et al., 2015) and the evolutionary model was selected for IQ-TREE with the -spp option (for detailed settings and the model, see Table S3). For the ML analysis, we assessed nodal support by performing ultrafast bootstrapping (Hoang et al., 2018)  We reconstructed both nuclear (H3+ITS2+28S) and combined (CO1+H3+ITS2+28S) phylogenies to assess phylogenetic positions among each species estimated by CO1 phylogeny. We used 81 individuals from 9 Radix species and 83 individuals from 11 Radix species for nuclear phylogeny and combined phylogeny, respectively (Table   S1). Racesina luteola was selected as an outgroup for both phylogenies. To eliminate the uncertainty of the ITS2 and 28S alignments, trimAl 1.2 (Capella-Gutiérrez et al., 2009) was used for subsequent phylogenetic analyses. All phylogenetic analyses were performed using the same approach as that used for the CO1 phylogeny (for detailed settings, see Table S3).

| Geographic reconstruction
To trace the historical changes in the distribution of taxa, we performed Bayesian ancestral distribution reconstruction using BayesTraits v.3.0.1 under a geographical model (Pagel, 2004(Pagel, , 2017. A Bayesian CO1 tree was used after removing outgroups and using L/L information (Table S1). The branch lengths were scaled to have a mean of 0.1 with reference to the BayesTraitsV3 manual (Maede & Pagel, 2016). Then, reconstruction was estimated under the following settings: ngen = 15,000,000, sample freq = 1000, burnin = 5,000,000, and rate dev = autotune. Then, the estimated L/L sites on each principal node were denoted using QGIS. Furthermore, to evaluate historical L/L colonization, the L/L immigration distances between pairs of sites in the sequenced nodes were calculated using the GRS80 model. Finally, we conducted an exact Wilcoxon signedrank test to examine the difference between the obtained distances of L/L immigration using the package "exactRankTests" in R version 3. 5. 1 (R Core Team, 2018; Torsten & Kurt, 2019).

| Ecological niche modeling
To clarify whether climatic variables with latitudinal gradients actually restrict the Radix distribution, we conducted ecological niche modeling (ENM) using MaxEnt version 3.4.1 (Phillips et al., 2004(Phillips et al., , 2006. We performed ENM with two species with a relatively F I G U R E 1 The map of localities where the materials of examined Radix spp. were collected. Solid lines denote sites obtained: the latitudinal and longitudinal (L/L) information from sampling data, published references, or Genbank registration data. Dashed lines denote sites estimated L/L information based on given locality names. This map was generated from Global Self-consistent, Hierarchical, Highresolution Geography Database, Version 2.3.5 (Wessel & Smith, 1996, 2016, using QGIS Version 2.18 (QGIS Development Team, 2016). See Table S1 for further information large sample size, R. auricularia and R. plicatula. The potentially dispersible landmass of the species was presumed to be −20 to 180° latitude and 0-90° longitude. First, 19 bioclimatic variables with a spatial resolution of 2.5 arcmin were obtained from WorldClim 2 (Fick & Hijmans, 2017), as proxies for the ecological traits of Radix.
To obtain a parsimonious and interpretable model (Merow et al., 2013), we eliminated spatially correlated bioclimatic variables of current climate data prior to the modeling using ENMTools (Warren et al., 2010(Warren et al., , 2019 in R (R Core Team, 2018). As MaxEnt is relatively robust for high collinearity (Elith et al., 2011), we only removed the variable pair with a correlation coefficient greater than |0.80|. In the modeling process, we used only linear and quadratic features, and the regularization multiplier was set to 3 to avoid overfitting (Merow et al., 2013;Radosavljevic & Anderson, 2014;Syfert et al., 2013).
In the pilot analyses, this setting had little effect on the area under the receiver operating characteristic (ROC) curve, which represents the predictive accuracy. MaxEnt runs were conducted under the following settings: maximum number of background points = 10,000, duplicate presence records = remove, maximum iterations = 5000, and output format = Cloglog. Models were evaluated using fivefold cross-validation and the AUC as an indicator. Moreover, to determine which environmental variables were important drivers of the distribution, we used three different approaches (Phillips, 2017).
First, we determined the percent contribution based on the increase or decrease in regularized gain under model construction with MaxEnt. Second, we used permutation importance based on the drop in training AUC when the values of that variable on training presence and background data were randomly permuted for each variable. Third, we used the jackknife test, that is, each variable was deleted sequentially and a model was created with the remaining variables. Next, a model was created with each variable in isolation and each corresponding model was compared. These approaches are considered popular metrics for the evaluation of variable contributions (Bradie & Leung, 2016).

| Phylogenetic relationship of Radix
Our mitochondrial phylogenetic analysis included 11 species  (Hoang et al., 2018;Nguyen et al., 2015). In this tree, these values of relatively well-supported branches are only shown. See Table S1 for further information. Each Radix species was determined based on a previous study (Aksenova et al., 2018). (b) The map of distribution areas of each Radix species. This map is illustrated in the Mollweide projection to show the correct area. This map is generated from Global Self-consistent, Hierarchical, High-resolution Geography Database, Version 2.3.5 (Wessel & Smith, 1996, 2016, using QGIS Version 2.18 (QGIS Development Team, 2016). For further information, see also Figure 1 and Table S1 (a)

(b)
distribution area, whereas other species had limited distribution areas ( Figure 1b and Table 1).

| Geographic reconstruction
Our geographic reconstruction estimated the ancestral locations of each node ( Table 2). The results suggest that the origin of Radix

| Ecological niche modeling
Nine bioclimatic variables remained after removing those that were strongly correlated, namely Bio1: annual mean temperature, Bio2:

F I G U R E 3
The Bayesian phylogenetic tree of the Radix inferred from the combined dataset (1945 bp). Racesina luteola was selected as an outgroup. Each operational taxonomic unit label represents the material ID. Scale bar indicates substitutions per site. Numbers on the branches indicate the Bayesian posterior probabilities (BPP) and the Maximum likelihood ultrafast bootstrapping value by IQ-TREE (Hoang et al., 2018;Nguyen et al., 2015). These values of low supported (BPP < 0.070) and terminal branches are not shown. See Table S1

| DISCUSS ION
Our geographical reconstruction estimated that the origin of Radix  (Dayrat, 2005;e.g., Aksenova et al., 2017;Bolotov et al., 2014Bolotov et al., , 2017Vinarski et al., 2016) is needed to clarify the actual status of these subclades; however, these lineages could be recognized as evolutionary units.
Radix auricularia, which was distributed in the most northern area of the 11 species, had the widest distribution area; R. natalensis, R.
Our geographic reconstruction showed distribution transitions throughout the evolutionary history of Radix (Figure 4). Most of the colonization routes had a longer immigration distance in the longitudinal direction than in the latitudinal direction ( Figure 4 and  (Barton, 2011;Clarke, 2003). Further researches, especially direct quantification of ecological traits of Radix, are needed to determine which mechanism is principal for the evolutionary history of Radix.
In any case, this indicates that the difficulty in climatic adaptation (i.e., NC) may have determined the distribution transitions in the evolutionary history of Radix. Many phylogeographical studies have shown that distribution transitions over climate zones are difficult (e.g., Economo et al., 2018Economo et al., , 2019Kerkhoff et al., 2014;Kozak & Wiens, 2010;Stephens & Wiens, 2009). In particular, many animal taxa with high active dispersal abilities have been found to have restricted historical transitions in their distribution areas due to climate factors (e.g., Bats: Buckley et al., 2010;Stevens, 2011;Birds: Duchêne & Cardillo, 2015;Hawkins et al., 2006;Butterflies: Hawkins & DeVries, 2009;Owens et al., 2017;Flies: Löwenberg-Neto et al., 2011). In contrast, few studies have focused on animal taxa with passive dispersals. However, our results showed that climate was also an important factor in determining the distribution of passive dispersers, which have a more probabilistic dispersal mode. In plants, which are critically different from animals in terms of ecology and genetics, many studies have shown the importance of climate to distribution (e.g., Kerkhoff et al., 2014). This may be due to the fact that many plants are passive dispersers. Regardless of the dispersal mode, colonization seems to be restricted by latitudinal climate when the distribution of an organism expands by dispersal.
Furthermore, our ecological niche models showed that the current distribution ranges of Radix species were restricted by climate ( Figure 5 and Figure S4.3). Based on our modeling, which had a high prediction accuracy despite the limited climate variables used, the annual mean F I G U R E 5 Maps show the occurrence probabilities of two Radix species by MaxEnt (Phillips et al., 2004(Phillips et al., , 2006. Maps for Radix plicatula are shown in a, and maps for R. auricularia are shown in (b). The probabilities were produced by a complementary log-log (cloglog) transform, which is considered to be given a stronger theoretical justification than the logistic transform . The average of replicate runs of the probabilities are shown in color on a per grid temperature (Bio1) and precipitation of driest month (Bio14) were relatively influential factors for both of the species analyzed under all the modeling approaches used (Table 4 and Figure S4.3). Air temperature is often a critical factor in determining distribution ranges (Calosi et al., 2010;Merriam, 1984), and it is also important for freshwater organisms (Calosi et al., 2010;Cordellier et al., 2012) because shallow water temperature is strongly correlated with air temperature. Moreover, precipitation of the driest month can be an important factor because droughts affect the population dynamics of freshwater mollusks that inhabit temporal inland water (Gérard, 2001;Woolhouse, 1992). The restriction to the distribution ranges, which caused by these climate variables, was clear; however, the potential distribution ranges seemingly had fairly wide ( Figure 5). This fits well with the high niche flexibility within the species (Cordellier & Pfenninger, 2009;Torres et al., 2018).
Our analyses showed two suggestions about the evolutionary history of Radix: their colonization is more likely to occur in the longitude direction than in the latitude direction, and one of the import-

ACK N OWLED G M ENTS
We are grateful to S. Uchida for suggestive advice on geographic analyses and supporting sample collection. We thank S. Isao, K.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
All obtained sequences were deposited in the GenBank