Integrating alpha, beta, and phylogenetic diversity to understand anuran fauna along environmental gradients of tropical forests in western Ecuador

Abstract The study of current distribution patterns of amphibian species in South America is of particular interest in areas such as evolutionary ecology and conservation biology. These patterns could be playing an important role in biological interactions, population size, and connectivity, and potential extinction risk in amphibians. Here, we tested the effects of spatial and environmental factors on the variation, turnover, and phylogenetic diversity of anuran amphibian species in tropical forests of western Ecuador. Data for presence/absence of 101 species of 34 genera and 10 families registered in 12 sites (nested in four biogeographic units) were obtained through fieldwork, museum collections, and literature records. We examined the influence of geographical, altitudinal, temperature, and precipitation distances on differences in anuran composition between sites. We found significant positive correlations among all of these variables with anuran distribution. The greatest alpha diversity (species richness) was found in the Equatorial Chocó biogeographic unit. Equatorial Pacific biogeographic unit could act as a transition zone between the Equatorial Chocó and Equatorial Tumbes. The western Andes (Western Cordillera biogeographic unit) was the most dissimilar and exhibited a higher species turnover rate than the other biogeographic units. Our results suggest that precipitation and elevation play a key role in maintaining the diversity of amphibian species in western Ecuador.

Studies from a global perspective have been carried out to analyze how the richness and turnover of amphibian species respond to different environmental and spatial gradients (Buckley & Jetz, 2007, 2008 and also to test the influence of phylogenetic history on the global patterns of amphibian species richness (Fritz & Rahbek, 2012).
Overall, studies of the spatial patterns of species assemblages are urgently required to delineate conservation strategies in ecosystems under strong anthropogenic pressures such as the evergreen forests of Chocó and Equatorial dry forests, two of the most fragmented forests in western Ecuador (Dodson & Gentry, 1991;Escribano-Avila et al., 2017).
Here, we analyzed the community composition and phylogenetic structure of anurans occurring in 12 sites of four biogeographic units of western Ecuador. For this purpose, we use the alpha, beta, and phylogenetic diversity in order to establish questions about the factors that determine the variation of the diversity of anurans between biogeographic units of western Ecuador. We hypothesized that diversity (alpha, beta, and phylogenetic diversity) would depend on local-level composition of species in the sites and their location along environmental gradients. In summary, this study provides a baseline of the phylogenetic diversity of anuran species in western Ecuador, allowing us to propose "hot spots" of amphibian diversity in this region.  Figure 1).
Species identification was performed using taxonomic keys and specialized literature, including the original descriptions of the species recorded (e.g., Coloma, 1995;Lynch & Duellman, 1997).

| Anuran assemblage alpha diversity and variation in species composition
The species richness (SR) was calculated for each site, the forest type and BU. Comparisons of S on BU's and forest type, respectively, were analyzed with linear regression models and perform ANOVA on the data; afterward, a Tukey test was used to determine which relationships were statistically significant.
These analyses were performed in R 3.3.2 software (R Core Team, 2016). In order to address which species are shared and which are distinct in the anuran assemblage, we calculated the Jaccard index for pairs of sites. To represent the ordering relationships among sites per BU in a reduced and predetermined number of axes, an ordination analysis (nonmetric multidimensional scaling analysis, 2D-NMDS) was performed on matrices constructed from Jaccard indices. To test for differences in species composition dissimilarity, a permutational multivariate analysis of variance (PERMANOVA, Anderson, 2001) was performed on the Jaccard similarity matrices using "Biogeographic Unit" as a fixed factor.
The probability value (p perm ) was calculated from a pseudo-F distribution with 10,000 permutations. All analyses were performed using PERMANOVA+ in the PRIMER v6 statistical package (Clarke & Gorley, 2006). In order to evaluate the effects of geographical distance on dissimilarity in species composition (distance-decay relationship) and to have another measure if species turnover and beta diversity, we calculated the distance in km between all pairs of sites and plotted the calculated Jaccard index. The distance-decay relationship was quantified, in the data set the linear relation of Jaccard similarity to geographic distance (on both logtransformed and original scales) was assessed using linear regression. This analysis was made using the Vegan package (Oksanen et al., 2016) in R 3.3.2 (R Core Team, 2016).

