Spatial relationships between fishes and amphibians: implications for conservation planning in a Neotropical Hotspot

Species distribution patterns are widely used to guide conservation planning and are a central issue in ecology. The usefulness of spatial correlation analysis has been highlighted in several ecological applications so far. However, spatial assumptions in ecology are highly scale-dependent, in which geographical relationships between species diversity and distributions can have different conservation concerns. Here, an integrative landscape planning was designed to show the spatial distribution patterns of taxonomic and functional diversity of amphibians and fishes, from multiple species traits regarding morphology, life history, and behavior. We used spatial, morphological, and ecological data of amphibians and fishes to calculate the functional diversity and the spatial correlation of species. Mapping results show that the higher taxonomic and functional diversity of fishes is concentrated in the West Atlantic Forest. Considering amphibians, are located in the East portion of the biome. The spatial correlation of species indicates the regions of the Serra do Mar and the extreme southern part of the Central Corridor as the main overlapped species distribution areas between both groups. New key conservation sites were reported within the Brazilian Atlantic Forest hotspot, revealing cross-taxon mismatches between terrestrial and freshwater ecosystems. This study offers useful spatial information integrating suitable habitats of fishes and amphibians to complement existing and future research based on terrestrial and freshwater conservation. New priorities for biodiversity conservation in rich-species regions highlight the importance of spatial pattern analysis to support land-use planning in a macroecological context.


Introduction
In the face of the global biodiversity crisis, multi-species conservation strategies are essential to achieve ecological outcomes (Mendenhall et al. 2012;Dinerstein et al. 2020;Jung et al. 2021). The use of spatial patterns of biodiversity components to support conservation actions has become one of the main research topics in ecology (Rodrigues et al. 2011;Westgate et al. 2014). There is an urgent need for reliable data on species composition and distribution for conserving biodiversity habitats and reducing biodiversity losses (Dinerstein et al. 2020;Mokany et al. 2020). However, time and money are running out for biodiversity strategies, which hampers the development of cost-effective conservation planning (O'Bryan et al. 2020;Strassburg et al. 2020;Yang et al. 2020). It is also necessary to consider that measuring the world's huge biodiversity patterns is difficult under the intense human pressure across species' ranges in natural ecosystems (Williams et al. 2005;O'Bryan et al. 2020). As an alternative, many conservation plans rely on cross-taxon congruencies to safeguard biodiversity patterns in balance with environmental emergency management Westgate et al. 2014;Yang et al. 2020).
There are several elaborated definitions of cross-taxon congruence for conservation planning (Westgate et al. 2014). According to Heino (2010), it can be defined as the strong correlation of biodiversity measurements between taxonomic groups across a set of localities. Usually, assemblages of well-known and mapped species are used to predict less-studied groups' presence (Lawler et al. 2015). Thus, cross-taxon congruence aims to accelerate the process of delimiting areas for conservation based on some species' requirements (Fraveau et al. 2006;Weerd and Haes 2010;Rodrigues et al. 2011). It is known that only cross-taxon congruence is not enough for one species group to be a good predictor for the distribution of another (Gaston 1996;Padial et al. 2012;Duflot et al. 2022). However, this relationship is essential, as it allows using a better-known taxon and managing another lesser-known taxonomic group. Lawler et al. 2015). Although this approach may introduce some practical limitations, it has been shown useful to investigate ecological and spatial relationships in aquatic and terrestrial ecosystems (Heino et al. 2005;Heino 2010;Westgate et al. 2014).
Amphibians and continental fishes are both dependent on humid environments (i.e., freshwater environments and humid forests) and have important roles in aquatic and terrestrial ecosystems (Holmlund and Hammer 1999;Wells 2007;Haddad et al. 2013;Villeger et al. 2017). The importance of these taxa in the ecosystem can be measured by their functional traits. Functional traits are the characteristics of physiological, morphological, or life history of species that are measured through functional diversity in an area (Petchey and Gaston 2006). The Brazilian Atlantic Forest represents a biodiversity hotspot with high species richness, high levels of endemism, and high rates of habitat loss due to human occupation (Myers et al. 2000;Oyakawa et al. 2006;Ribeiro et al. 2009). Fragmentation of the Atlantic Forest has promoted the success of generalists and the loss of specialist species leading to taxonomic and functional impoverishment of species (Sfair et al. 2016;Lourenço-de-Moraes et al. 2020). Within this heterogeneous ecosystem, some areas stand out as extremely important hotpoints (i.e. specific points of high species richness) due to the high potential for biodiversity corridors (Carnaval et al. 2009;Lourenço-de-Moraes et al. 2019a;Campos et al. 2020).
In this paper, the spatial patterns of taxonomic and functional diversity of fishes and amphibians were assessed to support efficient landscape planning for aquatic and terrestrial ecosystems in the Brazilian Atlantic Forest. Crosstaxon congruence of taxonomic and functional diversity patterns among fishes and amphibians was tested as a shortcut for mapping biodiversity components from small spatial scales. The results may help apply conservation strategies considering each group's functional role (i.e., fishes and amphibians) in an integrative way. This study offered an essential contribution to the conservation of terrestrial and freshwater ecosystems, having crucial implications for ecosystems under human pressure on different spatial scales.

