Connectivity and selfing drives population genetic structure in a patchy landscape: a comparative approach of four co-occurring freshwater snail species

The distribution of neutral genetic variation in subdivided populations is driven by the interplay between genetic drift, migration, local extinction and colonization. The inﬂuence of environmental and demographic factors has also been increasingly examined in empirical studies, but generally focusing on a single species. An open question is whether these factors will similarly, or idiosyncratically, aﬀect a guild of species occupying the same, though exhibiting diﬀerent traits, mating systems and histories. We addressed it by comparing the population genetic structure of the four most common species of hermaphroditic freshwater snails in Guadeloupe (Lesser Antilles), which occupy a network of patchily distributed sites experiencing temporal variation in water availability. We analyzed microsatellite variability in 21 to 43 populations per species, and built predictions on how several environmental and demographic variables, quantiﬁed from a long-term annual survey, as well as species traits, may aﬀect the distribution of genetic variation. These species displayed similarities, such as fairly high levels of variation, but with marked diﬀerences among sites, as well as strong genetic diﬀerentiation and limited isolation by distance, which can be explained by passive dispersal (strong role of site connectivity), extinction/colonization dynamics and variation in local population size. They also exhibit diﬀerences, largely due to the mating system with less genetic diversity and more genetic diﬀerentiation in the two selﬁng species when compared to the two outcrossing ones. These diﬀerences can also be attributed to interspeciﬁc interactions resulting from the ongoing invasion of Guadeloupe by one of the species studied, which aﬀects the demography of other species, and, to a limited extent, to local environmental factors. Our comparative approach shows both differences and uniqueness in the way species occupy the same landscape, and provides a possible entry to interspeciﬁc interactions in community assembly. a paired t -test to check for temporal homogeneity. The other variables were not considered, either because they do not vary in time (density of favorable habitats), or because the comparison is irrelevant (e.g., invasion age). All variables were compared between the populations studied for genetic variation and all other populations from the long-term survey using linear models. We calculated the correlation between pairs of variables over all populations studied using Pe arson’s correlation coefficient. We evaluated the effect of environmental and site variables on genetic variability ( R A in P. acuta and D. depressissimum and both R A and R G in A. marmorata and D. surinamense ) using linear


Introduction
(e.g., monogamous vs. polygamous reproduction) in shaping population genetic structure. For example, connectivity and habitat patchiness play a role in structuring population structure in orthopterans from central Spain (Ortego et al., 2015), while turtle species studied in the USA were sensitive to water availability (Reid et al., 2017; see Table S1). Note that only a fraction of these studies considers closely related species, and even fewer species sets occupying the same landscape.
Freshwater habitats, especially lentic ones (e.g., ponds), are appropriate systems for conducting multispecies (comparative) analysis of population genetic structure, because the available habitats are clearly delimited from the inhospitable matrix although temporal variation occurs as a function of precipitation / drought regimes, which may affect both population size and migration (Chase, 2007;Heino et al., 2015;Jarne & Delay, 1991). Whether dispersal is passive or active of course also affects migration. For example, snails disperse essentially passively while arthropods may disperse actively in the adult stage or as resting forms (see examples in Table S1). We focus here on a guild of freshwater snails (Hygrophila) occupying a patchy landscape in Guadeloupe, Lesser Antilles to address the influence of the three factors mentioned above on population genetic structure. This biological system has been well characterized through a long-term survey of snail distribution on an island scale, as well as population genetics, demographic and ecological studies (Chapuis et al., 2017;Dubart et al., 2019;Lamy et al., 2012;Lamy et al., 2013a;Pantel et al., in press). Two important characteristics of our study should be highlighted. First, it is conducted at the scale of the Grande-Terre island, meaning that migration takes place essentially at this scale, with limited 'natural' input from outside Grande-Terre, except through occasional introductions, essentially through aquarium trade or birds (see e.g. Duggan, 2010). Second, we focused on the four most common snail species, with reasonably similar occupancy at metapopulation scale (50-70% of occupied sites), certainly because sampling is easier, but also because genetic and demographic parameters can more precisely be estimated than in rare species (Heino et al., 2015).
We focused on four species, and characterized microsatellite variation in 21 to 43 populations per species. This allowed us to consider the influence of the three categories of factors mentioned above (environment, traits and demography), which are known to influence the distribution of genetic variation in snails in general (e.g., Jarne & Delay, 1991;Escobar et al., 2011;Figuerola & Green, 2002;Van Leeuwen et al., 2013), and in the Guadeloupe system in particular (Chapuis et al., 2017;Lamy et al., 2012). First, the environment (and landscape) was characterized based on several parameters, such as site size and connectivity, used in previous studies (Lamy et al., 2012;Dubart et al., 2019;Pantel et al., in press). Second, two species are mainly outcrossing, while the two others are mainly self-fertilizing. The four species also differ in their aestivation capacity when sites dry out, which occurs seasonally during the dry season in a fraction of sites. Third, they exhibit demographic differences, for example in extinction and colonization rates, but also similarities (e.g., similar metapopulation occupancy). Moreover, Physa acuta is an exotic species in the process of rapidly invading the whole Grande-Terre, while the other three are native and rather stable in terms of metapopulation occupancy. In order to characterize how these species react similarly vs. idiosyncratically to these factors (a comparative approach), we predicted how each factor should influence the distribution of genetic variation (e.g., site size should be associated with larger or more accessible populations and therefore more variation). Four species is certainly not a sufficient number to draw conclusions on the factors driving species differences in genetic structure in general. However, we suggest conducting more comparative studies, up to now rather rare, to evaluate the generality of these factors. We also discuss how such a comparative approach should help to progress towards a more integrated eco-evolutionary perspective on community building and species interaction (Hendry, 2017;Vellend, 2016).

