Predicting genetic biodiversity in salamanders using geographic, climatic, and life history traits

RESEARCH ARTICLE Predicting genetic biodiversity in salamanders using geographic, climatic, and life history traits Danielle J. Parsons1,2, Abigail E. Green3, Bryan C. CarstensID1,2*, Tara A. Pelletier3 1 Museum of Biological Diversity, The Ohio State University, Columbus, Ohio, United States of America, 2 Department of Evolution, Ecology and Organismal Biology, The Ohio State University, Columbus, Ohio, United States of America, 3 Department of Biology, Radford University, Radford, Virginia, United States of America * carstens.12@osu.edu Abstract The geographic distribution of genetic variation within a species reveals information about its evolutionary history, including responses to historical climate change and dispersal ability across various habitat types. We combine genetic data from salamander species with geo- graphic, climatic, and life history data collected from open-source online repositories to develop a machine learning model designed to identify the traits that are most predictive of unrecognized genetic lineages. We find evidence of hidden diversity distributed throughout the clade Caudata that is largely the result of variation in climatic variables. We highlight some of the difficulties in using machine-learning models on open-source data that are often messy and potentially taxonomically and geographically biased. Introduction Documenting biodiversity is an important first step in understanding both ecological and evo- lutionary processes [1], particularly the functional roles that act to connect processes function-

Editor: Christopher Nice, Texas State University, UNITED STATES OF AMERICA Received: February 22, 2024 Accepted: September 9, 2024 Published: October 18, 2024 Copyright: © 2024 Parsons et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: Data are uploaded to Dryad and can be accessed using the following link: https://doi.org/10.5061/dryad.m63xsj474. Code related to this manuscript, including data cleaning, imputation, predictive modeling, and significance testing has been deposited in GitHub (https://github.com/parsons463/ HiddenSalamanders). All remaining data are available in the manuscript and/or Supporting information.

or endangered (P. stormi). More recently, Parra Olea et al. [11] discovered five cryptic lineages in Chiropterotriton from Mexico, several of which are threatened due to their restricted ranges [11]. Species with small ranges and/or limited dispersal capabilities can be harder to protect because their distributions often do not fall within protected areas [12] and small ranges are often used as a factor in assigning conservation priorities [13]. Therefore, it is important to identify these hidden lineages, as they could easily go unnoticed and unprotected. Many other species of salamander that would have otherwise gone unnoticed and have been recognized using molecular data have small ranges and likely need protection [14–18]. The presence of cryptic diversity has been recently highlighted as a key component of undescribed biodiversity that requires greater attention [19,20]. Efforts to conserve undescribed genetic diversity can be facilitated using computational methods that identify genetic lineages representing potentially hidden diversity in need of fur- ther investigation. The use of data science techniques has allowed biodiversity studies to expand their geographic and taxonomic focus to explore broader patterns of evolution, which can be difficult to assess using traditional meta-analysis methods [21]. Macrogenetics, a rela- tively new field that merges biodiversity data with genetic data [22,23], has been used to explore how human impacts influence levels of intraspecific genetic diversity [24,25], to study

past and future climate refugia [26,27], and to quantify latitudinal biodiversity gradients [28– 31]. Macrogenetic methods, particularly in combination with predictive modeling, can be used to inform conservation policies by identifying species, taxonomic groups, or geographic areas in need of further investigation [32,33]. Machine learning models excel in processing and analyzing large datasets automatically, making them highly efficient at identifying com- plex, non-linear patterns within extensive data. While these models are powerful tools for making predictions with high accuracy, the process behind many machine learning models can be difficult to interpret [34]. Random forest models, a type of predictive machine learning model, enhance interpretability by providing clear insights into feature importance, thus rep- resenting a powerful tool for identifying factors influencing biodiversity. Recently, such analy- ses have expanded to taxonomic work. Parsons et al. [35] analyzed mitochondrial DNA sequences from over 4000 species of mam- mals, representing roughly 66% of currently described species, and found that mammal diver- sity is largely under-described using molecular species delimitation methods on publicly available barcode data. This is useful for several reasons. A comprehensive list of undescribed genetic lineages that may represent species now exists that can help focus taxonomic efforts. Parsons et al. [35] also found that taxa with small bodies, and large geographic distributions