| Phylogenetic diversity and phylogenetic structure
We used a phylogeny of anuran species present in the 12 sites, since some species did not have available sequences and other species have only been identified to the genus level, as is the case of several reported Pristimantis, we downloaded 70 sequences (761 base pairs in each sequence) of 16S mitochondrial gene available from Genbank (see Appendix S1). Phylogenetic relationships performed with 1,000 ultrafast bootstrap replicates and the most appropriate substitution model based on the Bayesian information criterion (BIC) were inferred using IQ-tree (Nguyen, Schmidt, von Haeseler, & Minh, 2015) and ModelFinder (Kalyaanamoorthy, Minh, Wong, Haeseler, & Jermiin, 2017), respectively. The sequences were analyzed under the TIM2 + I + G model, and the maximum likelihood tree was saved as Newick format for analysis. With this phylogeny and the community presence/absence matrix as input, we performed phylogenetic metrics for each site. We calculated two diversity measures, first the phylogenetic diversity (PD) index, defined as the sum of branch lengths between root and tips for a community (Faith, 1992) for each site, and we compared PD with forest type and BU's; first, we fit a linear regression models, and then, we perform an ANOVA on the data; afterward, a Tukey test was used to determine which two variables had significant differences. These analyses were performed in R 3.3.2 software (R Core Team, 2016). And then, we calculated the standardized effect size of Faith's PD (SES PD ) for all the sites. In order to assess how phylogenetically related are the average pair of species in a site, we use two indices proposed by Webb, Ackerly, McPeek, and Donoghue (2002)  Relatedness Index (NRI) and Nearest Taxon Index (NTI), respectively (Pearse, Jones, & Purvis, 2013). The community phylogenetic structure was calculated as follows: MPD calculate the mean pairwise distance between all species in each site. On the other hand, MNTD calculates the mean nearest taxon distance, the mean of the branch lengths connecting each species to its closest relative (Webb, 2000). We use a null model of randomly shuffling tip labels across the tips of the phylogeny with 1,000 Map of Ecuador showing the study sites. 1 = Río Canandé, 2 = Bilsa-Mache Chindul, 3 = Jama-Coaque, 4 = Ayampe-Machalilla, 5 = Chongón-Colonche, 6 = Churute, 7 = Buenaventura, 8 = Achiotes-El Faique, 9 = Cordillera Arañitas-La Ceiba, 10 = Río Guajalito, 11 = Río Faisanes, 12 = Quebrada Zapadores. Approximate distribution of terrestrial South American ecoregions modified from Olson et al. (2010) runs for each analysis (site). The reported p-value was calculated with a two-tailed test, thus, significance at the threshold α = 0.05 level is achieved when p ≤ .025 or p ≥ .975 (Cadotte & Davies, 2016).
Positive SES values and high p-values (p ≥ .975) indicate phylogenetic evenness and greater phylogenetic distance among co-occurring species than expected, and negative SES values and low p-values (p ≤ .025) indicate phylogenetic clustering and small phylogenetic distances among co-occurring species than expected (Kembel et al., 2010). The analyses were performed with PICANTE package (Kembel et al., 2010) in R 3.3.2 software (R Core Team, 2016).

| Effect of environment on anuran diversity
We evaluated the correlation between abiotic and biotic components. Correlation tests were performed between dissimilarity matrices (Bray-Curtis dissimilarity) of environmental variables (precipitation, temperature, and elevation between sites), and the inverse value of Jaccard similarity (Jdissim). These models were calculated using the package Vegan in R (Oksanen et al., 2016) and following the recommendation by Legendre, Borcard, and Peres-Neto (2005), Legendre, Fortin, and Borcard (2015) and Legendre and Legendre (1998). To detect multicollinearity of predictor variables, we used a statistic called the variance inflation factor (VIF) (Fox & Monette, 1992). The square root of the VIF indicates the degree to which the standard error is, comparing if a predictor variable was not correlated with other predictor variables in a model. As a gen-  (Table 1).
Also, we evaluated if there was an effect of the three environmental variables previously mentioned (see Table 1) on S and PD (response variables), in order to select the best-fit model that explains the maximum amount of variance, we created two multiple regression models to explain the two response variables, this was done in R 3.3.2 (R Core Team, 2016). To look for evidence of nonlinearity in the relationship between the dependent variable (S and PD) and the independent variables (precipitation, temperature, and elevation), we used component plus residual plots whit the crPlots function in the R package car (Fox & Weisberg, 2011). Similarly, we performed a global validation of linear model assumptions with the R package gvlma (Pena & Slate, 2006). To assess the relative importance, that is, the contribution of each of the predictor variables on the response variable in a multiple regression model; we use the R package relaimpo which provides several reasonable metrics, such as lmg that propose averaging sequential sums of squares over all orderings of regressors (Lindeman, Merenda, & Gold, 1980) for assessing relative importance (percent contribution) of each correlated predictor (regressor) in a linear model (Grömping, 2006).