Study area
Research design and analysis were conducted in the most threatened biodiversity hotspot on Earth (Myers et al. 2000)--the Brazilian Atlantic Forest. This biome had an original coverage of about 1.5 million square kilometers extended in the tropical and subtropical regions (Myers et al. 2000), of which only 12.9% remain in Brazil, Argentina, and Paraguay (Ribeiro et al. 2009). In the Brazilian territory, the biome comprises the coastal region from the Rio Grande do Sul to the Rio Grande do Norte states ( Fig. 1), retaining only around 11% of its original extension Spatial range, species richness, and functional traits Geographical distribution records of fishes and amphibian species were obtained according to three databases: (1) Global Biodiversity Information Facility--GBIF (https://www.gbif.org); (2) Species Link (http://www. splink.org.br) and (3) IUCN Red List database, version 2021.1 (IUCN 2021). To prepare the spatial data of the species, the ArcGIS Pro software (ESRI 2019) was used to build a presence/absence binary matrices for combining the species distribution databases. Therefore, the studied area was mapped in a grid system of 0.1 latitude/longitude degrees, creating a network with 10,359 grid cells. The Atlantic Forest area limits considered in this study followed the definitions of the Brazilian Ministry of Environment (MMA 2018).
For analyses, only species data regarding adult individuals of fishes and amphibians within the Brazilian Atlantic Forest were considered. According to specialized literature for each taxonomic group, all species distributions of fishes and amphibians were compiled in a dataset (see Supporting Information). Species nomenclature was standardized following the Catalog of fishes database (Fricke et al. 2022) and the Amphibians Species of the World database (Frost 2022). We did not include in the analysis species that have not yet been formally described and/or for which there is no accurate taxonomic identification.

Ecological niche modeling
Spatial occurrences of species were transformed in presence/absence matrices, using the "Spatial Join" ArcGIS toolbox in ArcGIS Pro software (ESRI 2019). Then, vector files based on expert knowledge of the species' ranges and forest remnant polygons were combined into an overall coverage for species distribution modeling. Only spatial occurrences of species distribution data intersected on at least one grid cell were considered (i.e.,~10 km 2 ). In addition, ecological niche models (ENMs) were used to predict and map the potential distribution area of each fish and amphibian species in the Atlantic Forest.
For ecological niche models (ENMs), we apply the species occurrence matrix and the layers of climaticenvironmental variables, resulting in a suitability matrix. We used the following bioclimatic variables in the modeling process: (1) annual mean temperature; (2) annual temperature range; (3) precipitation of the wettest month; (4) precipitation of the driest month; and (5) precipitation of the warmest quarter, which were determined as the most independent by a factorial analysis. We obtained these variables from CMIP6-Coupled Models Intercomparison Project Phase 6 (https://www.worldclim.org/data/cmip6/cmip6clim ate.html) and downscaled them to the resolution of 0.1 degrees. We used simulations provided by four Atmosphere-Ocean General Circulation Models (AOGCMs): CCSM (Community Climate System Model), CNRM (Centre National de Recherches Météorologiques), MIROC (Model for Interdisciplinary Research on Climate), and MRI (Meteorological Research Institute) for the consensus model. Original data resolution varied from 0.9°t o 2.8°(in longitude and latitude) and the climate variables were re-scaled to fit our grid resolution. We also used altitude and the number of rivers as a predictor of richness and dispersion for fishes and amphibians from the dataset available at WorldClim Global Climate Data (http://www.w orldclim.org).
We performed four conceptually and statistically different ENMs based on presence data (i.e., only occurrences are known, absences are unknown) using the algorithms: (1) Bioclim (BIO, Busby 1991) based on bioclimatic envelope logic; (2) Gower Distance and Euclidean Distance (GD, EUD, Carpenter et al. 1993) based on environmental distance approach; (3) Maximum Entropy (ME, Phillips et al. 2006) and Random Forest (RF, Breiman 2001) based on machine learning technique; and (4) Ecological Niche Factor Analysis (ENFA, Hirzel et al. 2002) based on multivariate analysis, and Genetic Algorithm for Ruleset Production (GARP, Stockwell and Noble 1992). Given each model's particularities, they provided different predictions, generating uncertainties about which model is more appropriate to represent the geographical distribution of species (Diniz-Filho et al. 2009). To overcome this uncertainty and minimize errors, we employed the ensemble forecasting approach, which offers a consensus of multiple models (Araújo and New 2006). The main idea of ensemble forecasting is that different sources of errors will affect each niche model in different ways and, by obtaining a consensus result of these models, errors will tend to cancel each other out and produce a more trustworthy and conservative solution . Assuming that the species richness consensus model (CONS) reduces uncertainty and error associated with alternative ENMs, we interpreted only the CONS model's range sizes.
We randomly partitioned the presence and absence data of each species into 75% for calibration (or training) and 25% for evaluation (or test); repeating this process 10 times by cross-validation for all models. For each ENM, we converted the continuous predictions of suitability into a binary vector of 1/0 (presence/absence in each cell), finding the threshold that maximizes sensitivity and specificity values in the receiver operating characteristic (ROC). The ROC curve is generated by plotting the fraction of true positives versus the fraction of false positives at various threshold settings. The distribution areas were estimated by obtaining 280 predictions (7 models × 10 randomizations × 4 AOGCMs) for each species. This allowed us to generate a frequency of projections in the ensemble. Then, we generated the frequency of projections weighted by true skill statistics (TSS) (the best models according to this metric have more weight in our consensus projections). The TSS range from −1 to +1, where values equal to +1 is a perfect prediction and values equal to or less than zero is a prediction no better than random (Allouche et al. 2006;Eskildsen et al. 2013). We considered the species present only in cells where at least 50% of models retained in the ensemble point out the species as the present. In our analyses, we obtained the CONS for each AOGCM. Thus, we obtained the final maps of species richness through the average of values projected by CONS for each grid cellconsidering the different GCMs.
All models were performed using the computational platform Bioensembles (Diniz-Filho et al. 2009) and mapped results using the software SAM v.4.0 (Rangel et al. 2010). To determine the species richness patterns of fishes and amphibians of the Atlantic Forest, the modeling strategy at the community level of "predict first, assemble later" were employed (Overton et al. 2002), where the ranges of individual species are modeled one at a time as a function of environmental predictors and then overlapped to obtain the species richness. This approach is called stacked Species Distribution Models--SDMs (Dubuis et al. 2011;Hof et al. 2012;Mata et al. 2017).

Data analyses
After the establishment of Taxonomic Diversity (TD), the Functional Diversity (FD) from the generated functional traits database for fishes and amphibians were estimated according to the FD index (Petchey and Gaston 2006). The FD index measures the extent of complementarity between species characteristics values (Petchey and Gaston 2006), relating directly to the ecological niche concept (Cianciaruso et al. 2009). The protocol proposed by Petchey and Gaston (2006) was followed to calculate functional diversity (FD): (1) construction of a species-trait matrix; (2) conversion of the species-trait matrix into a distance matrix; (3) clustering distance matrix into a dendrogram (UPGMA); and (4) calculating functional diversity by summing dendrogram branch lengths of species community. Besides, the FD considers the dependence of distances between species in the n-dimensional space (Petchey and Gaston 2006). The method proposed by Pavoine et al. (2009) was used to create the distance matrices through the distance of Gower (1971). The Gower distance was selected due to being the most recommended distance for simultaneously assessing continuous and categorical variables (Gower 1971).
Independent swap null models (Gotelli and Entsminger 2001) were used to verify whether FD was influenced by TD (Devictor et al. 2010), according to the protocol proposed by Swenson (2014). The null model is independent of the species richness of an assemblage. In this way, it ensures that the patterns of assembly of characteristics do not simply reflect particular species' differential occurrence (Swenson 2014). Thus, every~10 km² of the grid cells was tested to verify whether the values attributed to functional diversity were higher, equal, or lower than expected in a randomized and non-random distribution, assuming a random distribution in which every species could occupy any grid cell in the biome (for further details, see Swenson 2014). For analysis, were computed 1,000 replicates of FD, obtaining a p-value of predicted FD as compared to the distribution of the random replicates. All analyses were performed using the packages "ade4" (Dray and Dufour 2007), "picante" (Kembel et al. 2010), "FD" (Laliberté and Legendre 2010), and "vegan" (Oksanen et al. 2018) through software R (R Core Team 2020).
We then used simple linear regression models (testing normality through the Shapiro-Wilk test) to evaluate TD (CONS) and FD correlations in each grid cell. Correlation matrices were used to compare the topographic patterns and spatial references (altitude and latitude) to the values obtained by the richness consensus model for TD and FD in each grid cell. Thus, we correlated the values obtained for TD and FD with altitude and latitude using simple linear regression models. We also calculated the values overlap of TD and FD between groups in each cell. For this, we considered only areas with values ≥50% of the total value obtained for each index (i.e., TD and FD). From these areas containing the highest values, we build overlap maps between TD and FD. Thus, we avoided the possibility of obtaining statistical significance that could emerge from high congruence in the order of areas with low TD and FD values (Fattorini et al. 2012). We performed all analyses on the R program using the "vegan" package (R Core Team 2020).

Results
The results from TSS for most species of fishes had a mean and standard deviation of 0.65 ± 0.15 and amphibians 0.61 ± 0.11, both the results indicating a relatively high fit model. With the maximum mapped values of 295 fish species (TD) and 36.6 FD index per grid in the 10,359 cells evaluated. The highest indices of TD and FD of fish occurred in the western Atlantic Forest (Fig. 2a, c). This region corresponds to the areas that make up the Paraná River basin. We also identified high TD and FD values for fish species in the eastern Atlantic Forest, in the Serra do Mar region. This area is composed of a small part of the Paraná River basin, the southeastern and eastern Atlantic basins. In addition, the linear correlation analysis of the data showed a high correlation between the TD and FD index of fish species (R² = 0.98; P < 0.001) (Fig. 2b). This means that the regions that contain a high species richness also harbor a greater diversity of ecosystem functions.
For amphibians, the maximum values were 155 species (TD) and 17.8 FD index per grid in the 10,359 cells evaluated. The highest concentrations of these indices are in the eastern of the biome (Fig. 2d, f), mainly in the Serra do Mar, Central Corridor of Atlantic Forest, and Endemism Area of Pernambuco (altitude areas) in ascending order from south to north (For details about the regions see Fig. S1). Amphibian TD and FD values are also highly correlated (R² = 0.95; P < 0.001) (Fig. 2e).
In general, mapping the assemblages' relationships, important spatial mismatches, and low congruencies were revealed among the biodiversity components. The results showed no correlation between TD of fish and amphibians (R² = 0.006; P < 0.001), as well as for FD (R² = 0.003; P < 0.001 - Fig. 3). On the cross-taxon congruence of points ≥50% of the total value, fishes and amphibians also presented low spatial overlap in TD (0.19%) and FD (2.58%) (Fig. 4). Overlapping areas between TD and FD of fishes and amphibians correspond mainly to the regions of the Serra do Mar, the southern end portion of the Central Corridor of the Atlantic Forest, and the Pernambuco Endemism Centre. Further, when considered only the areas of TD and FD indexes ≥50% of the total value, it was observed a higher correlation for TD (R² = 0.69; P < 0.001) and for FD (R² = 0.95; P < 0.001).

Discussion
The cross-taxon congruence patterns highlight three main areas of interest ('hotpoints') aquatic and terrestrial conservation within the Brazilian Atlantic Forest (see Fig. 4). Besides, Serra do Mar and Central Corridor also refuge (Carnaval et al. 2009;Lourenço-de-Moraes et al 2019a, b) the only regions in the biome that displays high congruence between the taxonomic and functional diversity of fishes and amphibians. Considering only fishes, the most important hotpoint is located in the west of the biome, in the Paraná River basin--Upper Paraná River Floodplain. This is the last stretch of the Paraná River that is free of dams and received some species of fish after a set of waterfalls that acted as natural geographic barriers to be submerged by Itaipu Reservoir in 1983 (Agostinho and Júlio 2002;Agostinho et al. 2015). Some invasions of species have changed the trophic relationships (Benedito-Cecílio and Agostinho 1999; Alves et al. 2017) due to the recent impact of the Itaipu Reservoir. Thus, it is highly recommended the use of cross-taxon congruence between these groups to define conservation strategies and biomonitoring studies in the Brazilian Atlantic Forest.
Most of the Atlantic Forest ecosystems display a highly heterogeneous TD and FD distribution between fishes and amphibians. For fish, the main hotpoints are located in the west of the biome, along the Paraná River basin. This basin is among the world's five largest water systems (Rosa 1983;Agostinho et al. 1995). Several studies point to the importance of this region as a high biodiversity source (Bini et al. 2001;Agostinho et al. 2007;Graça and Pavanelli 2007;Lourenço-de-Moraes et al. 2021). In the last decades, the full extension of the basin in the Brazilian territory has been impacted by the construction of dams to produce hydroelectric power (Rosa 1983;Agostinho et al. 1995;Stevaux et al. 2009). Considering only the Upper Paraná River, for instance, there are at least 150 dams built along its course (Stevaux et al. 2009).
The impoundment of rivers creates a strong perturbation of fluvial dynamics (Stevaux et al. 2009). Among these disturbances is the alteration of flood magnitude (Souza Filho 2009), the turbidity reduction by retention of suspended load (Grimshaw and Lewin 1980), bank erosion (Rocha et al 2013), and reduction in bed form sizes (Dos Santos et al 2017). These alterations can induce significant changes in all biodiversity (Stevaux et al. 2009). However, this area still stands out as the region with the highest species richness and fish ecosystem services. Thus, the whole of the Paraná River basin must be recognized as a focus on the conservation of freshwater fish in this biome.
On the other hand, the high priority regions for amphibians are located in east Brazil, comprising the Serra do Mar, the Serra da Mantiqueira, the Central Corridor of the Atlantic Forest, and the Endemism Area of Pernambuco (altitude areas) (Carnaval and Moritz 2008;Carnaval et al. 2009;Campos et al. 2017;Lourenço-de-Moraes et al. 2019a, 2019bCampos et al. 2020). These important Atlantic Forest remnants have also suffered many environmental impacts from human pressures (Ribeiro et al. 2009;Rezende et al. 2018). The main factors responsible for amphibian biodiversity declines are habitat loss and fragmentation, global warming, and diseases (Young et al. 2001;Pounds et al. 2006;Ferreira et al. 2016). For amphibians, this happens because their habitat-specialist patterns in the Atlantic Forest are favored by the richness of habitats and microhabitats (Figueiredo et al. 2019), and Among assemblages, the low degree of congruence indicates that selecting one taxon as an indicator for conservation planning of the other would lead to a significant loss of species and ecosystem functions (Lawler and White 2008;Lawler et al. 2015). Considering that the selection of protected areas is often based on cross-taxon congruence among different groups (Campos et al. 2014), this strategy should be cautiously applied in conservation plans. This is because our study showed a negative result and the lack of settled down standards for this methodology appliance. Results gathered over time indicate that the patterns found are dependent on the context of the study, ranging widely with the grain, the extent, and the region of analysis .
The only region where fishes and amphibians feature high cross-taxon congruence values (≥50%) between TD and FD is the Serra do Mar and the southern part of the Central Corridor. Thus, despite the smaller number of fish species compared to the Paraná River basin, these areas can also be characterized as biodiversity hotpoints, bringing opportunities to integrate conservation strategies of terrestrial and aquatic ecosystems. Furthermore, considering that there is still much to be studied in the region (Menezes et al. 2007), the number of species that occur in these areas of the Atlantic Forest may be much larger than the one already gathered. The great concentration of biodiversity that characterizes the Serra do Mar as a zone of high conservation value is because this region harbor most of the forest remnants of the entire biome (Oyakawa et al. 2006;Ribeiro et al. 2009). The preserved forests are essential for maintaining not only amphibians' communities but also fish communities. This is because fish also depend on forest integrity for protection and food (Oyakawa et al. 2006;Menezes et al. 2007). The Serra do Mar's high preservation index is due to the slope and poor productive soils, which make it unsuitable for agricultural practices (Aguiar et al. 2003;Oyakawa et al. 2006). Therefore, despite being located between the two largest metropolitan centers in Brazil, Serra do Mar remains well preserved (Aguiar et al. 2003), containing most of the forest fragments in the entire Atlantic Forest (Ribeiro et al. 2009).
Highly correlated relationships were gathered between TD and FD within each assemblage. These results match what is usually expected for each group separately. Considering that the greater the number of species, the greater the probability of different biological characteristics, FD is expected to increase with TD (Tilman et al. 1997;Balvanera et al. 2006;Cadotte et al. 2011). This high degree of congruence among diversity measures inside each group indicates a high impact of the composition on the ecosystem processes (Tilman et al. 1997;Naeem and Wright 2003). This is because species' functional characteristics determine the contributions that a group will have to ecosystem processes, affecting their ecological performance (Cadotte et al. 2011). In this way, actions aimed at the conservation of the most species-rich areas will also protect high ecosystem values.
Conservation scientists have been continuously searching for the best methods for identifying patterns that can be successfully applied in conservation strategies (Myers 2000;Williams et al. 2005). However, when it comes to spatial relationships, isolated positive or negative results should not be used as a basis for generalization (Rodrigues and Brooks 2007;Trindade-Filho and Loyola 2011). One of the biggest challenges for using this approach is that often the strength of correlations can vary among different locations according to independent environmental factors (i.e., historical and biogeographical events, and methodological shortcomings) (Sanders 2002;Hess et al. 2006;Wolters et al. 2006;Araújo et al. 2008;Lourenço-de-Moraes et al. 2020). This high variability can potentially undermine the usefulness of this approach in conservation actions (Trindade-Filho and Loyola 2011). Considering this variation, it could gather different conclusions if only part of the Atlantic Forest territory were selected for the present assessment (e.g. in the Serra do Mar). This choice would result in different landscape management recommendations about the cross-taxon congruence patterns evaluated.
It is worth mentioning that fishes and amphibians perform varied and significant ecosystem services (e.g. energy flow regulation, nutrient cycling, and soil bioturbation; see Whiles et al. 2006;Hocking and Babbitt 2014). In addition to performing other functions directly related to the maintenance of human well-being, such as the provision of food resources, medicines, and disease control (Whiles et al. 2006;Hocking and Babbitt 2014). Thus, considering the role of these species in community functioning, their disappearance may have potentially systemic implications (Hocking and Babbitt 2014), leading to the loss of terrestrial, freshwater ecosystem services and changes in flow between these ecosystems (Holmlund and Hammer 1999;Whiles et al. 2006). These losses threaten not only human health but also their livelihoods (Egoh et al. 2009).
In conclusion, before the need to contain increased biodiversity loss, spatial relationships have been widely considered in conservation strategies (Campos et al. 2014(Campos et al. , 2017Lourenço-de-Moraes et al. 2021). However, the number of researches that show the inconsistency of generalizing the results beyond their original place of study is increasing Trindade-Filho and Loyola 2011). It is important to note that not all groups are good predictors of biodiversity, and even good predictors may be efficient for some groups and inefficient for others. In fact, this would lead to a loss of relevant information on ecosystem functioning and structure. Moreover, this study highlights the positive relationship between biodiversity and ecosystem functioning. Thus, it is expected that functions performed by different species are compromised as biodiversity loss occurs. Therefore, taken together, our results cast doubt on the use of cross-taxon congruence in a different ecosystem and from local to macro scales.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.