with variation in precipitation and isothermality, were more likely to contain cryptic diversity. While some of this might seem obvious (morphological differences are harder to observe in small-bodied animals and these animals may be harder to find), it does allow researchers to document characteristics of species, higher taxonomic groups, or even geographic regions that contribute to diversification and therefore biodiversity patterns. When done in disparate taxo- nomic groups (e.g., vertebrates, invertebrates, plants, and fungi) and at different levels (e.g., Class, Order, Family) this furthers our understanding of core evolutionary processes. A similar approach was taken in birds. Using a tree-based molecular species delimitation method, Smith et al. [36] found that latitude explained variation in phylogeographic breaks, while other traits pertaining to habitat and life history explained very little. In this case, phylo- geographic structure was higher in the tropics. Conversely, in other organisms, isolation-by- distance within species is often higher at higher latitudes (multiple taxonomic groups: [32]; amphibians: [37]). Further, genetic variation within amphibians was best explained by range size and elevation, rather than latitude, in the neotropics [37], while latitude was an important predictor of genetic diversity in the nearctic [30]. This suggests that differences exist in how PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 2 / 20 Funding: This work was supported by the

Directorate for Biological Sciences (DEB-1910623 to BCC; DEB-1911293 to TAP). Competing interests: The authors have declared that no competing interests exist.

genetic variation is distributed within species depending on which taxonomic groups are being examined, and at what spatial scale. In order to expand these approaches, we conducted an assessment of genetic lineages in roughly 100 described salamander species using the phylogatR database [38]. PhylogatR aggre- gates DNA sequence data from both GenBank and BOLD into sequence alignments, providing associated GBIF occurrence records (i.e., GPS coordinates) for each sequence. There are over 700 nominal species of salamanders belonging to nine families [39], most located in the north- ern hemisphere. While salamanders contain a wide variety of life history strategies and habi- tats, they are likely to have high levels of cryptic diversity due to their moisture requirements and similar body forms. However, their eco-evolutionary processes can vary from species to species and sometimes oppose our expectations [40–45]. We follow methods from Parsons et al. [35] and use molecular species delimitation methods to estimate the number of genetic lineages present in previously collected data that is both openly available and easily tractable. We utilize these delimitation results to identify species that are likely to harbor undescribed diversity. Species for which delimitation reveals multiple genetic lineages are classified as hid- den species. The individual genetic lineages that comprise these hidden species are referred to as hidden genetic lineages. We then use a random forest classification to determine whether