| Anuran assemblage alpha diversity and variation in species composition
A total of 101 species of frogs were recorded (Table S1)  The sites with lowest S were in the Equatorial Tumbes BU; these sites included Cordillera Arañitas-La Ceiba (site 9) and Achiotes-El    (Table S2).
In relation to species composition variation in biogeographic units, the nMDS indicated that the sites form two main clusters of low similarity. The most distinct group included the Western Cordillera sites while the second group included the Chocó, Pacific, and Equatorial Tumbes sites ( Figure S3). This result was complemented by high goodness of fit resulting from repeated optimization; F I G U R E 2 Boxplot showing species richness (SR) (a and b) and phylogenetic diversity (PD) (c and d) per biogeographic unit and forest type. Outliers are shown with black circles. Thick horizontal black lines indicate means the stress function of the nMDS was 0.069, which indicates that the scaling was properly adjusted. From a total of 44 species recorded in the Equatorial Chocó, only eight species were also present in W.
Cordillera; similar variation occurs with Equatorial Pacific (43 species), and only eight were also in W. Cordillera. BU's that shared more species were E. Chocó-E. Pacific (18 species shared) and E.
The results of the PERMANOVA analyses showed significant  (Table S5).
When we evaluate the geographical distance with the Jaccard dissimilarity in species composition, we found a statistically significant positive correlation coefficient (r = .461, p < .001), indicating that there is a distance decay of similarity (communities far away from each other have more different species compositions) (Figure 3).

| Phylogenetic diversity and phylogenetic structure
The relationship between SR and PD for the community data showed that PD is strongly correlated with SR (p < .001, R 2 = 0.86) ( Figure S1). As expected, PD was found to be the highest in moist forests of Equatorial Chocó and the lowest in dry forest of Equatorial Tumbes. The highest SES PD was found in Buenaventura (Equatorial Tumbes, transition forest), while the lowest SES PD in Quebrada Zapadores (Western Cordillera, montane forest) ( Table 2). There were differences in PD among different forest types (F (3,8) = 7.12, p = .012) and among different BU's as predictor (F (3,8) = 6.06, p = .019) (Figure 2 (Table S2).
Phylogenetic structure of the anuran communities varied across the spatial extent of the study area (Table 2). There were no significant differences of SES MPD when this standardized effect size was calculated for the different sites, forest types, and biogeographic units. Similar results were found with SES MNTD , there were no significant differences in the same three previous levels, only the site Q. Zapadores had a p-value = .015 calculated with a two-tailed test ( Figure S2, Table 2).

| Effect of environment on anuran diversity
We do not find multicollinearity evidence of variable predictors in this model, elevation  (Table S4). However, a significant relationship was found between S and precipitation alone (p = .031), yet there was no significant effect of elevation and temperature on S (p > .05). PD was not significantly correlated with each of the three environmental variables examined (F (3,8) = 3.38, p = .075), as in S precipitation was significant (p = .046). However, the p-value for elevation and temperature (0.324 and 0.885, respectively) is greater than the common alpha level of 0.05, which indicates that were not statistically significant (Table S4).
According the relative importance of three environment variables on PD and S, we use the method called lmg (Grömping, 2006).
Precipitation had the highest relative importance or regressor contribution (R 2 ) on S (lmg = 79.8%) and PD (lmg = 72.5%) (Table 3). When comparing lmg with other methods to measure relative importance, similar results were obtained ( Figure S4).

