Patterns of richness of freshwater mollusks from Chile: predictions of its distribution based on null models

Faunistic similarity patterns

To determine the faunistic similarity patterns, a cluster analysis was done on a presence and absence matrix, from which a similarity dendrogram was built, using the Jaccard Index coefficient of similarity based on presence/absence records of the molluscs taxa. For this analysis, we used the unweighted pair/group method with arithmetic averages (UPGMA) as agglomeration algorithm (Sokal & Rohlf, 1995) through PRIMER-6 software (Clarke & Gorley, 2006).

Analysis of the hierarchical cluster and SIMPROF analysis (similarity profile permutation test) were carried out to evaluate the biogeographic zones of freshwater molluscs based on one-degree bands. The SIMPROF permutation test (5% level) was used to determine the clusters with a significant internal structure, using 50,000 permutations (Clarke, Somerfield & Gorley, 2008).

Rapoport’s rule and geometric restrictions

To evaluate Rapoport’s rule effect, the extension of the latitudinal range was calculated as the difference between the maximum and minimum latitude of the distribution range of each species and the mean calculated as the mean among these values (Gaston, Blackburn & Spicer, 1998). This relation was evaluated by making a regression analysis with a randomization approximation (50,000 random matrices), implemented in the EcoSim 7.71 software (Gotelli & Entsminger, 2004). When the relation between both variables is positive, the hypothesis that the geographical range of the species increases at higher latitudes is upheld, reflecting the possible adaptations of the tolerant species to extreme climatic conditions. Additionally, to explain the distribution of the species richness, the mid-domain effect (MDE), a hypothesis which predicts the existence of a limit or geometric restriction in the distribution of the species, was evaluated. The fact of being a null model, directly excludes any environmental and evolutionary influence over the richness of species (Colwell, Rahbek & Gotelli, 2004). To evaluate the MDE, the specific richness data for each latitudinal band, was compared with a null model, through Monte Carlo simulations of the species richness curves. The simulated curves were based on empirical range sizes within a limit domain using the stochastic analytical model by Colwell & Hurt (1994) and Colwell (1999). In order to do this, we used 50,000 simulations with and without replacement to calculate the amplitude at a 95% confidence interval of the simulated curves (McCain, 2003; McCain, 2004). Sampling without replacement corresponds to a randomization technique, while sampling with replacement corresponds to a bootstrap method (Manly, 1997; McCain, 2004). The analysis was made using the Mid-Domain Null software (McCain, 2004).

Historical dynamic of the richness distribution

To evaluate the sink-source hypothesis, the degree of nestedness of a presence-absence matrix was estimated, where the columns represent the cells or grids and the rows, the species present in them. After this, it was proceeded to order the matrix according to the total sum by rows and columns, where the common species were placed on the top rows and the cells richest in species were placed on the columns on the left (see Valencia-Pacheco et al., 2011; Vallejos-Garrido et al., 2017). The nestedness metric based on overlap and decreasing fill index (NODF) was used (Almeida-Neto et al., 2008). The values varied from 0 to 100, where the highest values indicate an increase in the degree of nesting (Almeida-Neto et al., 2008; Valencia-Pacheco et al., 2011). The significance was calculated through null models using the Monte Carlo algorithm, contrasting values observed with a random probability distribution. A null model was used with fixed rows and equiprobable columns, where the totals observed by rows are kept, but the totals of the columns varied randomly (Patterson & Atmar, 1986; Gotelli, 2000). This null model conserves the species occurrence frequency, allowing that the specific richness varies equiprobably between the latitudinal bands (Valencia-Pacheco et al., 2011; Vallejos-Garrido et al., 2017). This model was chosen, considering a priori, that all the latitudinal bands could be occupied by the species. A total of 50,000 iterations were made to generate the frequency distribution. The nesting analysis was done in the NODF software (Almeida-Neto & Ulrich, 2011). This methodological approximation corresponds to a type of indirect test on the biogeographical hypothesis evaluation in a historical-geographical context, where the metacommunity dynamic is analyzed in a hierarchized or nested spatial system, implicitly incorporating a temporary component on considering large spatial scales (Holling, 1992; Ulrich, Almeida-Neto & Gotelli, 2009). As a result, a nested system will provide, non-random distribution patterns (Connor & Simberloff, 1979; Connor & Simberloff, 1986; Jackson, Somers & Harvey, 1992), where species that are rare in terms of occurrence, will only be present in the groups with greater richness, whereas the most common ones will be present in all the sites (Ulrich, Almeida-Neto & Gotelli, 2009).

Pattern of co-occurrence of species and regression analysis