any variables pertaining to geography, the environment, or life history traits contribute to the presence of hidden genetic lineages within species. We also discuss some of the difficulties in using open-source data that are often messy and potentially taxonomically and geographically biased. Materials and methods Collection of genetic and geographic data We downloaded all available data from the phylogatR database (https://phylogatr.org/) using the search term ‘Caudata’ on 2/4/22. The uncleaned data represented four families, 93 different species, and 14 loci with a total of 3768 DNA sequences. To begin cleaning the data, we calcu- lated nucleotide diversity (pi) values for each locus in every species and found outliers by set- ting lower and upper bounds of 2.5% (0) and 97.5% (0.2193634) respectively. For each of the four outliers and two species with missing pi values, we opened the DNA sequence file in Mes- quite v3.7 [46] and removed any extremely short or non-overlapping sequences (S1 Data). Additionally, we discovered a typo for the species Batrachuperus karlschmidti causing there to be two different species folders for the same species. Both the sequence and occurrence files were merged for the species and the sequence files were realigned to correct the error. Two species complexes were present in the dataset, and these were kept named as downloaded: Tri- turus cristatus x dobrogicus macrosomus and Ambystoma laterale jeffersonianum complex. A review of the available loci indicated that two genes, COI and cytb, were the most well-repre-

sented in both total number of sequences and species coverage. Consequently, we opted to uti- lize these two genes for downstream analysis. Species alignments from the download for both the mitochondrial genes Cytochrome oxi- dase I (COI) and Cytochrome b (cytb) were merged for all salamander species and aligned using MAFFT v7.5 [47] with the default settings and including the–adjustdirection command to account for reverse complement sequences. We visually inspected alignment files for both genes and removed all short sequences, which we classified as those missing 50% or more of the second half of the sequence. Twenty-one sequences were removed from the COI alignment and 99 were removed from the cytb alignment, leaving totals of 768 and 908 sequences for COI and cytb, respectively. The sequences for seven species were completely removed from further analysis due to their short length (missing 50% or more of the second half of the sequences). In PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 3 / 20

total, eighty-three species remained with an average of approximately 20 sequences per nomi- nal species (see S2 Data for a list of identifiers corresponding to the sequences used in this study). Species delimitation We used three methods of species delimitation to determine the number of genetic lineages present in our samples. The GMYC is a tree-based method that takes a phylogenetic tree as input and finds a point in the tree where branching changes from within to between species [48]. The ABGD [49] and ASAP [50] methods are distance-based delimitation methods that use pairwise genetic distances to establish the threshold between intra- and inter-species diver- gence. Because each method is based on a specific set of assumptions, it is best to use multiple methods and compare their results in order to achieve a more accurate delimitation [51]. By looking for concordance across methods, we can increase our confidence in the identified line- age boundaries and minimize the potential impact of bias introduced by any single method. While we report delimitation results from the genes COI and cytb for all methods, we used a consensus of these results—reflecting agreement among the GMYC, ABGD, and ASAP delimi- tation methods for both COI and cytb—for assessing the influence of geography, environment, and life history traits on predicting salamander genetic diversity. To estimate a species tree for input into the GMYC, we used BEAST v2.5.1 [52]. We used the default parameters except for conducting 100,000,000 million generations, sampling every

5,000, and setting the model of sequence evolution to GTR+I+G [53]. The log files were checked by eye using Tracer v1.7.2 [54]. ESS values were all over 1000 for both cytb and COI. We removed 10% as burnin and retained the maximum clade credibility tree using TreeAnno- tator. After checking that the tree was binary and ultrametric, we used the R package splits [55] to conduct GMYC analyses. In each case we used the single threshold model and all other default settings. We conducted both ABGD and ASAP delimitation analyses via their web por- tals (https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html and https://bioinfo.mnhn.fr/abi/ public/asap/asapweb.html, respectively) using the default parameter settings. Predictor variables For each nominal salamander species, we examined numerous geographic, environmental, morphological, and life history variables to identify traits predictive of undescribed salamander diversity using a classification model based on our species delimitation results. A variety of predictor variables were collected, including geographic and environmental values derived from georeferenced locality data (see S3 Data). In addition, three life history traits were avail- able from AmphiBIO, a global database for amphibian ecological traits [56], for most of the species in our study: reproductive strategy (direct developing, larval phase), habitat (terrestrial, fossorial, aquatic, or some combination of these), and body size (total length). To supplement this dataset and fill in any missing trait values, we used AmphibiaWeb [57] and other online

sources (S4 Data). To extract species specific data related to its environmental distribution, we utilized 42 GIS data layers (see S4 Data for data layer details), meant to capture various aspects of the ecology and habitat of each species. These include all 19 BIOCLIM layers from the CHELSA database [58,59] at 1 km resolution, elevation [60], terrestrial habitat heterogeneity [61], global land cover classification [62], global river classification [63], disaster risk [64], and various indica- tors of seasonal growth 58–59]. In addition to traits meant to capture various ecological factors we also gathered data for several traits relating to human impact, in order to measure levels of human disturbance and activity to the species environment. While ecological factors directly PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 4 / 20

affect levels of biodiversity by influencing species biology, anthropogenic factors can influence the way we find and describe this diversity (e.g., increased sampling in wealthier, more popu- lated locations). We extracted species specific data from several GIS layers, including anthro- pogenic biome [65], human population density [66], and gross domestic product [67] in order to evaluate how anthropogenic factors impact undescribed diversity. We utilized the R packages ’raster’ [68], ’rgdal’ [69], ’geosphere’ [70], and ’plyr’ [71] to extract species specific information from each layer using geographic occurrence records obtained from phylogatR. To represent the environmental variation within the occupied range of each species, we extracted the value of each environmental layer for each GPS coordinate associated with each species. We then took the mean and standard deviation for each environmental variable. To obtain species specific data related to geographic distribution we extracted the minimum, maximum, mean, and length of latitude and longitude from the GPS points of each species. We used the R package ‘mice’ [72] to impute trait values missing from our dataset (see S1 Fig in S1 File for distribution of missing data and specific trait values imputed). The imputa- tion method ‘pmm’ was used for all numeric variables and ‘polyreg’ was used for categorical variables (i.e., reproductive strategy and habitat). We ran the imputation 15 times (S2 Fig in S1 File) and then pooled the iterations to generate the final imputed values. The final database

containing all trait values (both imputed and original) is available in S4 Data. Predictive modeling We used the R package ’caret’ [73] to generate a random forest classification model [34] based on our previously generated database of predictor variables and a consensus of our species delimitation results. Two separate sets of consensus models were generated to assess the role of geography, environment, and life history traits on the presence of hidden diversity (Fig 1A). The first model (all agree) represents a strict consensus of delimitation results from species in which results from all methods of species delimitation agree (Fig 1B). Any species with con- flicting delimitation results were excluded from analysis. The second model (majority rules) represents a majority rule consensus in which species are assigned to a response category based on relative support of delimitation results (Fig 1C). For each model, we used 70% of the data to train the model and the remaining 30% was set aside as a test set. Models were gener- ated using 10-fold cross validation with five repeats to tune the parameter ’mtry’, the number of variables randomly sampled at each split, and optimize the area under the receiver operating characteristic curve, ROC. After training, we extracted the variable importance measures mean decrease accuracy (MDA) and Gini impurity (Gini) from the final models. We then used the final models on the test set data to evaluate model performance. Model performance was evaluated across a variety of metrics including model accuracy, which reflects how well the

predicted classifications agree with the observed classifications, and both positive and negative predictive value, which indicate the how the model performs on observations from each class. Additionally, we calculated the no information rate (NIR), the proportion of observations that fall into the majority class, and the p-value [Accuracy>NIR], to test for model significance. The top important predictor variables from our best model were compared using a Kruskal- Wallis test to determine if these variables are significantly different between species that do or do not contain hidden diversity. Results Genetic and geographic dataset Our final dataset consisted of 1676 DNA barcoding sequences (Fig 2). Of these, 768 sequences were from the Cytochrome oxidase I gene (COI), and 908 sequences were from the PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 5 / 20

Cytochrome b gene (cytb). These sequences were derived from 83 nominal species of salaman- ders, which were distributed among 26 distinct genera occurring across the globe. The dataset contained 13 species with sequences from the gene cytb. Comparatively, COI exhibited notably broader taxonomic coverage, with 77 nominal species represented. Out of the 83 species ana- lyzed, only seven were shared between COI and cytb. Of the remaining 76 species, 70 were unique to COI and six were unique to cytb. To supplement the genetic data collected, a total of 1676 georeferenced occurrence records from phylogatR were utilized to collect a combination of geographic, environmental, and life history trait values for each nominal species present in the dataset. Species delimitation and consensus assignment Species delimitation results were generated by analyzing COI and cytb sequences from each nominal species under three different delimitation methods, ABGD, ASAP, and GMYC. We classified each nominal species as either containing undescribed genetic lineages or not con- taining undescribed genetic lineages based on the number of genetic groups predicted by each delimitation analysis. While taxonomic overlap between COI and cytb was narrow, delimita- tion results for species shared by both loci were mostly congruent with respect to species classi- fication. Of the seven species with sequences from both genes, only two species produced conflicting results regarding the presence of undescribed genetic lineages within a specific

Fig 1. Consensus classification of species delimitation results. A, Flowchart describing the process of generating a consensus of delimitation results (among different methods and loci). B, C, Pipeline for classifying nominal species as either containing or not containing hidden genetic diversity in each consensus analysis (all agree and majority rules, respectively). https://doi.org/10.1371/journal.pone.0310932.g001 PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 6 / 20

taxon based on loci. Delimitation results across different methods showed slightly less agree- ment. Classifications resulting from the GMYC and ASAP methods were similar across spe- cies. These methods, on average, resulted in slightly fewer predicted species per nominal species than the ABGD method (see Fig 3 for predicted species numbers). To account for this variation in our final random forest classification models, we generated two consensus classifications to evaluate concordance between delimitation results from dif- ferent methods and loci. The results of our consensus models indicate that roughly 2/3rds of the nominal salamander species used in this analysis are likely to contain genetic lineages that may be unexplored diversity. The strictest of these classifications produced a consensus model (all agree) consisting of 51 total species, 41 of which were classified as containing hidden diver- sity and 10 of which were classified as not containing hidden diversity. The remaining consen- sus model (majority rules) consisted of 83 total species, of which 51 were classified as containing undescribed genetic lineages and 32 were not (Fig 3). Predictive modeling For our majority rules and all agree consensus classifications, we developed random forest clas- sification models using all available predictor data. To assess potential correlation between var- iables in our dataset we used the R package ’corrplot’ [75] to generate a correlation matrix of

Fig 2. Geographic spread of salamander data. Map shows geographic distribution of salamander occurrences pulled from phylogatR [38] and used in these analyses. Pie charts show the total number of cytb and COI sequences used (left) and the number of species represented by those cytb and COI sequences (right). Basemap created with world map data from the public domain Natural Earth project (http://www.naturalearthdata.com). Salamander figures in black were obtained from Phylopic [74] and are licensed under public domain. https://doi.org/10.1371/journal.pone.0310932.g002 PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 7 / 20

Fig 3. Species delimitation results. A, Graphs show the results of ABGD, ASAP, and GMYC species delimitation analyses of the genes cytb and COI for each nominal species. Numbers represent the predicted genetic lineages from each analysis. Results highlighted in red indicate no hidden genetic lineages were predicted (i.e., number of genetic lineages = 1). Results highlighted in green indicate hidden genetic lineages were predicted (i.e., number of genetic lineages > 1). Grey highlighting indicates that specific analysis was not performed due to a lack of data. B, Pie charts display the number of nominal species classified as either containing or not containing hidden diversity in each consensus analysis (i.e., all agree and majority rules). https://doi.org/10.1371/journal.pone.0310932.g003 PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 8 / 20

our predictor variables (S3 Fig in S1 File). Due to the presence of strong correlations between several of the geographic and environmental variables in our dataset we performed multiple random forest models with progressive sets of correlated variables removed at different cutoff values (i.e., |correlation coefficient| > 0.75; 0.85; 0.9). The results of these random forest mod- els are presented below (Table 1). All random forest models were found to have high predictive accuracy, with the majority rules and all agree models achieving accuracies of 75–85% and 87–93%, respectively, in identi- fying nominal species likely to contain hidden diversity. Although these results may initially seem to suggest that all our models are able to make meaningful predictions, further examina- tion of additional model evaluation metrics reveals potential overfitting and inflation of pre- dictive power. For example, despite the high accuracy of the models, the 95% confidence intervals for these values are broad with an average length of nearly 40% for most of the models (Tables 1 and 2). Additionally, the no information rates (NIRs), a measure of prediction signif- icance based on the underlying dataset that needs to be exceeded in order for model results to be significant, are particularly high for the all agree consensus models, where the class frequen- cies are more skewed towards species predicted to harbor hidden diversity. The high NIR val- ues combined with wide confidence intervals result in a p-value [Accuracy > NIR] greater

than 0.05 in all models, except for the majority rules consensus using a correlation cutoff of 0.90. While all our models show high accuracy, when the additional model evaluation metrics are considered only one has strong predictive power. Therefore, we only used the majority rules consensus using a correlation cutoff of 0.90 for interpreting variable importance of our data. Evaluation of variable importance We extracted variable importance measurements from each random forest classification using the variable importance metrics MDA and Gini. While there was some overlap of top predictors between different models (Fig 4; S4 Fig in S1 File), no specific predictors were Table 1. Results of majority rules consensus random forest models. Model metrics for each random forest predictive model generated using the majority rules consen- sus classifications are shown. Majority Rules Models Original |Correlation| > 0.75 |Correlation| > 0.85 |Correlation| > 0.90 Accuracy 0.75 0.75 0.75 0.8333 Accuracy (95% CI) (0.5329, 0.9023) (0.5329, 0.9023) (0.5329, 0.9023) (0.6262, 0.9526) No Information Rate 0.625 0.625 0.625 0.625 Pos Pred Value 0.7368a 0.8 0.7647 0.7895 Neg Pred Value 0.8 0.6667 0.7143 1 P-Value [Acc > NIR] 0.1453 0.1453 0.1453 0.02435 https://doi.org/10.1371/journal.pone.0310932.t001 Table 2. Results of all agree consensus random forest models. Model metrics for each random forest predictive model generated using the all agree consensus classifica- tions are shown.

All Agree Models Original |Correlation| > 0.75 |Correlation| > 0.85 |Correlation| > 0.90 Accuracy 0.8667 0.9333 0.8667 0.8667 Accuracy (95% CI) (0.5954, 0.9834) (0.6805, 0.9983) (0.5954, 0.9834) (0.5954, 0.9834) No Information Rate 0.8 0.8 0.8 0.8 Pos Pred Value 0.8571 0.9231 0.8571 0.8571 Neg Pred Value 1 1 1 1 P-Value [Acc > NIR] 0.398 0.1671 0.398 0.398 https://doi.org/10.1371/journal.pone.0310932.t002 PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 9 / 20

Fig 4. Variable importance for random forest classification models generated using the majority rules consensus. Variables ranked among the top ten most important variables (based on MDA and Gini) from the classification model generated at different correlation cut-offs are included. Blue highlighting indicates the best consensus model (majority rules–correlation cutoff 0.90). https://doi.org/10.1371/journal.pone.0310932.g004 PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 10 / 20

consistently predicted to be of significantly higher importance than other predictors in the model. Instead, importance was split across numerous predictors that were found to be unstable between models. This instability supports previous indications that many of the random forest models are likely prone to overfitting. Despite the lack of a strong set of stand- out predictors across models, one pattern does emerge that is applicable to the species in our dataset. Of the top ten most important predictors in each model, approximately 85% are measurements of standard deviation (vs. measurements of mean values or life history traits) (S5 Data). This is supported by further examination of our one model that was able to predict significantly better than random, the majority rules consensus with a correlation coefficient cutoff of 0.90, in which the top five most important predictors are measurements of standard deviation. Significance testing indicates that species identified as containing hidden genetic lineages often have ranges characterized by a larger variance in annual and seasonal precipi- tation, isothermality, and net primary productivity than species not identified as harboring hidden genetic lineages (Fig 5). Fig 5. Comparison of hidden vs not hidden trait values for the top five most important predictors of the best consensus model (majority rules– correlation cutoff 0.90). A, Columns 1–2 of the table identifys the specific model and predictors (i.e., traits). Columns 3–4 show the median trait values for

each group (i.e., hidden vs not hidden). Columns 5–6 show the results of Kruskal-Wallis significance tests, which determine if the difference in median trait values for each group is statistically significant. B, Corresponding boxplots of the median trait values for the top five most important predictors show a significant difference in the range of values between hidden and non-hidden genetic lineages. https://doi.org/10.1371/journal.pone.0310932.g005 PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 11 / 20

Discussion When identifying genetic lineages or delimiting species, it is important to recognize that spe- cies concepts are complex and often differ based on various factors, such as geographic loca- tion, reproductive isolating mechanisms, genetic markers, and taxonomic practices. Therefore, it is essential to approach species delimitations with caution and to recognize that they repre- sent a hypothesis or starting point rather than a definitive answer [76]. In addition, while mito- chondrial data can be suitable for preliminary assessments of species diversity [77], these assessments should be considered in tandem with other species information and relevant data when describing species boundaries. However, with recent advances in technology rapidly increasing the quantity of publicly accessible genetic and geographic datasets, these data offer a cost effective and efficient way to explore large-scale patterns and predictors of intraspecific genetic variation (e.g., [24,29,78]). Our results suggest that there are undescribed genetic lineages that may warrant further investigation distributed within Caudata. Adequately documenting biodiversity, both at the species and population level, is a first step in understanding the eco-evolutionary processes generating this diversity. However, in most clades, the Linnean shortfall is likely to influence broad scale patterns detected using macrogenetic approaches [13], making it essential to con- sider how the taxonomic designations used to inform these approaches influence the patterns

detected. This is particularly important when dealing with clades suspected of harboring high levels of cryptic diversity. For example, Miraldo et al. [24] generated the first global map of genetic diversity within species of mammals and amphibians. One of their main conclusions was that amphibians displayed lower levels of genetic variation in areas with higher human impact. Similarly, in amphibians, several recent studies have found within species genetic diversity to be lower in temperate regions in species with smaller ranges and at higher eleva- tions [30,37]. The methods used to detect these patterns are based on current taxonomic knowledge, and as such, rely on the assumption that the species designations used are accurate. However, if species descriptions inaccurately reflect biological diversity, nominal species that contain cryptic species will display higher levels of genetic diversity, while not reflecting true within species variation, potentially skewing our interpretation of any patterns that result. Evaluating support for identified genetic lineages While our delimitation of genetic lineages are a starting point, or hypothesis generation step, for evaluating a species in nature where complex processes, such as hybrid zones, and adequate sampling must be considered [e.g., 75–77], we believe these computational approaches are use- ful for targeting species in further need of examination. We conducted a literature search to explore whether the nominal species in our dataset have been previously explored from a spe-

cies delimitation approach. We used the online American Museum of Natural History taxo- nomic and nomenclatural database, Amphibian Species of the World [79]), to evaluate current taxonomic research in each nominal species of salamander predicted to contain hidden diver- sity in our consensus model. Species in which we were able to identify research-based support for the potential of undescribed diversity were recorded, along with the related articles in which the diversity was described as well as the type of data used (see S7 Data). Nearly 70% of species the majority rules consensus suggests harbor hidden lineages, contain results that also support the potential splitting of species into separate lineages. Out of these about 38% were explored using mt DNA only, 10% with nuclear DNA only, 35% using a combination of both nuclear and mt DNA and 17% using mt DNA, nuclear DNA and morphology. Just under 10% of the species display a complex history of hybridization, making delimitations difficult, a situ- ation not uncommon in salamanders [e.g., 44,80,81]. We were unable to find results for PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 12 / 20

roughly 25% of our species data. We encountered 5 species in which the results of previous delimitation work were either unclear or considered highly contested (e.g., Ichthyosaura alpes- tris, Batrachuperus karlschmidti, Batrachuperus taibaiensis, and Salamandrella schrenckii). Taxonomy is dynamic field [33] and given our search, it can be difficult to use current open- source data relying solely on species names. However, the current literature largely supports the delimitation results found here and suggests a number of species in further need of investi- gation (see citations in S7 Data, formal name changes, and an ability to update current open- source databases to reflect these changes). Additionally, even though there are limitations to using current open-source data that might not keep up to date with current taxonomy, we can still determine what factors might predict the presence of species likely to possess undescribed genetic diversity. Significant predictors of diversity Significance testing of the most important predictors from our best model (majority rules con- sensus with a correlation coefficient cutoff of 0.90) indicates that the species which our analysis identified as containing hidden genetic lineages often have ranges characterized by a larger variance in annual and seasonal precipitation, isothermality, and net primary productivity when compared to species that were not identified as containing hidden genetic lineages by our analysis (Fig 5B). And while the order of the most important traits is unstable across differ-

ent models, across all models most of the traits found to be important were measurements of standard deviation (vs. measurements of mean values or life history traits) (S5 Data). This sug- gests that the presence of variation in climate, rather than any species-specific trait or charac- teristic is the most identifiable driving force of within species genetic diversity for salamanders at this scale. Our findings align with similar studies of amphibians using a different measure of genetic variation within species (nucleotide diversity), which concluded that species traits were not a predictor of intraspecific genetic diversity [30,37]. Using similar methods, our results in salamanders differ from that found in mammals, where body size and range size were the most important predictors [35]. These findings are somewhat consistent with other studies of salamander diversification. Reproductive mode (larval stages, direct development) and habitat (combinations of terres- trial, aquatic, arboreal) vary across species and have evolved multiple times but have not been found to directly correlate with speciation, though being a direct developer might increase diversification rates [82]. In vertebrate clades, terrestrial organisms tend to have higher diversi- fication rates than aquatic organisms [83], but we did not have a large number of fully terres- trial species in our dataset, which might have limited our ability to detect this as an important predictor. Alternatively, in one species which has intraspecific variation in habit, Salamandra

salamandra, terrestrial-breeding individuals exhibited greater geographic genetic differentia- tion compared to aquatic-breeding individuals [84]. Not surprisingly, this species showed con- flicting results in our delimitation analyses. Because various species delimitation methods are not similarly sensitive to differing levels of population structure, we would expect these meth- ods to perform more inconsistently within species with highly variable genetic and geographic distance across different life histories [80,85]. Given that salamanders are relatively constrained in body form and ecological niches, varia- tion in climatic variables seems like a reasonable explanation for species containing cryptic diversity. This follows the suggestion that change in climatic niche variables increases diversifi- cation rates in plethodontid salamanders [86]). Diversification rates in frogs and salamanders have been shown to be higher near the tropics [83], so one might expect latitude to be an PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 13 / 20

