We reviewed the most relevant Salvia taxonomic papers from Michoacán [Espejo-Serna and Ramamoorthy (1993), Bedolla-García et al. (2011), Cornejo-Tenorio and Ibarra-Manríquez (2011), Iltis et al. (2012), Martínez-Gordillo et al. (2013), González-Gallegos (2014), González-Gallegos et al. (2016), Lara-Cabrera et al. (2016), Martínez-Gordillo et al. (2017) and Zamudio and Bedolla-García (2018)] and obtained a summarised species list for this study.

Given the taxonomic complexities of the group, with frequent synonymy changes, we only considered the species that were reported in at least two of those studies and that could be corroborated with herbarium specimens. Discrepancies in the Salvia list include the number of species, 78 in three papers (Cornejo-Tenorio and Ibarra-Manríquez 2011, Lara-Cabrera et al. 2016, Martínez-Gordillo et al. 2017), although this number is slightly overestimated; for example, Salvia arbuscula Fernald is recognised by Cornejo-Tenorio and Ibarra-Manríquez (2011), but considered as a synonym to S. iodantha Fernald by Lara-Cabrera et al. (2016). Martínez-Gordillo et al. (2017) do not even include the species for the State. Finally, we updated the list including a recently-described species S. madrigalii Zamudio and Bedolla-García (2018).

Salvia occurrences information for Michoacán was obtained from databases of the National Biodiversity Information System of Mexico (SNIB-REMIB) of the Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO), the digital repository of the National Herbarium of Mexico (MEXU-UNIBIO) of the Instituto de Biología, Universidad Nacional Autónoma de México (UNAM) and through the review of taxonomic studies. All specimens deposited in the National Herbarium (MEXU) were reviewed, verifying species identification (see Suppl. material 1).

