Assessing the effect of fish size on species distribution model performance in southern Chilean rivers

Study area and modeled species

The study area covers the Valdivia and the Bueno River basins located in the southern zone of Chile between 39.33° and 41.08°S (Fig. 1). The Valdivia River basin has a pluvial hydrological regime and is characterized by a chain of interconnected lakes at higher altitudes. The upper section of the Bueno River basin has a pluvio-nival regime, while the middle and lower parts of the basin are governed by a pluvial regime (Errázuriz et al., 1998).

Our study examined nine native freshwater fish species (Table 1): Aplochiton taeniatus Jenyns, 1842, Aplochiton zebra Jenyns, 1842, Basilichthys microlepidotus (Jenyns, 1841), Brachygalaxias bullocki (Regan, 1908), Cheirodon australe, Diplomystes camposensis (Arratia, 1987), Galaxias maculatus (Jenyns, 1842), Percilia gillissi Girard, 1855, and Trichomycterus areolatus (Valenciennes in Cuvier & Valenciennes, 1840). Statistical analyses of the effect of body size were carried out using theoretical species maximum length, which is available for all these species. Maximum length estimates were obtained from official species descriptions provided the conservation assessment of each species, developed by the Chilean Ministry of the Environment, with the only exceptions being B.  bullocki and B. australis (Table 1).

Modeling methods

Model training and evaluation

We used neural network (NNET) algorithms to estimate SDMs (Stern, 1996). This method was chosen based on their good performance with presence and absence or pseudoabsencesfor species-distribution data (Mastrorillo et al., 1997; Elith & Leathwick, 2009). NNET is derived from a simple model that mimics the structure and function of the brain and maximizes the prediction during the model-training phase by comparing actual outputs with desired outputs (Manel, Dias & Ormerod, 1999). All analyses were performed in R (v 3.5.0) using the Caret package (Kuhn, 2008).

Models were trained using a 3- or 2-fold cross-validation scheme, according to records by species (Table 1). Pseudoabsences were set at twice the number of observed species presences and were randomly selected. We designated as the final model the ensemble average of each k-fold model by its weighting each individual model by their Area Under the Curve statistic (AUC) of the Receiver Operating Characteristic (ROC).

The Caret R package was used to fit NNET models and tune the model’s two hyperparameters, namely, the weight decay for successive neural layers (“decay”) and the number of hidden units (“size”). The grid search procedure examined weight decay values ranging between 4 and 6, while the number of hidden units was allowed to vary between 0.05 and 0.9. Both hyperparameter ranges were calibrated by a trial and error process, optimizing the model performance.

Occurrence probabilities were categorized as presence/absence for all models. Thresholds were determined to maximize the sum of sensitivity and specificity (MaxSens+Spec; PresenceAbsence package in R v 3.5.0) (R Core Team, 2017). This criterion is independent of the theoretical prevalence (Manel, Dias & Ormerod, 1999; Allouche, Tsoar & Kadmon, 2006), causing the distribution of rare species to be overpredicted. In our case, the theoretical prevalence in the study area for all the species is close to 0.5, but the presence of our study species is low, requiring a relaxation of this criterion when defining the threshold that allows for the definition of each of the species distribution across the studied watersheds. All the presences and environmental data by river segment in supplemental material (Appendix S1).

Relationship between fish size and models

As proxy of body size, log₁₀—transformations of maximum length (max. length) were calculated for each species. We assessed the relationship between max. length and the number of historical records in both basins, with balanced accuracy, sensitivity, specificity, true skill statistic (TSS), and area under curve (AUC), all different performance metrics (Altman & Bland, 1994; Allouche, Tsoar & Kadmon, 2006; Velez et al., 2007; Kuhn, 2008), and with variable importance to each predictor variable. Sensitivity evaluates the proportion of actual presences that are correctly classified and specificity evaluates the proportion of actual negatives that are correctly identified (Allouche, Tsoar & Kadmon, 2006). AUC, TSS and balanced accuracy evaluate overall performance using different approximations, being valuable contrast them all (Velez et al., 2007; Kuhn, 2008). In all this cases, obtain a value of 1 to this metrics represent the best performance, and in the AUC case, 0.5 represent the worst performance. Variable importance is calculated with “Weights” method (Gevrey, Dimopoulos & Lek, 2003), and the highest importance is 100.