important predictor. However, latitude was not included in the list of predictor variables that were likely to be important (Fig 4). Predictive modeling as a tool to address the Linnean shortfall Recently, Parsons et al. [35] used publicly available genetic barcoding data to develop a predic- tive framework to identify mammalian clades most likely to contain hidden species and deter- mine specific trait complexes that indicate where hidden mammal diversity is likely to exist. We adopted a similar approach to evaluate undescribed genetic lineages in the clade Caudata, a group which differs from mammals in several key aspects, including species richness and sampling intensity. We focused on a lower taxonomic level so there are fewer recognized spe- cies of salamanders (<1000; [57]) compared to the mammal dataset, making the ability to pro- duce robust predictive models more challenging. Additionally, there was a smaller proportion of available data for salamanders than mammals (~10% compared to 60% of described spe- cies). However, these smaller datasets might be more realistic in that they are more representa- tive of the type of data most likely to be available for the taxonomic groups that are in greatest need of attention from taxonomists. While the random forest models generated in this study actually have a higher overall accu- racy than those used in Parsons et al. [35] (see Table 3), relying on this metric alone to evaluate the performance of predictive models can be misleading [87–89]. For classification models,

model accuracy depends on how well the predicted classifications match the observed classifi- cations. While seemingly straightforward, accuracy does not account for other model charac- teristics that may be influencing model behavior, such as the class frequencies of the underlying dataset [87]). In cases where one class occurs at a much higher frequency than the other, a predictive model can attain a high accuracy by simply always predicting the higher class. Therefore, an important benchmark to consider when interpreting overall model accu- racy is the frequency at which the majority class occurs, the no information rate (NIR) [88]. If a model’s accuracy is not significantly higher than the NIR (i.e., p-value [Accuracy > NIR]), it can remain unclear whether the model is making meaningful decisions. In our models, the overall accuracy was found to be high, but the 95% confidence intervals for the accuracy values are very wide for most of the models. In addition, because the dataset is skewed towards spe- cies classified as containing hidden diversity, the p-value [Accuracy > NIR] was found to be significant in only one model. This is important to point out because even though there are large datasets available, choosing the right analytical tools can remain challenging depending on the use of the predictive models. Beyond analytical tools, it’s also important to consider your dataset, and how the characteristics of your dataset are affecting the results you obtain. Considering the scale of not only the dataset, but also the analytical methods used and the