Long-term survey of snail communities
The biological system, surveys and species studied have been described in papers spanning the last five decades (see Pointier & Combes, 1976;Pointier & David, 2004;Pointier, 2008;Chapuis et al., 2017;Dubart et al., 2019;Lamy et al., 2012;Lamy et al., 2013a;Pantel et al., in press), and we provide here the relevant information for the current analysis. We initiated in 2001 a long-term yearly study of the guild of freshwater snails from the geographically close islands of Grande-Terre and Marie-Galante (Guadeloupe archipelago, Lesser Antilles). Grande-Terre is a plateau of about 800 km 2 harboring a network of ponds (ca. 2000) of various sizes, and a few small rivers (mostly temporary) and swamp grasslands connected to mangroves (back-mangroves). Marie-Galante is a smaller island (160 km 2 ) with a similar topography and pond network. Water availability varies widely in time with the highest water level at the end of the rainy season (August to November or December). A fraction of sites completely dry out either yearly, or more irregularly, and remains dry for up to several months. This long-term survey includes almost 280 sites, 245 of which are in Grande-Terre (Lamy et al., 2012;Lamy et al., 2013a). Field sampling took place in January-February when these communities peak in abundance. They include 24 mollusc species, which makes up the major part of the macrobenthos. Each site was explored by three persons for at least 15 min at each visit (total searching time: 45 min). The presence and density of snails were estimated by foraging both the sediment and plant strata using a scoop, and rock surfaces or floating debris were visually surveyed. Snail density was estimated on a semilog scale (no individuals, < 1 ind. / m 2 , 1 to 5 ind. / m 2 , 5 to 10 ind. / m 2 …). Several environmental parameters were characterized per site (detailed in Table S2) including size, aquatic vegetation cover, water connectivity to neighboring watersheds, hydrological regime, and density of favorable habitats. Site stability was defined based on the first four variables.

The species studied and population sampling
Our analysis bears on two pairs of closely related species: P. acuta and Aplexa marmorata from the Physidae family, and Drepanotrema depressissimum and Drepanotrema surinamense from the Planorbidae family (Morgan et al., 2002;Wethington and Lydeard, 2007). In each pair, species are ecologically and morphologically similar, but do not hybridize. These four species are the most common ones in the Guadeloupe metacommunity, occupying ca. 50-70% of sites on average, and co-occur in a large number of sites (Dubart et al., 2019;Pantel et al., in press;Pointier & David, 2004). They differ in mating system, since P. acuta and D. depressissimum are outcrossers, while A. marmorata and D. surinamense are selfers (Escobar et al., 2011). The two Drepanotrema species also resist better in dry sites than the two Physids, aestivating under rocks or vegetation (Pointier & Combes, 1976). They also exhibit demographic differences, for example in rates of site colonization and extinction (Dubart et al., 2019;Pantel et al., in press). Moreover, P. acuta originates from the Eastern-central USA and has markedly increased its range over the last two centuries to now be cosmopolitan (Bousset et al., 2014;Ebbs et al., 2018;Vinarsky, 2017). We observed it in two sites from Grande-Terre in 2001, but this figure increased to ca. 60% of sites in 2017. The three other species are neotropical and are considered native in the Lesser Antilles (Pointier 2008). The two Physid species have been shown to interact in Guadeloupe, most likely through direct reciprocal competition, with a depressing effect of P. acuta on the demography of A. marmorata, but a predicted stable co-existence (Dubart et al., 2019) and a character displacement (life-history traits) towards more 'rstrategy' in A. marmorata within sites, following the arrival of P. acuta (Chapuis et al., 2017). Demography was characterized based on two parameters (Table S2). The first is 'long-term population size', based on both snail density and site size. The second refers to the time in years since the first detection of the invasive P. acuta and sampling for the genetic analysis -referred to as 'invasion age' in this species and 'time since P. acuta invasion' in the other species (Table S2). We considered two classes, not or recently invaded by P. acuta vs. invaded for more than two years, as in Chapuis et al. (2017). In A. marmorata and Drepanotrema spp., populations of the first class are unlikely to show a detectable evolutionary response to the presence of P. acuta because of the limited number of generations (Chapuis et al., 2017).
The four species were sampled in Grande-Terre during the yearly field surveys (January-February) in 74 sites in which one to three of the four species were co-occurring (Figure 1, Tables S3 and S4; 1.57 species per site). Four sites from Basse-Terre and Marie-Galante were also considered in the study of P. acuta (three sites) and A. marmorata (one site). Physa acuta, A. marmorata, D. depressissimum and D. surinamense were sampled in 21, 43, 25 and 27 sites respectively, from 2007 to 2010. Individuals were killed just after sampling in 80ºC water for 1 minute and preserved in 95º ethanol. The genetic analysis was conducted on at least 10 individuals per population for a total of 459 individuals in P. acuta (mean per population ± s.d.: 21.86 ± 5.51), 781 in A. marmorata (18.14 ± 6.42), 749 in D. depressissimum (29.96 ± 2.73), and 365 in D. surinamense (13.52 ± 6.3). Less than 10 individuals were available in three populations (Table S4).  Table S3. The four sites sampled in Basse Terre of Guadeloupe and Marie-Galante (Table S3) are not represented. 0.1 min (geographic coordinates) = ~ 10 km.
Philippe Jarne et al. Table 1. Predictions on the influence of mating system, site (environmental) characteristics and demographic factors on population genetic diversity and structure in the four species studied. The rationale for these predictions is given in the main text. Although predictions for several characteristics are identical (e.g., size and connectivity), they are presented separately for the sake of clarity. The effects are given on genetic diversity within populations (Div) and differentiation among populations (Diff). Predicted effect: positive (+), negative (-), limited, and unclear (?).