The statistical test of the relationship between max. length and historical records was a classical linear model (lm) (Chambers, 1992), while the other comparisons were evaluated with beta regression analysis (betareg) (Cribari-Neto & Zeileis, 2010). This test assumes response variables to be in the standard unit interval (0, 1) and to be beta distributed, and it is more precise in performance metrics cases or variable importance cases where the values are at a fixed interval. Rescaled Moment test (RMtest) and Q–Q plots were applied to all residual models to evaluate model validity (Das & Imon, 2016; Kozak & Piepho, 2018). While, the normality test is used on observed data to evaluate parametricity, in regressions the true errors are unobserved, being a common practice to use residuals as substitutes for observed data in tests for normality (Das & Imon, 2016). In this context, we used RMtest to evaluate normality of residuals and validate betareg and lm models. RM test correct shrinkage and superimposed normality effect of the residual data, improve its normality evaluation (Imon, 2003; Rana, Midi & Imon, 2009).

To avoid a potential bias from conservation status by specie on all statistical tests, we introduced this category as a covariable. Given the hierarchical relationship between conservation status categories, we transformed this to a numeric covariable, with Least Concern (LC) being 1, Vulnerable (VU) as 2, and Endangered (EN) as 3. Finally, an a priori betareg test was developed to discard the relationship between the number of presences by model and metric performance models, allowing for directly linking fish size with model indicators.

Historical records

The number of presences and conservation status did not correlate with fish size (lm test: Conservation status, coefficient = −4.7, std. error = 7.9, t value = −0.56, p value = 0.57; Historical records, coefficient = −8.1, std. error = 9.67, t value = -0.84, p value = 0.43). Model residuals show a normal distribution (Chisq = 1.83, p value = 0.59). The species with the lowest presence was A. zebra, with 22 records and 28 cm maximum length and is one of the largest native species found in the basins. The species with most historical records was T. areolatus, with 72 presences and a maximum length of 15 cm. However, many of the T. areolatus presences were spatially clustered, and the model was trained using only 36 presence records (Table 1).

Metric performance

Of all the evaluated performance metrics, only specificity showed a significant positive relationship with fish size (Fig. 2). In most modeled species, we observed residuals being normally distributed considering Q–Q plots and RMtest (Tables 3A and 3B). The only one exception is Source-of-flow: Mountains, with a non significant RMtest. Specificity had a positive relationship with fish size (Fig. 2). Other metrics did not show significant relationships with fish size (Table 3).

Variable importance

When we examined the significance of individual predictor variables, we observed that only Riparian Vegetation and Within-Channel Structures showed a significant relationship with fish size. Again, we observed residuals being normally distributed, thus validating the statistical analysis (Tables 3A and 3B). In these cases, the fitted models for larger fish tended to give less importance to these variables than models for small fish, showing a negative relationship (Fig. 2). Other variables did not show significant relationships with fish size (Tables 3A and 3B).

Records by fish species

Although many authors mention the relationship between body size and data availability (Boone & Krohn, 1999; McPherson & Jetz, 2007; França & Cabral, 2016), we found no studies in which this relationship was statistically evaluated. The premise that larger fish species are more detectable than smaller fish species depends on the sampling method and species characteristics (diurnal or nocturnal, color, habitat selection, etc.) (Boone & Krohn, 1999). Thus, a detection bias may occur if sampling is conducted at the community or ensemble level using a single sampling methodology. However, if sampling is stratified or designed for a given species population using adequate methodologies, we would expect this bias to be corrected. In this study, presence information was obtained from a governmental database generated from published scientific studies, with sampling methodologies having been determined by studies focused on independently assessing species abundances and distributions. This situation explains the lack of relationship between fish size and presence records in our study.

Fish size and model fit

Our results regarding the body-size effects on SDM performanceare relevant to the unresolved debate about this expected theoretical relationship. While Morán-Ordóñez et al. (2017) found no significant relationship between body size and model performance for trees and birds, França & Cabral (2016) successfully related model performance to species feeding mode and estuarine functional groups, with body size having a marginally significant contribution to model performance. In studies aimed at river fish, while Radinger et al. (2017) and Filipe, Cowx & Collares-Pereira (2002) found that fish size increased model performance, Markovic, Freyhof & Wolter (2012) did not observe this pattern.