Table 3. Summary of results of mammal random forest classification models presented in Parsons et al. (Parsons et al., 2022 [35]). Model metrics for each random forest classification model generated using data from the class Mammalia are shown. Mammal Models ABGD COI ABGD cytb GMYC COI GMYC cytb consensus Accuracy 0.737 0.68 0.6429 0.6517 0.781 Accuracy (95% CI) (0.6802, 0.7885) (0.6333, 0.7241) (0.5821, 0.7004) (0.6014, 0.6996) (0.7273, 0.8285) No Information Rate 0.7222 0.6235 0.6128 0.5488 0.6533 Pos Pred Value 0.56667 0.6304 0.17271 0.6624 2.85E-06 Neg Pred Value 0.75833 0.6937 0.5571 0.6345 0.807 P-Value [Acc > NIR] 0.32 0.008792 0.6735 3.00E-05 2.85E-06 https://doi.org/10.1371/journal.pone.0310932.t003 PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 14 / 20

pattern one is attempting to examine is especially important in meta-analyses, as different pat- terns emerge at different scales [89]. Conclusions Here, we chose to utilize biodiversity data from phylogatR (i.e., genetic data for which directly associated specimen locality information is available) to avoid potential discrepancies between the distribution of the genetic and geographic data analyzed. By doing so we hoped to gain a more fine-grain understanding of how species genetic diversity is influenced by geographic and environmental factors [23]. However, making this choice significantly decreased the amount of data available and led to a greatly reduced dataset. Our study included 1676 DNA barcoding sequences from the genes COI and cytb (768 and 908 sequences each, respectively). However, a 3/31/23 search of GenBank for salamander barcoding sequences from the genes COI and cytb returned a total of 17097 sequences (4468 and 12629 sequences each, respec- tively; see S6 Data). Similarly, while we were able to obtain 1676 occurrence records tied to the genetic sequences used in this study, a GBIF search for geographic occurrences tied to sala- mander preserved specimens and material samples returned 675243 records (see S6 Data). This study highlights the lack of genetic data with easily-associated geographic information. Despite limitations in dataset size and geographic coverage, our framework effectively iden- tified salamander species likely to contain undescribed genetic diversity with agreement across