Physa acuta (outcrosser)
Aplexa marmorata (selfer) Note that the population genetic data used in this paper have partly been published: a full analysis of population structure in D. depressissimum in Guadeloupe is reported in Lamy et al. (2012) and population genetic richness of both D. depressissimum and A. marmorata is used in Lamy et al. (2013b). However, we will present all the genetic data for the sake of clarity, and refer to previous work when required.

Genetic diversity and differentiation
The analyses were conducted in the same way in the four species, except when otherwise stated. We computed the number of polymorphic loci, allelic richness (rarefied to the smallest sample size per species (R A , Petit et al., 1998;8, 10, 14 and 7 diploid individuals in P. acuta, A. marmorata, D. surinamense and D. depressissimum respectively) per population and species using FSTAT 2.9.3.2 (Goudet, 2001). Gene diversity (unbiased H e ;Nei, 1987) and genotypic disequilibrium for locus pairs (Weir, 1996) were computed using GENEPOP 4.4 (Rousset, 2008). In the two selfing species (A. marmorata and D. surinamense), we also estimated the genotypic diversity (R G ; i.e. the number of genotypes) rarified to sample size of 10 and 7 respectively. We calculated the correlation between R A and H e in each species, as well as with R G in the selfing species, using Pearson's correlation coefficient.
Deviations from Hardy-Weinberg equilibrium (HWE) were tested at each locus using exact tests (Raymond & Rousset, 1995), and a global P-value for all loci was obtained using Fisher's exact test and a Markov Chain algorithm with GENEPOP 4.4. A multilocus estimate of Wright's inbreeding coefficient (F IS ) was estimated following Weir and Cockerham (1984). The selfing rate (S) was calculated based on David et al. (2007)'s method, which is insensitive to null alleles, an acute issue in outcrossing species, and less of a problem in selfing species since null alleles are readily detected when occurring in the homozygous state (see . S was here derived from the g 2 method using RMES (David et al., 2007).
Genetic differentiation between population pairs and overall differentiation was assessed in each species based on exact tests using GENEPOP 4.4, and their magnitude was estimated based on the estimator θ of F ST , which was calculated following Weir and Cockerham (1984). Hereafter, we refer to the estimates between population pairs as pairwise F ST . We also estimated the genetic differentiation between each local population and all other sampled populations, hereafter site-specific F ST , using the Bayesian method implemented in GESTE 2.0 (Foll & Gaggiotti, 2006). We calculated Pearson's correlation coefficient between site-specific F ST and genetic diversity (H e ) in each species. The geographic distribution of genetic diversity was evaluated first by testing for isolation by distance in each species based on pairwise distance matrices (Euclidean geographical distance and F ST /(1-F ST )) using a Mantel test (1000 permutations) in Vegan R package (Oksanen et al., 2015;R Development Core Team, 2005 version 3.6.3). We also examined visually the pairwise F ST matrix after they have been restructured based on a seriation approach implemented with seriation in R (R Development Core Team, 2005) -the occurrence of population clusters would be evidenced by population blocks with lower intra-block F ST than with other populations. This also allows building cladograms to better visualize population relationships using the "complete" method of hclust (R Development Core Team, 2005). We also used a more formal clustering approach, sNMF, proposed by Frichot et al. (2014). This method reconstitutes population clusters based on ancestry coefficients. This method is not sensitive to strong inbreeding, and was previously successfully used in the highly selfing snail Biomphalaria pfeifferi (Tian-Bi et al., 2019) -this is the reason why we preferred it to popular approaches such as STRUCTURE (Pritchard et al., 2000). sNMF was implemented using the R package LEA (Frichot & François, 2015;R Development Core Team, 2005). We used a target number of clusters varying from 1 to 10, a regularization parameter of 50, 50 runs were performed per analysis, and several analyses were performed per species. Based on the optimal number of clusters, we drew barplots of average ancestry coefficients per population, and reported these values as piecharts on the Grande-Terre maps using the gdal R package (R Development Core Team, 2005).

Drivers of genetic diversity
We assessed the influence of several site characteristics, the mating system and demographic factors (described above; see Table S2) on allelic richness (R A ) and on genetic differentiation (site-specific F ST ) in the four species, but also on genotypic richness (R G ) in the selfing A. marmorata and D. surinamense. The site characteristics were size, vegetation cover, connectivity, stability and density of favorable habitats, and the demographic factors were invasion age and long-term population size. Given that our first genetic sampling was performed in 2007, we used the value per site of these variables over the 2001-2007 period (averaged for those variables for which we got a value per year; Table S5). The spatial distribution of two variables, averaged over years, is provided in Table S1. Because environmental data were missing, the following populations were excluded from the analyses: Sainte Rose (Basse-Terre) for P. acuta, Ravine Vinty (Basse-Terre) for A. marmorata and Pointe des Chateaux (Grande-Terre) for D. depressissimum.
We built general predictions on genetic diversity and structure, based on (i) environmental parameters, (ii) the mating system, and (iii) demographic parameters (Table 1). These predictions are qualitative: we indeed know that the mating system should have a strong effect, but the precise relationship between environmental parameters and effective population size is not known. Moreover, all these parameters may interact, and the effect of each one may not easily be single out (see Discussion). (i) The predictions depend on the model of population genetic structure underlying species dynamics, especially on how the environmental variables are connected to species demography. However, we expect that site size, site connectivity, and density of favorable habitats should contribute positively to intrapopulation variability and negatively to differentiation among populations, because they positively affect effective population size and gene flow. Stability might have the same effect as the first three parameters in the two physids. The two Drepanotrema species are dwelling well in unstable sites, because of their aestivating capacity (Pointier & Combes, 1976), which should limit (perhaps annihilate) the effect. These predictions can be refined on the basis of previously obtained estimates of extinction and colonization rates per species at metapopulation scale (Dubart et al., 2019;Pantel et al., in press). In most cases, these environmental parameters did not affect these rates, and when they did, this was in the direction of the predictions above (Table S6). For example, site size positively affects colonization in three species, which should contribute positively to genetic diversity. Vegetation cover should also have the same effect in the four species, because of its negative influence on extinction and its positive influence on colonization in both P. acuta and D. surinamense (Table S6). (ii) We expect less diversity and more divergence in selfing than in outcrossing populations, as detailed in Introduction. Note also that given the near-absence of variation in selfing rates within species, the mating system was not considered in the statistical models presented below. (iii) Invasion age is predicted to affect the two Physa species. In the invasive P. acuta, young populations should be less diverse than old ones reflecting temporary founder effects at the colonization front. Both demographic population size and colonization rates are negatively affected in A. marmorata by the competition with P. acuta (Chapuis et al., 2017;Dubart et al., 2019), which should depress diversity in sites that have been invaded by the latter for a longer time. The expectation is less clear in the two Drepanotrema species, since the arrival of P. acuta only slightly affects their density (P = 0.051 and 0.038 respectively vs. 1.2*10 -4 in A. marmorata; see Chapuis et al., 2017 for details). Long-term population size should have the same effect as site size.
We compared size, vegetation cover, connectivity and stability over all sites studied between the 2001-2007 and 2008-2015 periods using a paired t-test to check for temporal homogeneity. The other variables were not considered, either because they do not vary in time (density of favorable habitats), or because the comparison is irrelevant (e.g., invasion age). All variables were compared between the populations studied for genetic variation and all other populations from the long-term survey using linear models. We calculated the correlation between pairs of variables over all populations studied using Pearson's correlation coefficient. We evaluated the effect of environmental and site variables on genetic variability (R A in P. acuta and D. depressissimum and both R A and R G in A. marmorata and D. surinamense) using linear models. We first built models with the five site variables (size, vegetation cover, connectivity, stability and density of favorable habitats), and then added invasion age and long-term population size in turn in order to avoid overfitting. The most parsimonious combination of explanatory variables was retained based on AIC (Akaike, 1974). We then looked for interactions between these variables (when more than one was significant), but no interaction was significant (see Results). Explanatory and response variables were standardized in order to obtain comparable regression coefficients in the linear models. All statistical analyses were performed using R (R Development Core Team, 2005).

Genetic diversity within populations and species
Summary results on genetic diversity in the four species are reported in Table 2, and full results in Table  S4. For the sake of clarity, we tried, as much as possible, to present the results in the same order: P. acuta, A. marmorata, D. depressissimum and D. surinamense. Much more variation was detected in the two outcrossing species (especially in D. depressissimum, due to highly variable tetranucleotide loci), with all loci polymorphic in all populations, than in the two selfing species. Some populations of these two species even showed no variation at all, i.e. exhibited a single, fully homozygous, multilocus genotype. The average allelic richness (R A ) per population was 3.287 (minimum: 1.982 / maximum: 4.426), 1.588 (1.000 / 2.555), 11.052 (6.260 / 14.160) and 2.039 (1.000 / 3.193) in P. acuta, A. marmorata, D. depressissimum and D. surinamense respectively. As a comparison, the mean number of alleles per species, averaged over loci, was 9.29, 6.00, 37.20, and 9.08 respectively. Genotypic diversity ranged from 1.000 to 8.816 in A. marmorata, and from 1. Significant deviations from HWE (P < 0.05) were observed in most populations of P. acuta (18 out of 21), and heterozygote deficiencies were low to moderate (mean F IS = 0.159; Tables 2 and S4A). However, the selfing rates estimated with RMES were not significantly different from zero (mean S(g 2 ) = 0.030), suggesting that the significant heterozygote deficiencies are due to null alleles, as reported previously in this species (Bousset et al., 2004(Bousset et al., , 2014Escobar et al., 2008). The situation was very different in the other Physid species, A. marmorata: most tests of departures from Hardy-Weinberg equilibrium were highly significant (Table S4B). The estimates of F IS were accordingly very high (0.970; Table 2), and the selfing rates as well (0.984). The selfing rates estimated with RMES were high and significantly different from zero in the subset of populations in which it was possible to successfully run the analysis (Table S4B), and similar to those derived from F IS . Results in D. depressissimum are similar to those in P. acuta: we detected low F IS values (mean = 0.036) with only four populations in which Hardy-Weinberg proportions were rejected at P < 10 -3 (Lamy et al., 2012). The selfing rates estimated with RMES were similarly very low in D. depressissimum (mean S(g 2 ) = 0.012), significantly departing from zero in a single population out of 25 (Table S4C). In D. surinamense, as in A. marmorata, most tests of departures from Hardy-Weinberg equilibrium were highly significant (Table S4D) and estimates of both F IS and selfing rates were very high (0.922 and 0.950 respectively; Table 2). The selfing rates estimated with RMES were also high, and consistent, except in two populations, with those derived from F IS (Table S4D).

Population genetic differentiation
Highly significant genetic differentiation was detected between most population pairs in P. acuta and the two Drepanotrema species (Table S7). In both P. acuta and D. surinamense, a single P-value was higher than 0.01, and none in D. depressissimum. The situation is somewhat different in A. marmorata with 142 tests out of 841 with P > 0.05 (Table S7). Fourteen populations out of 43 were involved in more than 10 such tests (out of 42 per population), especially Favreau, Grand Maison, and Sainte Rose Est (18, 22, and 17 such results). Interestingly, eight of these fourteen populations are located in the Eastern part of Grande-Terre (Figure 1). The average estimate of pairwise genetic differentiation (θ) between P. acuta populations was 0.213 (range: 0.007 / 0.533; Table 2 Table S7. The average site-specific F ST was 0.249 (range: 0.100 / 0.562) in P. acuta, 0.547 (0.215 / 0.872) in A. marmorata, 0.071 (0.012 / 0.243) in D. depressissimum, and 0.506 (0.121 / 0.954) in D. surinamense, and are not very different from the mean pairwise values. The average genetic diversity predicted quite accurately the site-specific differentiation, with more scattered data in A. marmorata than in the other species (Table 2, Figure 2).
Mantel test revealed significant isolation by distance in P. acuta populations (r = 0.710, P = 0.001; Figure  2), even when populations from Marie-Galante were excluded (r = 0.301, P = 0.011), but not in A. marmorata (r = 0.037, P = 0.236; r = 0.074, P = 0.130 when excluding the single population outside Grande-Terre). Mantel test was also significant in D. depressissimum (r = 0.318, P = 0.020), but not in D. surinamense (r = -0.077, P = 0.857; Figure 2). The seriation and clustering approach did not allow detecting clear patterns of population differentiation (results not shown). The few indications are consistent with the results of the more quantitative sNMF approach, and we therefore rely on the latter, considering each species in turn. In P. acuta, the optimal number of ancestral populations is four ( Figure S2A). The two populations from Marie-Galante belong to a single cluster (V2). One cluster (V1) is largely represented in the North of Grande-Terre, and one (V3) in the South-West, while V4 is widespread. A lower optimal number of clusters was detected in A. marmorata (K = 2; Figure S2B), with the smallest cluster rather distributed in the South-West of Grande-Terre, but not only. A single cluster was detected in D. depressissimum, suggestive, together with the limited isolation by distance, of an even distribution of genetic diversity at Grande-Terre scale. The pattern is much more complex in D. surinamense with eight clusters ( Figure S2C), few populations (e.g., Belin or Poucet) belonging to a single cluster and no obvious population structure. Moreover, this analysis did not suggest common geographic patterns of clustering among species. Even when enforcing a fixed number of clusters for all species (e.g., K = 2, as in A. marmorata), no common pattern emerged (data not shown).

Environmental and demographic drivers of genetic diversity
The mean values of site characteristics (Table S4) for the 2001-2007 period did not differ significantly from those of the 2008-2015 period in the same set of sites (Table S8), suggesting temporal stability over our survey. It should be added that no temporal variation in cumulative rainfall over both the rainy season and the little rainy season has been detected over the same period (Pantel et al., in press). We also compared values for these sites with those of all other sites in our survey over the 2001-2007 period, and detected no difference, except for stability: the genotyped populations were in more stable sites than the others (Table S8), presumably because unstable sites are often dry by definition, limiting the possibility of sampling. With this exception, the sites in which genetic analyses were conducted therefore seem representative of all sites in our survey. We conducted the same two series of tests for each species separately. No differences were detected (all P > 0.05; results not shown), except when comparing the sites studied in A. marmorata and D. depressissimum and all other sites for vegetation cover (t = -3.039, P = 0.004 and t = -4.009, P = 3*10 -3 respectively): the sites sampled for genetic analyses showed less vegetation cover. We detected correlations between some of these variables when considering all genotyped populations (Table S9). Unsurprisingly, stability was positively correlated with size and vegetation cover, and connectivity was positively correlated with invasion age. We observed two unexpected significant correlations, namely between size and connectivity (negative) and between stability and the density of favorable habitats (positive). Note also that none of the variables was correlated to all others, suggesting that they bring partially independent information and can therefore be used in the statistical models below.
In both P. acuta and D. depressissimum (the two outcrossing species), all descriptors of genetic diversity and differentiation (R A , H E and site-specific F ST ) were highly correlated, F ST being negatively correlated to the other two descriptors (R 2 of linear regression of F ST on H E was 0.882 in P. acuta, and 0.909 in D. depressissimum, and the corresponding value for the regression of F ST on R A was 0.926 and 0.910; Figure 2 and Table S10). The same was true for R G , H E and F ST in A. marmorata and D. surinamense (R 2 of F ST on H E was 0.529 in A. marmorata, and 0.784 in D. surinamense, and the corresponding value for the regression of F ST on R G was 0.539 and 0.530; Figure 2 and Table S10). When two species co-occurred in a large-enough number of populations, we tested for correlations between species pairs for these descriptors: the only significant correlation was found for P. acuta and A. marmorata (R 2 was between 0.300 and 0.326 for the three descriptors; N = 11 sites).
We analyzed how genetic diversity and differentiation was affected by the environmental and population variables (Table S2) using linear models (results in Tables 3 and S11). The results were not qualitatively affected by considering, or not, the populations outside Grande-Terre in P. acuta and A. marmorata. Site size, stability and density of favorable habitats had no influence. Connectivity was the most influential variable with a positive effect on genetic diversity in three species (marginally significant for R A in P. acuta; Figure 3), except in D. surinamense, and a negative one on population differentiation in A. marmorata and D. depressissimum. Vegetation cover had a significant positive effect on diversity, and a negative one, on differentiation in D. depressissimum. Allelic and genotypic diversity were higher in A. marmorata populations recently or never invaded by P. acuta than in those that have been invaded for longer. However, invasion age had no effect on P. acuta itself. Long-term population size had (surprisingly) a positive effect on differentiation in A. marmorata. It also had the same effect than connectivity in D. depressissimum, but that on differentiation was no more significant when vegetation cover and connectivity were included in the model. These results can be compared to the expectations proposed in Table 1, and most were not fulfilled with a majority of non-significant effects. However, significant effects (15 situations) went in the expected direction, except that of long-term population size on differentiation in A. marmorata. Table 3. Effect of environmental and population variables on genetic diversity and differentiation in the four species studied. Results are reported for R A and site-specific F ST in P. acuta and D. depressissimum, and for R A , R G and F ST in A. marmorata and D. surinamense. Results are from linear models (see text for explanation), and the best models were evaluated based on Akaike Information Criterion. The significant (positive or negative) effects are indicated by + or -(P < 0.05), and by ++ or --(P < 0.01). NS indicates no effect. The P-value for connectivity in P. acuta was indeed 0.08, and 0.07 for invasion age in A. marmorata. Full models and effect sizes are reported in Appendix Table S11. All significant effects were in the Direction predicted in  Figure 3. Relationship between allelic richness and connectivity in the four species studied. Each point is a population. The black lines are simple linear regressions (with 95% SE as grey overlay). Connectivity is measured on an increasing scale (from 0 to 3).

Discussion
A comparative multispecies approach to population genetic structure The most classical approach for deciphering the neutral population genetic structure has been, and remains, to study several populations of a single species sampled at a given time, and to infer population processes based on signatures in molecular markers (Charlesworth & Charlesworth, 2010;Lewontin, 1974). These signatures have been in many cases related to environmental variables, traits (e.g., the mating system), and to spatial distribution and population demography, using more and more sophisticated statistical approaches (Hanski & Gaggiotti, 2004;Foll & Gaggiotti, 2006;Balkenhol et al., 2015;Manel & Holderegger, 2013). We here analyzed genetic variation in several species at once, a comparative approach that remains relatively rare (Table S1). The four species studied display rather high levels of genetic variation within populations, but with strong variance among populations, suggesting a variable effect of genetic drift within sites, as already observed at similar spatial scale in P. acuta (Bousset et al., 2004;Escobar et al., 2008). The extremely high variation in D. depressissimum is largely due to highly variable (tetranucleotide) loci (Lamy et al. 2012). In the two other species, we do not know of previous studies to which our values can be compared. The marked differentiation among populations is associated in the four species with no to limited isolation by distance and an absence of clear geographic patterns. Moreover, most of genetic differentiation, especially in the two selfing species, is due to differences in withinpopulation diversity, irrespective of geographic position, and therefore probably driven by populationspecific factors. We can therefore suggest that the genetic structure of the four species at Guadeloupe scale corresponds more to an asymmetric island or a metapopulation model than to other models (e.g., a single panmictic population), as may be expected for freshwater snails occupying ponds (Jarne & Delay, 1991;Lamy et al., 2012). In the former models, we indeed expect marked variation in genetic diversity among populations and no marked geographic structure. Note, however, that cautiousness is required with P. acuta since quasi-equilibrium has probably not been reached because of its expanding dynamics. We previously show that the asymmetric island model should be favored over a metapopulation model in D. depressissimum, based on both spatial and temporal genetic data (Lamy et al. 2012). Temporal data are not available in the three other species, preventing further distinction in models of population structurenoting that the conditions under which the metapopulation model actually applies are rather narrow (Fronhoefer et al., 2012). In our study, we also observed large differences in genetic diversity and differentiation, lower and higher respectively in A. marmorata and D. surinamense than in P. acuta and D. depressissimum. These similarities and differences will be discussed below with regard to connectivity and the mating system, the two variables with the strongest effects, and then to other environmental parameters and demography.

Connectivity is a major driver of the distribution of genetic variability
We proposed in Table 1 how several environmental factors may affect genetic variation following the logic of Lamy et al. (2012). We found essentially no effect in P. acuta and D. surinamense, and two to three significant effects in the two other species. The most marked effect is that of connectivity, detected in three species, affecting both genetic diversity and population differentiation. Our definition of connectivity is based on local landscape structure, i.e. site accessibility linked to local topography, and therefore falls in the 'structural connectivity' category (see e.g., Storfer et al., 2007). In a comparative study of freshwater fish and invertebrates from Australia, structural connectivity corresponds to realized connectivity, estimated with genetic markers, in 50% of the 110 studies considered (Hughes et al., 2013). It can therefore be considered as a rather good indicator of the distribution of genetic variation, as observed in our study. The central effect of connectivity on genetic variation was already highlighted by Lamy et al. (2012) in D. depressissimum in Guadeloupe, but also by Bousset et al. (2004) and Escobar et al. (2008) in southern France populations of P. acuta. Sites with high connectivity index are indeed connected during the rainy season, and the effective size of local populations is probably much enlarged by these connections. It is perhaps less obvious why D. surinamense is not affected by connectivity. This species might disperse less efficiently through flood or more actively through bird transportation. However, this should apply to D. depressissimum, since the two species have similar ecology and aptitude to colonize empty sites and an extinction/colonization dynamics driven by the same environmental factors (Pantel et al., in press). Interestingly, the colonization rate is not affected by connectivity, except in P. acuta, even if it is higher in physids than in planorbids (Table S6; neither is extinction with no difference among species). This suggests that the effect of connectivity is not mediated by the extinction/colonization dynamics -we already noticed that the latter has little influence on the distribution of genetic variability in D. depressissimum (Lamy et al., 2012). Moreover, there was no correlation in genetic diversity among species, except between the two physids. Overall, these results might be explained by the fact that dispersal is mostly passive (through water currents or birds), and leads to colonization or migration events that are both occasional and occurs at random with regard to both species status and geographic scale.

The impact of the mating system
Our analysis confirms that the observed genetic patterns are consistent with an outcrossing mating system in P. acuta and D. depressissimum and a selfing one in A. marmorata and D. surinamense, with limited variation in selfing rate within species, as previously recognized (Escobar et al., 2011;Lamy et al., 2012). Genetic diversity is extremely similar in the two selfers, and comparable to what has been observed in other selfing snails (Charbonnel et al., 2002;Lounnas et al., 2017;Tian Bi et al., 2019). The results on population structure (F ST ) mirror those obtained with within-population variation, unsurprisingly given the strong correlations between within-and among-population variation in the four species, especially in D. depressissimum. We observed both lower diversity and higher differentiation in the two selfing species than in the two outcrossing ones. This was expected because of reduced effective population size and increased selective interference in selfers compared to outcrossers (review in Charlesworth, 2003;Ingvarsson, 2002;Otto, 2021). Interestingly, if we assume a neutral island model and an effective population size divided by two in selfers compared to outcrossers, based on the differentiation measured in the selfing species, F ST in the outcrossers should be that observed in P. acuta (see below the role of invasion in this species). Even if this reasoning is only approximate, the mating system alone might explain this difference in genetic differentiation. However, we observed strikingly large intrapopulation diversity and low differentiation in D. depressissimum. Both could be due to a higher mutation rate, consistent with the extremely large number of alleles detected, but also to a higher regional demographic size. This latter can be estimated as the product of the number of occupied sites multiplied by the number of individuals per site. The first does not differ markedly among the four species with 50-70% of occupied sites per species (Pantel et al., in press). The second can be approximated through long-term population size (Table S2), which, however, was lower in D. depressissimum than in the other species. Demographic size therefore does not seem to be an explanation to the higher diversity / lower differentiation, except if the 'genetic size / demographic size' ratio is much higher in D. depressissimum. One cause might be aestivation since this species survives very well in dry sites (Lamy et al., 2013a), but this is also true of D. surinamense (Pantel et al., in press). We also note that long-term population size did not affect genetic variation within species, suggesting that this estimator of demographic size does not capture adequately effective population size. As mentioned above, the four species differ in their colonization rate (not in extinction rates), but the difference is between the two families, not between selfers and outcrossers (Pantel et al., in press). As a consequence, it is likely that demographic extinction and colonization do not affect genetic variation strongly, as already shown in D. depressissimum (Lamy et al., 2012), but rather modulate the influence of the mating system. Overall, this latter seems to strongly influence genetic variation in the four species studied, due to its effects on effective population size.

Which other factors affect genetic variation?
We tested the influence of factors other than the mating system and connectivity with mostly identical expectations in the four species (Table 1). Most are not fulfilled (no effect), but when they are, it is generally in the right direction (Table 3). Vegetation cover shows the expected effect in D. depressissimum only. It should in principle favors large population size; it indeed runs against population extinction in all species, except P. acuta (no effect; Pantel et al., in press), and favors colonization in D. surinamense. A possible explanation is that the 'vegetation cover' variable aggregates several plant and algae species (Pantel et al., in press), and that the four species react to different types of vegetation. It is also somewhat unexpected that site size, site stability and the density of favorable habitats did not affect genetic variation, all the more that stability is positively correlated with both site size and density of favorable habitats. Site size might be a poor indicator of population size, as suggested by the fact that long-term population size, which beyond size also includes population density, is connected to genetic variation in two species out of four. Stability is perhaps of limited impact, because the three native species studied have developed strategies allowing coping with drought (Pointier & Combes, 1976;Lamy et al., 2013a). The absence of effect of density of favorable habitats (an indicator of possible immigration at the scale of a few km) is consistent with the weak spatial patterns, both estimated through isolation by distance and the clustering analysis -P. acuta is an exception, but this might be due to the invasive dynamics. Overall, this is suggestive of a limited role for passive dispersal through birds compared to passive dispersal through water flow, given the strong effect of connectivity (the two variables are not correlated).
The distribution of genetic variation may also be affected by species invasion, especially in the invasive species P. acuta. In terms of site occupancy at landscape scale, P. acuta has been spreading since our first survey in 2000, and still shows an invasive dynamic in Grande-Terre (Dubart et al., 2019;Pantel et al., in press). This dynamic might be associated with variation loss (Dlugosch & Parker 2008;Estoup et al., 2016;Wares et al., 2005;), but limited loss has been observed in invasive populations of this species at worldwide scale (Bousset et al., 2014;Ebbs et al., 2018). Moreover, the average variability of Guadeloupe populations (e.g., H e = 0.530) is similar to that observed in the two studies conducted at similar scale (Bousset et al., 2004;Escobar et al., 2008). It therefore seems that the invasion of P. acuta in Guadeloupe is ancient enough -the species presence has been noted in some restricted areas since the early 1970's (JP Pointier, pers. obs.), and involved a large number of propagules as perhaps the result of several invasion origins. Theory also predicts a loss of variation at invasion front as a consequence of bottlenecks (Excoffier & Ray, 2008;Nei, 1975). However, we did not detect an effect of invasion age on genetic variability (except perhaps in the two populations of Marie-Galante, which are poorly variable), suggesting that site colonization involves large number of propagules.
A last point is whether this invasion may affect genetic variation in native species. The two physid species indeed showed competitive interactions through modified extinction and colonization when the other species is present (Dubart et al., 2019). Moreover, Chapuis et al. (2017) showed that, although the number of sites occupied by A. marmorata is stable, population density has strongly been affected by the invasion of P. acuta. Moreover, this study highlighted an evolution of life-history traits in A. marmorata in response to this invasion. This was mirrored in the present study by a loss of neutral genetic variation in A. marmorata. Such an effect was not detected in the two Drepanotrema species, although their population density is slightly affected by the invasion of P. acuta. This suggests that the effect might occur in Drepanotrema with some delay. Invasive species may have a negative effect on population densities of native species in introduction areas, which in turn may affect the dynamics of genes and induce rapid evolution (e.g., Strauss et al., 2006;Stuart et al., 2014). A rapid shift in life-history traits has previously been documented in A. marmorata as a result of the invasion of P. acuta (Chapuis et al., 2017). Whether the reduction in neutral genetic variation documented here is related to this shift remains an open question, but documenting such decreases in other bioinvasions might be a first approach to analyzing their evolutionary impact.
Overall, our work indicates that both common and idiosyncratic effects of the environment, space and interspecific interactions on genetic variability can be detected, calling for detailed studies. For example, Haag, Riek, Hottinger, Pajunen, & Ebert (2005) showed that genotypic richness increases with population age in two Daphnia species, but that other factors (e.g., presence of competitor species) did not. Other examples are mentioned in Introduction and can be found in Table S1.

Conclusion and perspectives
The comparative approach developed here leads towards a more general perspective on population genetic structure than proposed in studies focusing on single species. Indeed, it allows detecting common and idiosyncratic effects of the environment, space and interspecific interactions on the distribution of genetic variation in sets of species occupying the same space with limited migration from outside. The four species studied here displayed striking similarities in genetic structure, with marked variation in genetic diversity among sites and strong genetic differentiation, due to passive dispersal, extinction/colonization dynamics and erratic variation in local population size. Differences are due to the mating system that amplifies local drift in selfing populations, interspecific interactions, and to a very limited extent local environmental factors. We also detected little geographic coherence in the distribution of genetic variation among species. In studies conducted at the same (not too large) geographic scale (Table S1), the effect of traits on genetic variation was detected in some cases (e.g., lifespan in Reid et al., 2017), but landscape structure (e.g., habitat fragmentation in Ortego et al., 2015) and demographic factors (e.g., population age in Haag et al., 2006) also play a role. However, none of these studies, with the possible exception of the last two studies cited, carry out a comparative analysis to detect common genetic patterns and the role of idiosyncratic vs. common influences among species. We suggest extending these comparative approaches, moving towards comparative landscape genetics (see e.g. Burkhart et al., 2016;Ortego et al., 2015).
We also suggest to include in such studies the impact of interspecific interactions, as we did here and in previous genetic and demographic studies (Chapuis et al., 2017;Dubart et al., 2019;Lamy et al., 2013a). Our approach constitutes a first entry into the processes shaping these interactions, i.e. upscaling from population ecology to community ecology. The similarity between forces acting at species and community level has long been recognized (e.g. Bell, 2001;Hubbell, 2001), particularly when considering reasonably similar species occupying the same matrix (i.e. a metacommunity), and is at the heart of community genetics (Agrawal, 2003;Vellend, 2010;Vellend & Geber, 2005). Studying the relationships between species and genetic diversity and assessing how they are affected by environmental parameters and space is one way of highlighting these similarities Vellend & Geber, 2005). In the metacommunity of Guadeloupe, this approach highlighted (again) the influence of connectivity (Lamy et al., 2013b;Lamy et al., 2017). Taking into account eco-evolutionary dynamics would be an alternative pathat this stage they have essentially been studied at the intraspecific level (Hendry, 2017). The evolution of life-history traits and the demographic decline of A. marmorata following the invasion by P. acuta is an example of such eco-evolutionary dynamics (Chapuis et al., 2017). Integrative approaches, including longterm genetic, demographic and ecological studies, have been carried out in other species groups, such as guppies (Bassar et al., 2017) or anoles (Pringle et al., 2019) in the Caribbean, to name just a few well-known examples, thus providing an insight into the links between population biology and community ecology. They should be continued and generalized to other model systems, including more complex communities.