| D ISCUSS I ON
We found heterogeneity in the alpha, beta, and phylogenetic diversity among the four Ecuadorian biogeographic units. The Equatorial Chocó was the unit with the highest species richness, which can be mainly explained by climatic factors such as high average annual F I G U R E 3 Distance-decay plot, the distance was measured in km between all pairs of sites (Spatial Euclidean distance) and plotted the calculated Jaccard index (Species composition dissimilarity) rainfall (2,000 mm; Sierra, Cerón, Palacios, & Valencia, 1999) and dry months with rainfall <100 mm per month (Gentry, 1995;Mooney, Bullock, & Medina, 1995;Pennington, Lavin, & Oliveira-Filho, 2009); likely as a result, amphibian richness was lower in Equatorial Tumbes.
As in species richness, the differences in phylogenetic diversity in community assemblage of anurans are related to differences in precipitation; this is relevant to understand the turnover across different sites, forest types, and biogeographic units. The effect of precipitation may be due to the fact that alpha and phylogenetic diversity are inherently positively correlated, since a greater number of species almost always correlates with a greater genetic divergence summarized in a phylogeny (Cadotte & Davies, 2016;Venail et al., 2015).
The similarity/dissimilarity among the anuran communities studied here can be explained in part by environmental or climatic factors (Lynch & Suárez-Mayorga, 2002). Here, we used ordination analysis to identify significant relationships between biological (e.g., number of species, species turnover and phylogenetic diversity) and environmental variables (e.g., temperature, precipitation) among sites. We found that biogeographic dissimilarity measured as spe- In this study, the genus Pristimantis (Craugastoridae) had the highest number of recorded species (36.6% of the total species recorded). Lynch and Duellman (1997) show that Pristimantis species from the lowlands of western Ecuador have wider distribution ranges than congenerics from the Andes; this could explain the high number of species of this genus recorded in the Western Cordillera (Pristimantis richness was much higher than that of other genera in the same unit).
According to the phylogenetic structure in the assembly of the communities, most of the communities are phylogenetically grouped (e.g., La Ceiba-Cordillera Arañitas, Quebrada Zapadores, Bilsa-Mache Chindul, Jama-Coaque), however, no significant differences were found. Only two sites (Buenaventura and Río Faisanes) presented are high phylogenetically overdispersal or a greater phylogenetic distance between coexisting species than expected. These results could be explained given that the community assemblages consist mainly of species that have diverged relatively recently.
Because species richness and distribution patterns at local scales are the result of complex biotic and abiotic interactions at many spatial and temporal scales (Wisz et al., 2013), there is no single cause of these patterns. On the other hand, environmental factors such as precipitation or elevation can influence ecological processes in organisms, and therefore their capacity for dispersion and persistence in different environments (Brown & Lomolino, 1998). By analyzing precipitation, temperature, and elevation, we have sought to widen our inference of the factors affecting the distribution of amphibians in western Ecuador. Others have found that the species richness of amphibians is influenced by factors such as temperature, geography, and precipitation (e.g., Ortiz-Yusty et al., 2013;Soares & Brito, 2007).
Nonetheless, these factors are not the only studied, other studies have also found that anuran diversity has been determined as a response to either different types of vegetation, distance to water bodies or environmental heterogeneity (e.g., Goncalves, Crivellari, & Conte, 2015;Ribeiro et al., 2017). Here, from the regression analysis, we also found that precipitation could have a strong effect on the diversity of amphibians (SR-PD; Figure 5, Figure S4).
The results of this work may suggest that Equatorial Pacific would act as a transition zone between Equatorial Chocó (wet/ moist northern forests) and Equatorial Tumbes (dry southern forests) (see Figure 1) in terms of anuran species composition. Overall, our results support this suggestion, as has been previously defined for the area (Valverde, 1991;. These three biogeographic units mentioned above, could be characterized by high species turnover, which would follow a latitudinal gradient. For example, among these three units the species composition of some dendrobatids varies latitudinal and ecologically (see Coloma, 1995;Grant et al., 2006Grant et al., , 2017Santos et al., 2009;Tarvin, Powell, Santos, Ron, & Cannatella, 2017, for distribution data). As in the case of replacement, in the biogeographic We conclude that environmental factors such as precipitation, elevation, and temperature could affect the diversity of anurans in Western Ecuador. For example, the composition of anuro-fauna in the forests of the Western Cordillera, sites that present low temperatures on average, is markedly different from the composition found in three other biogeographic units, presenting a high species richness but belonging to a few taxonomic groups (e.g., rainfrogs of genus Pristimantis or glassfrogs of family Centrolenidae). On the other hand, the high rainfall in the Ecuadorian Chocó lead to the sites in this biogeographic unit maintain a constant humidity throughout the year, which make available numerous ideal microhabitats for the persistence of several amphibian species distributed in different clades within a phylogenetic tree.
Finally, these ecosystems in the coast and western Andes of Ecuador have already been categorized as high priority areas for conservation and as high exposure risk zones (Cuesta et al., 2017;Sierra, Campos, & Chamberlin, 2002). Furthermore, given the high phylogenetic diversity of amphibians and even because most of the sites in this study do not have a formal declaration of forest protection, which could allow these forests to be considered as conservation areas of biodiversity (see also Arteaga et al., 2013;Cuesta et al., 2017;Lessmann et al., 2014).

ACK N OWLED G M ENTS
We would like to thank Diogo B. Provete, Guillermo D'Elía, Emily Marcial Quiroga for their help with some files.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
LA and JMG conceived the study concepts; LA, JMG, and MS-G designed the methodology; LA collected and analyzed the data; LA, MS-G, and JMG led the writing of the manuscript. All authors contributed critically to the drafts and provided final approval for publication.