multiple delimitation methods. These species likely represent good candidates for further taxo- nomic evaluation. While we were unable to pinpoint a specific predictor variable as the most important for predicting undescribed diversity, our findings suggest that hidden diversity in salamanders is likely higher in species with broad geographic ranges characterized by signifi- cant climatic variability. This insight serves as a starting point for future integrative taxonomic work and underscores that much diversity remains undiscovered. Although our study was constrained by data availability, the framework we used could further elucidate these relation- ships with access to more comprehensive genetic and geographic data, highlighting the crucial importance of such data in biodiversity studies. The numerous benefits of making biological data more broadly available have been repeat- edly demonstrated [90]. And recent years have seen a significant increase in the amount of available specimen and biodiversity data. The utility of these data to address large scale pat- terns of biodiversity, such as those examined in this study, is enhanced by our ability to inte- grate and synthesize data across different data sources, types, and taxonomic groups [91]. Our study highlights the importance of not just making these data available but making them avail- able in a way that is standardized and will facilitate integration and re-use for future genera- tions to come (e.g., [92,93]). Supporting information

S1 File. (Contains Figure S1: Distribution of missing data in the salamander trait database; Figure S2: Distribution of imputed trait data; Figure S3: Correlation matrix of predictor variables; Figure S4: Variable importance for predictive models). (DOCX) S1 Table. Comparison of salamander data available for each loci on the phylogatR data- base. (DOCX) PLOS ONE Predicting hidden lineages in salamanders PLOS ONE | https://doi.org/10.1371/journal.pone.0310932 October 18, 2024 15 / 20