The database was corrected following the recommendations by Chapman (2005) and Castillo et al. (2014) with the following procedure: 1) correct identification of the specimens as verified by specialists, 2) eliminate synonym and duplicates, 3) geo-referencing the localities lacking this information in the specimen label, in Google Earth (https://www.google.com/earth/), using locality name and description and finally 4) eliminate correlated records through pattern analysis.

Pattern analysis was applied to the data for all species using the Ilwis v.3.4 programme (http://52north.org/ilwis). This analysis allowed us to estimate the distance at which the collecting points are not correlated (Cruz-Cárdenas et al. 2014a). The distance obtained from the pattern analysis was used to filter the uncorrelated data through the R spThin package (Aiello-Lammens et al. 2015, R Development Core Team 2016), integrated in the Data module Process Occurrence in Wallace v.1.0.6 (Kass et al. 2018). This analysis reduced sampling bias by eliminating records per geographic distance, which significantly influences the models performance (Boria et al. 2014).

In addition to the review of the specimens in the herbarium, specimens cited in taxonomic works of the study group were also considered. For the species and specimens with taxonomic circumscription issues, the specimens were individually evaluated to determine whether or not to include them in the analyses. The outliers identified in the pattern analysis were reviewed to confirm the identity from a genus specialist; the specimens that were not approved by the specialist, were excluded from the database.

The polygon for the state of Michoacán was divided into 118 grid-squares of 15' latitude and 15' longitude; a finer scale grid would have resulted in no species being represented. Thus 15' x 15' grid was a compromise between the distribution, collection effort and the number of grids. Geographic information was analysed in ArcGIS 10.2 (ESRI 2013) to produce a known distribution map and intersecting the grids with the number of species recorded in each grid. Grids were subsequently grouped per species number.

Fifty eight environmental variables were considered, with a resolution of 30 arc seconds, a pixel size of about 1 km² (Table 2); 26 of them were climatic (Hijmans et al. 2005), nine topographic, nine edaphic and 14 included remote sensing data (Cruz-Cárdenas et al. 2014b). However, it is known that the use of a large number of variables influences the final results of the models, making it necessary to select only those that are most significant (Peterson and Nakazawa 2008, Varela et al. 2014). To reduce the total number of variables included in models and to minimise multicollinearity amongst the 58 variables, their correlation values (Pearson) were evaluated using NicheToolsBox (Osorio-Olvera et al. 2018); uncorrelated variables were included in the analysis and correlated variables (r > 0.85) were excluded from the analysis.

Species distribution result from several factors; amongst the most important ones are environmental variables (A), biotic component (B) and the set of sites that has been accessible to the species (M). Soberón and Peterson (2005) illustrated the interaction of these factors in a diagram called BAM. We used the WWF Ecoregions following Olson et al. (2001) to delimit the accessible area for all species (M in the BAM diagram, Soberón and Peterson 2005), considering only the ecoregions with occurrence records of Salvia in Michoacán. Selected ecoregions were then cut, using the State's polygon as a limit.

Evaluating models using spatially independent data improves model configurations and balances their complexity. There are several methods for identifying optimal model configuration (Muscarella et al. 2014, Hijmans et al. 2017). Here, we used the package ENMeval v.0.1.2 in R (Muscarella et al. 2014) to determine the best performance models, run with multiple regularisation (MR) values, from 1 to 4 (in a 1 increment), with four different feature class combinations (L, LQ, H, LQH, where L = Linear, Q = Quadratic, H = Hinge). ENMeval provides multiple evaluations that allow the identification of the optimal model configuration (Muscarella et al. 2014).

Maxent 3.4.1 software (Phillips et al. 2017) was used to run the ENM. Maxent has gained popularity for ENM of many plant and animal species due to its simplicity in configuration. Models were constructed by changing the regularisation values and the feature class, following the results from ENMeval analysis. The default options of clamping and extrapolate were omitted, in order to avoid extrapolation of the extreme values of variables (Elith et al. 2011). Software was set to use 75% of the data to run the analysis and 25% to validate the model (Phillips et al. 2006), only for species with more than five records, the minimum number required in Maxent (Hernandez et al. 2006). For a few records species (5-20), we decided to construct the models using the basic configuration of Maxen (Phillips et al. 2017) and not the feature classes and regularisation multiplier, since these generally result in low Area under the Curve (AUC) values and overprediction. Species with less than five records were considered for generating the map of known richness and are considered amongst the rare species.

Results from Maxent were then processed in ArcGIS 10.2 (ESRI 2013) to obtain binary maps, using the maximum training sensitivity plus specificity as the cutoff threshold (Liu et al. 2016). These binary maps were added with the Algebra map tool to obtain a richness map of Salvia species. In addition, the known richness in the State was obtained from collection records in digitised databases, compared with the potential richness obtained through the models. Mexdem (digital elevation model) altitude layer in ArcGIS was used to extract the altitude values of each record of Salvia species.

We partitioned the data into training and testing groups (k-fold cross validation), allowing us to evaluate models’ performance (Peterson et al. 2011). We chose to use "block (k-4)", one of the options implemented in ENMeval (Muscarella et al. 2014), which divides the presence data according to its longitude and latitude. This method is one of the most recommended, since it provides better spatial independence of the data (Radosavljevic and Anderson 2014, Fourcade et al. 2018).

The models were evaluated using the AUC, a value that is part of Maxent's results. AUC values range from 0 to 1, where values close to 1 indicate models with perfect discrimination ability and values less than 0.5 indicate that the model is no better than a randomly-generated model (Peterson et al. 2011, Fourcade et al. 2018). Individual models were also evaluated by the binomial test analysis and partial-ROC (receiver operating characteristic), considering a significance value < 0.05. This allowed the evaluation of the models’ ability to predict the largest possible number of independent points amongst those that were not used to create the model. The tests were carried out using the R statistical package and NicheToolBox (Osorio-Olvera et al. 2018).

To assess the conservation status of the genus Salvia, habitat loss and the area currently within a Protected Area (PA) were analysed along the potential distribution of the genus in the State of Michoacán. To analyse the habitat loss of the genus, a land use map of Series I and VI of the National Institute of Geography and Informatics (INEGI 1997, INEGI 2017) was used. Maps were then reclassified into two categories: vegetated and non-vegetated (Staude et al. 2019); (1) Specifically, those areas that represent some forest cover within the State polygon were selected (e.g. oak forest, coniferous forest, evergreen forest, deciduous forest, sub-deciduous forest and xerophytic scrub). (2) To estimate the degree of protection of Salvia, the map of potential distribution was overlapped with the map of protected natural areas of Mexico (CONANP 2010).

A list of 66 species was obtained for the Mexican State of Michoacán (Suppl. material 2) from bibliographic review of Salvia. Fifty two species are reported as endemic to Mexico and three species are restricted to the State (S. madrigalii, S. subobscura and S. synodonta). A complete dataset record of 3,093 was later reduced to 1,925 records after data cleaning. Of these records, 1,836 correspond to the 42 Salvia species with pattern analysis and, finally, 404 records were used to build the ENMs.

A total of 26 uncorrelated variables (of the initial 58) were considered to run the ENMs for each Salvia species with sufficient (> 5) records, including eight climatic variables, seven topographic, eight edaphic and three from remote sensing data (Table 2). The variable of highest contribution value (42% of the models) was the Annual mean temperature (bio01) based on the Jackknife test provided by Maxent, with contribution percentages between 36 and 77.5%. The second most important variable was the Magnesium content (mexmg), important in 14% of the Salvia models, with contribution percentages from 19 to 69% and the least informative variable was the Precipitation of the coldest quarter (bio19) with 0.5%.

We obtained 42 ENMs of Salvia species in Michoacán (Suppl. material 3) out of 66 species, based on 404 spatially-uncorrelated records (Table 3); 10 of these models are considered exploratory having less than the five records required to be validated statistically. In addition, only three models were run using the regularisation values and the feature class obtained from the evaluation with ENMeval (Suppl. material 4); all the other models were obtained using the default parameters as specified in Maxent.

All models (made by default setting and configured) showed good performance, AUC values above 0.84 were considered good, 19% (AUC > 0.8) and 81% excellent (AUC > 0.9, Table 3). The validation of the models with the binomial test indicated that 42% models had statistically-significant prediction capacity (0.53-0.93, p-value = 0.03-3.46 x 10^-14 respectively), while 58% of the models were not statistically significant. In contrast, the Partial-ROC tests were statistically significant (p < 0.05) for all species, with values between 1.13 and 1.99 (mean 1.58, Table 3).

The geographic area obtained with the assembly of models for the Salvia species has been profoundly affected by land use change. It has been estimated that 22% of the area has been lost in the last 20 years (1997 - 2017); TMVB ecoregion is the most severely affected, with a reduction of about 43% loss of its primary vegetation. The least affected primary vegetation loss ecoregions are JDF and the SMS (less than 10%).

Ten Protected areas in Michoacán harbour potential distribution areas of Salvia species. Overall, these PAs include 3.3% of the potential distribution predicted. Eight of these PAs are located in the TMVB ecoregion; unfortunately, the SMS ecoregion does not register any PA. A total of 55% Salvia species in the State are represented in the PAs; of these, 25 species are distributed in the Monarch Butterfly Reserve (the Monarch Butterfly Hibernation Protection Zone). Only 6% of the rare species include populations within the State's PAs.

The Salvia’s database of Michoacán State includes more than 3,000 records, reduced to 1,925 unique records after eliminating duplicates. Although the database includes a significant number of records for the State, only 39.6% are cited in the referred floristic-taxonomic studies of relevance to the State. The two projects, focused on Salvia in Michoacán (Cornejo-Tenorio and Ibarra-Manríquez 2011, Lara-Cabrera et al. 2016), combine up to 701 specimens representing 36.4% of the total database; 685 specimens were collected before 2010 and are currently already included in online databases. Both studies cited a total of 217 reviewed specimens, representing 44% and 31.2%, respectively of those here reviewed. Although they do not explicitly indicate whether the number of specimens cited represent a total sample or a sample of the reviewed specimens, only 39.6% of the total collections for Salvia are cited.

Undoubtedly, 1,925 records representing 66 of the 69 species occurring in the State (Martínez-Gordillo et al. 2017), are a good sample of the State's generic diversity, with broad habitat diversity where the genus thrives. Many additional records here documented were reviewed by the group's specialists, although they were not cited in their floristic-taxonomic works. The 1,925 unique records included here are a summary of more than 3,000 specimens kept in about 20 different herbaria, both in Mexico and abroad (Table 1). The National Herbarium (MEXU) at the Instituto de Biología, Universidad Nacional Autónoma de México stands out, housing 64.5% of the total records (1,251); 1,159 of the 1,836 records used for the elaboration and validation of ENMs were at MEXU, constituting the main supply of information and headquarters of the study carried out here.

Our results indicate that Salvia species in Michoacán preferentially occur in temperate or seasonally-dry forests, the predominating biomes in the State (Gopar-Merino and Velázquez 2016); that is why more than 50% of the surface of the State is suitable for the species. Our results agree with those reported by Lara-Cabrera et al. (2016), which document higher species richness in the TMVB, where coniferous, oak and humid mountain forests abound. Similarly, Cornejo-Tenorio and Ibarra-Manríquez (2011) report that 40 Salvia species are distributed in the pine-oak and oak forests. Thus, our ENM are consistent with the previously reported environmental preferences for the genus.

The majority of Salvia species act as generalists to environmental conditions throughout their distribution in Michoacán. Only temperature had the highest variable contribution amongst all biomes; temperature is frequently documented as the dominant abiotic driver in determining plant distributions (Körner and Hiltbrunner 2018, Scherrer and Guisan 2019). Mean temperature is associated with plant growth, acting mainly in a gradual manner (Körner et al. 2016), determining whether a plant species will be present (Guisan et al. 1998), particularly in the mountainous landscape, such as is characteristic in the TMVB (Cornejo-Tenorio and Ibarra-Manríquez 2011, Lara-Cabrera et al. 2016). Oddly, precipitation was not a significant variable; a correlation between precipitation and species richness at regional scale in these environments has been well documented, the greater the precipitation, the greater the species richness (Pau et al. 2012). However, precipitation's contribution here is low, probably due to scarce records from the Balsas Dry Forests and JDF (Fig. 2), where rainfall is below 1,000 mm.

Rare species (< 3 records), although not modelled, show a preference for similar environmental conditions as compared to modelled species or widely-distributed species (> 3 records). For example, the mean temperature of rare species (18°C) is slightly higher than that of widely-distributed species (15.7°C), while in precipitation, the mean values are very similar, 1,057 mm for rare species and 1,058 mm for widely-distributed species. Furthermore, several state endemic Salvia species require very specific conditions; S. madrigalii, S. synodonta and S. cyanantha, are known only from two sites with practically the same climatic characteristics. The specificity of their habitat, as a result of their restricted distribution, makes these species particularly susceptible to extinction (Isik 2011). Despite having information on the optimal environmental conditions of this species group, very little is known about the biological or physiological aspects that regulate it.

Environmental suitability is indicated for the genus along the SMS (Fig. 2); however, further fieldwork for the region could add important floristic information (Cornejo-Tenorio and Ibarra-Manríquez 2011, Lara-Cabrera et al. 2016), as far as the optimal environmental conditions where they thrive. The above is a particular way for endemics or rare species that would allow us to confirm or improve the distribution patterns reported in this study.

Cornejo-Tenorio and Ibarra-Manríquez (2011) mention altitude as an important variable in the distribution from Salvia. The inclusion of this variable in the potential geographic range of Salvia in Michoacán confirms its importance in explaining the species distribution. The genus inhabits altitudes from 200 to 3,600 m, as reported by Lara-Cabrera et al. (2016); however, species richness is not evently distributed throughout this altitudinal range, with middle altitudes (2,000-2,500 m) as the richest in species (52 species) and decreasing richness both at lower and higher altitudes. These data coincide with those published by Cornejo-Tenorio and Ibarra-Manríquez (2011), who mention that the highest richness was found above 2,000 m altitude (35 species).

Although the large amount of information in the databases would seem sufficient for ENMs and thus inferring many species distribution, this was not always the case. When performing statistical analyses to develop more robust and precise models, such as eliminating spatial autocorrelation of points to reduce sampling bias or to optimise the specific configuration for each model adjusted for each species data (Boria et al. 2014, Radosavljevic and Anderson 2014, Warren et al. 2014), it became clear that we did not have enough independent records. In the case of Salvia in Michoacán, despite several works specifically focused on these species (Cornejo-Tenorio and Ibarra-Manríquez 2011, Lara-Cabrera et al. 2016, Zamudio and Bedolla-García 2018), we still did not have sufficient data for most species to develop ENMs under specific configurations. We found that 78% of the purged data for Salvia are spatially correlated, which is likely the result of sampling bias to areas of easy access or of particular interest, resulting in several records from a single site for the same species. Due to this strong point correlation, we were able to obtain only three configured models (models whose rigour was increased by adjusting their configuration or balancing their complexity). We corroborate the guidelines by other authors, at least 20 records are necessary to obtain and reach the predictive capacity for models (Stockwell and Peterson 2002, Mateo et al. 2010, van Proosdij et al. 2016). In our case, using few records to obtain regularisation values, we obtained models with low predictive value, resulting in high rates of model overprediction and AUC values below 0.5. For this reason, we too recommended that the tests to obtain regularisation values and make configured models be used only when more than 20 records are available; the larger the number of records, the better the model performance (Stockwell and Peterson 2002) as shown in Salvia.

The AUC has been criticised as a tool to evaluate only the performance of presence models, because of its dependence on prevalence; therefore, it is not considered a precise performance index (Lobo et al. 2008, van Proosdij et al. 2016). For this reason, other performance evaluation metrics, such as Partial ROC for continuous data (Peterson and Nakazawa 2008) and the binomial test for binary models (Anderson et al. 2002), were implemented. We included these test results as additional evaluation measures (Table 3). The number of non-significant binomial tests (58% of the models) correspond to the number of species with less than five records. Blair et al. (2013) obtained similar results when using ten records; below this sample size, their binomial tests were not significant (p > 0.05). We recommend reaching a minimum number of 20 independent points of occurrence to construct models with better predictive power that pass the test of sensitivity to small samples.

High deforestation rates in the north-central region of Michoacán, is the main cause of habitat loss. The most vulnerable vegetation types are temperate forests, constantly being replaced by avocado orchards (Mas et al. 2017) amongst other threats, thus decreasing the original distribution of many species, as we found here for Salvia.