To evaluate if the assemblage patterns of species co-occur or are randomly associated, we use the “C score” index (Stone & Roberts, 1990). We calculated a Checkerboard score (“C-score”), which is a quantitative index of occurrence that measures the extent to which species co-occur less frequently than expected by chance (Gotelli, 2000). A community is structured by competition when the C-score is significantly larger than expected by chance (Gotelli, 2000). After we compared co-occurrence patterns with null expectations via simulation (Gotelli, Hart & Ellison, 2015) using the SIM 9 algorithm, which considers the rows and column sums of the matrix are preserved. This means that each random community contains the same number of species as the original community (fixed column), and each species appears with the same frequency as in the original community (fixed row). The analyses were carried out in the Ecosim R package (Gotelli, Hart & Ellison, 2015).

We perform regression models to evaluate climatic, hydrological and topographic predictors as environmental factors that modulate the species richness pattern. Models were generated using the average UV radiation, average annual temperature, temperature seasonality, annual rainfall and seasonality of precipitation as climatic variables, topographic heterogeneity (range, Shannon, average and altitudinal range) as topographic variables and surface of water bodies for each latitudinal band (in Km²), and density of water courses as hydrographic variables. The climatic information was obtained from Wordclim database (Hijmans et al., 2005). The solar radiation was obtained from the glUV database (Beckmann et al., 2014), available at http://www.ufz.de/gluv. The variable altitude was obtained from SRTM database 90 m Digital Elevation Database v4.1 (Jarvis et al., 2008), available at https://cgiarcsi.community/data/srtm-90m-digital-elevation-database-v4-1. The topographic heterogeneity was calculated as the altitude range for each 1 degree band. Shannon’s index for topographic heterogeneity is based on textural features of the EVI index (Enhanced Vegetation Index), and was obtained from the EarthEnv database (Tuanmu & Jetz, 2015), available at https://www.earthenv.org/. The descriptive statistics of climatic, hydrological and topographic variables are shown in Table 1).

To evaluate the relationship between species richness and different environmental predictors, we use generalized linear models (GLM) using a Poisson distribution and a “log” link. Since the determination coefficients (R²) are not provided by the GLM models, a pseudo-R² was calculated which was obtained through the “rcompanion” package (Mangiafico, 2017). The selection of models was made through the Akaike information criterion (AICc) (Burnham & Anderson, 2002). The analyses were carried out through the MuMin package (Barton, 2017).

All statistical analysis were carried out in R (R Core Team, 2017).

Faunistic similarity patterns and geographic ranges

According to the concentration of species by latitudinal band, a greater specific richness is seen between 37° and 42°S, with a maximum richness at 41°S. It is also possible to highlight a lower specific richness at extreme latitudes (Fig. 1). The cluster analysis showed ten groups statistically significant, which correspond to: (1) 18–24°; (2) 25–30°; (3) 32–33°; (4) 34–36°; (5) 37–38°; (6) 41–44°; (7) 45–51°; (8) 46–47°; (9) 48–50°; (10) 53–54° (Fig. 2).

The geographic ranges of freshwater molluscs show a greater frequency of reduced geographical ranges (<2°). In general, the size distribution showed asymmetric bias towards the right (see Supplementary Material 1).

Effect and mid-domain effect

The species richness curve for each latitudinal band shows a marked MDE. The species diversity pattern observed, compared with the curves simulation at 95% for a stochastic null model without replacement, showed that 98% of the points of diversity are located within the signaled range, a result that is similar to what was obtained for the stochastic null model with replacement. Only one point of diversity does not fit the MDE model, this includes the band located at 40°–41°S (Fig. 3).

Rapoport’s rule

The regression analysis between the latitudinal range and the midpoint range showed a non-significant positive slope (p = 0.128) and a low value of the correlation coefficient (r = 0.014) (Fig. 4). These results indicate that the distribution of species richness is not fitted to Rapoport’s latitudinal rule.

Nestedness analysis of the richness distribution pattern

Results indicate a significantly nested matrix (except for the BR index), with values lower than expected by chance for the T index, and higher than expected by chance for the NODF index (Table 2). In addition, a significant nestedness was observed between rows (i.e., incidence of species) and columns (i.e., species composition) (Table 2). These results indicate that there is a significant nested pattern, which suggests a sink-source dynamic among the communities.

Patterns of species co-occurrence

The simulations of the null co-occurrence model indicated aggregate patterns, indicating that the species tend to be grouped greater than expected by chance (Table 3), indicating the existence of possible ecological mechanisms that underlie this pattern. The histogram of the simulated values and observed value is shown in Supplementary Material 3.

Climatic, hydrographic and topographic factors that explain the mollusc richness patterns