S2 Table. Comparison of model accuracy confidence intervals between salamander and mammal predictive models. (DOCX) S1 Data. Nucleotide diversity of Caudata sequences from phylogatR. (XLSX) S2 Data. PhylogatR identification numbers (accession_sourceID; see Pelletier et al., 2022 [38]) for records analyzed. (XLSX) S3 Data. Final dataset of response and predictor variables. (XLSX) S4 Data. Variable specifics and source information. (XLSX) S5 Data. Variable importance extended results. (XLSX) S6 Data. Results of search for publicly available genetic and geographic salamander data. (XLSX) S7 Data. Results of literature search for genetic lineages in recognized salamander species. (XLSX) Acknowledgments We thank Radford University Office of Undergraduate Research for supporting AEG through the Accelerated Research Opportunities and Research Rookies programs. We thank members of the Carstens lab for their comments on a draft of this manuscript. Author Contributions Conceptualization: Danielle J. Parsons, Bryan C. Carstens, Tara A. Pelletier. Data curation: Danielle J. Parsons, Abigail E. Green, Bryan C. Carstens, Tara A. Pelletier. Formal analysis: Danielle J. Parsons, Abigail E. Green, Bryan C. Carstens, Tara A. Pelletier. Funding acquisition: Bryan C. Carstens, Tara A. Pelletier. Methodology: Tara A. Pelletier. Project administration: Tara A. Pelletier. Resources: Bryan C. Carstens. Software: Danielle J. Parsons, Abigail E. Green, Tara A. Pelletier. Supervision: Tara A. Pelletier. Validation: Danielle J. Parsons. Visualization: Danielle J. Parsons, Tara A. Pelletier.