Our results show that the specificity performance metric increases with fish size, suggesting that the ability to predict true absences increases with fish size. This outcome means that models for smaller fishes tend to overpredict presences. These results reflect those of Lobo & Tognelli (2011), who found using virtual species, models with randomly determined presences tend to have lower specificity relative to models with spatially biased presences. We suggest that in our study, smaller fish may have a statistically random distribution as a result of a potential mismatch between the study grain and their home range. This outcome would be the consequence of a random distribution of predictor variables used to characterize fish habitat (McPherson & Jetz, 2007).

This finding has important implications for conservation and management. First, it suggests that the environmental resolution used in the evaluation should be commensurate with the modeled species’ home range. However, when this outcome is not possible (due to financial or logistic resources available, for example), then explicit assessment is necessary to detect any potential biases in model performance for species of different body sizes. In this way, it is also important to assess both metrics that represent general performance, such as AUC, TSS or Accuracy, and to include other specific metrics that may provide relevant information, such as any likelihood of overprediction of smaller species.

Our use of theoretical maximum length as the best available estimate of species body size does not necessarily reflect actual maximum body size in the study areas, and certainly may result in biases or estimation errors. Explicit assessment of possible biases stemming from the use of theoretical or literature based maximum lengths versus empirical estimates of maximum body length would certainly improve our current understanding of the interaction between body size and SDM model performance. This could be achieved by explicitly describing observed median and maximum body size of the set of recorded presences within study hydrological basins, allowing body size to reflect current ecological factors modifying this trait. Such a study design would allow the assessment of either fine scale habitat alterations, or possible impacts of invasive species abundance. It must be noted that such a study is not free of logistic and bioethical constraints. Most important, most of the fish species addressed in this study face important threats, and wide scale sampling and individual monitoring may pose additional stress to these individuals. As a result, an important source of data for such a study could be the samples collected in the development of environmental impact assessment (EIA) studies. However, under Chilean legislation, EIA sampling designs respond the project developer specific goals, and body size sampling is not often included among the variables registered.

Participation by predictor variable

The predictor variables that show trends in model importance according to fish size are both obtained at smaller scales (reach scale) (Frissell et al., 1986). The participation of riparian Vegetation and Within-Channel Structures decreases with fish size, suggesting that these variables are not relevant to estimating the ecological niche of large fish species. For Within-Channel Structures, Radinger et al. (2017) found that these may cause larger habitat reduction for larger fish species by creating stream barriers that limit their dispersal. This findingis consistent with our results, which show how the largest three species have low occurrence probabilities on reaches with almost one structure detected.

The relationship between riparian vegetation cover and fish size may be determined by the association between riparian vegetation and river slope. Rivers with higher slopes have short alluvial plains that allow for more riparian vegetation (Stefunkova, Neruda & Vasekova, 2019). However, the morphologies of Chilean fish species are not adapted to high-slope habitats (Arratia, 1987; Arratia, 1992); therefore, we presume that riparian vegetation acts as a proxy of habitat conditions adequate for smaller species. Functionally, this notion means that the occurrence probability of smaller species increases in rivers with less riparian vegetation cover.

The current study presumed that no interspecific effects occurred between species that impacted SDMs, particularly the potential role of the invasive species. Invasive species have been shown to influence native fish body size (Blanchet et al., 2010), which can have significant effects on SDM performance. In Chile, invasive species such as rainbow trout (Oncorhynchus mykiss) and brown trout (Salmo trutta) have strong impacts on fishes (Pardo, Vila & Capella, 2009; Arismendi et al., 2014) and may benefit from riparian vegetation conservation, especially in upland areas (Lacy, Ugalde & Mao, 2018). Regrettably, in the current study, we were unable to evaluate the effect of invasive species on body size. Future studies should seek to integrate ecological effects, especially effects caused by invasive species, into the development of SDMs.

Finally, we want to highlight the use of source-of-flow as a predictor variable in our study; we have not found this use in prior researchon river-species modeling, and it is especially important in torrential river systems, such as those found in Chile. These systems have short runs, with relatively large lakes, glaciers, or salt pans that significantly affect hydrological and hydraulic conditions. The source-of-flow variable is implemented in river evaluations in New Zealand (Snelder & Biggs, 2002) and Chile (Peredo-Parada et al., 2011), which would facilitate its use in SDMs.