The environmental factors that explain the mollusc richness patterns indicate that variables related to environmental energy, such as seasonality temperature, showed a significant relation with richness, explaining 35% of the variability. Another important variable was the average altitude with 28% of the variance explained (Table 4). The other environmental, topographic and hydrographic predictors showed a low percentage of variance explained (Table 4). The full model, that is, with significant climatological, hydrographic and topographical predictors, explains 40% of the variability of specific richness. In this model, the variables related to topography and water availability are the most important according to the value of the estimated coefficient (Table 5).

Causal mechanisms of latitudinal richness patterns

Rapoport’s rule is not supported by our results because there was no increase in diversity towards the poles, as has been recorded in marine molluscs (Valdovinos, Navarrete & Marquet, 2003). However, there would be a significant MDE as a result of the existence of hard limits at both ends of the continent. In this sense, the geometry of the continent can impose powerful forces which limits the size of the geographic range, being able to counter the ecological and evolutionary processes used to explain Rapoport’s rule effect (Gaston, Blackburn & Spicer, 1998). While there is an increase in diversity towards higher latitudes, pattern reported also for other taxa (Hernández, Moreno & Rozbaczylo, 2005), only one point is greater than expected by chance at −52°S, so in that geographical area there would be a factor or mechanism different from the geometric restriction that would be modulating the richness.

The barriers that limit the dispersion of the organism, have also been preliminary explored in marine taxa (Hernández, Moreno & Rozbaczylo, 2005; Moreno et al., 2008), but not in freshwater organisms. In the MDE case, the first limiting barrier of the spread of freshwater molluscs is located at 18°S, caused by the arid diagonal, causing low water availability basins, most sporadically and under climatic conditions that are highlighted by the extreme aridity caused by the interaction of factors, such as the Pacific Anticyclone, Humboldt Current and elevation of the Andes Mountains (Villagrán & Hinojosa, 1997; Villagrán et al., 1994). Between 41°S and 56°S the continental fragmentation, determining the formation of fjords and glaciers (Camus, 2001), stands out as a barrier to the spreading, where although there is higher water availability, extreme weather conditions for the freshwater fauna stand out. This includes basins with less stability with predominating estuary conditions and the entry of sub-Antarctic waters (Silva, Calvete & Sievers, 1997). Studies made on zooplanktonic fauna from fjords show that the high environmental heterogeneity, limits the colonization of species, intensifying the isolation (Palma, 2006).

Source–sink dynamics and patterns of co-occurrence

For freshwater molluscs, there is a metacommunity sink-source dynamics, where the zones of greatest specific richness would act as a dispersion source towards other latitudes. According to the results obtained, the latitudinal areas comprising 38–40°S (Araucano Lakes region), would be a source of species, which may colonize extreme latitudes (i.e., sinks). The dynamics of nestedness can be explained by simultaneous balance processes between extinction (Patterson & Atmar, 1986; Bruun & Moen, 2003; Whethered & Lawes, 2005) and selective colonization events (Ulrich, Almeida-Neto & Gotelli, 2009; Valencia-Pacheco et al., 2011), especially in fragmented habitats, current situation of our hydrological resources or even their tolerance, depending on environmental filters (Greve et al., 2005; Driscoll, 2008), as a direct consequence of the geomorphological, climatic and historical events that have determined the current condition of the drainage basins (Fuentealba, Morrone & Figueroa, 2010). Regarding the distributional patterns, species that are rare in terms of incidence are only present in the groups with greatest richness, whereas the most common ones are present throughout the entire latitudinal extension (Ulrich, Almeida-Neto & Gotelli, 2009), represented in our study by the geographic extremes of our study area, previously characterized by extreme conditions and that depend of the tolerance-adaptation of the colonizing species. Another important finding was that the patterns of co-occurrence of freshwater molluscs are not random, or in other words, an aggregate pattern, namely, the species tend to co-occur in a greater way than expected by chance, so potential processes that cause this co-occurrence might be due to the similarity in habitat requirements (Gotelli & McCabe, 2002) or existence of positive biotic interactions (Ovaskainen, Hottola & Siitonen, 2010). This aggregation of species can be linked to a large extent to historical dynamics, such as Pleistocene glaciations and postglacial tectonic events that affected the hydrography of Chile (Fuentealba, Morrone & Figueroa, 2010).

Otherwise, contemporary mechanisms recognize temperature as the main predictor of the distribution of the richness of freshwater fish (Habit et al., 2012), corroborating that this climate variable is a powerful descriptor of freshwater fauna, and for molluscs is able to explain by itself 35% of the variability. On the other hand, our regression models with multiple predictors, indicated that the variables of habitat heterogeneity, topography and water availability (see Table 5) revealed a greater variability explained (40%); therefore, the configuration of the landscape (i.e., heterogeneity) (Stein, Gerstner & Kreft, 2014), as well as the availability of habitat and water regimes are efficient modulators of the richness of freshwater molluscs.