INTRODUCTION
Deyeuxia Clarion ex P. Beauv. (Poaceae) is one of the richest genera of native grass species in South America (Rúgolo de Agrasar, 2006) and the Andes are its center of diversification. The validity of this genus is still under discussion. According to a recent treatment, most of the South American species of Deyeuxia within the Koeleriinae clade B of the Aveninae J. Presl (Saarela et al., 2017) has been transferred to the genus Cinnagrostis Griseb. (Peterson et al., 2019). Their taxonomical decision was based on a molecular phylogeny that considers three plastid and one nuclear DNA regions though the type species C. polygama Griseb. had not been included (Saarela et al., 2017); also, they considered for this decision an unpublished study of Romaschenko et al. (Peterson et al., 2019). Furthermore, new combinations and the establishment of new genera were based on weakly supported clades (Saarela et al. 2017). On this basis, we decided to maintain the validity of the genus Deyeuxia.
Endemism is commonly found within Deyeuxia (Rúgolo de Agrasar, 2006). This genus is characteristic of the “Altoandina” phytogeographic province (Weberbauer, 1945; Cabrera, 1957). The weather is cold throughout the year, and snow and hail are the most common type of precipitation being less than 400 mm annually. In this harsh environment, xerophytic grasses and cushioned or crawling dicotyledons are prevalent. The vegetation is poor and Deyeuxia is mentioned as a characteristic grass of this province (Cabrera, 1978).
Deyeuxia velutina Nees & Meyen is a native species that inhabits high elevation areas of the Andes (Peru, Bolivia, Chile, and Northwestern Argentina; Rúgolo de Agrasar, 2006; Villavicencio & Rúgolo, 2014). The presence of this species in Peru and Bolivia is based on only one specimen of each country [the type specimen of D. scabriuscula Wedd., a synonym of D. velutina var. nardifolia (Weddell 441, P; Rúgolo de Agrasar, 2006) and Ruthsatz & Vogt 10507 (LPB), Villavicencio & Rúgolo, 2014, respectively]. Deyeuxia velutina has only been considered in a molecular phylogenetic study based on plastid and nuclear markers (Saarela et al., 2017). These authors found that D. velutina specimens were grouped in a clade with low support along with other taxa of Deyeuxia [D. breviaristata Wedd., D. brevifolia J. Presl, D. cabrerae (Parodi) Parodi, D. chrysophylla Phil., D. curvula Wedd., D. deserticola Phil., D. fiebrigii (Pilg.) Rúgolo, D. filifolia Wedd., D. malamalensis (Hack.) Parodi, D. trichodonta Wedd. and D. violacea Wedd.].
Based on Rúgolo de Agrasar (2006), two varieties of D. velutina were recognized: D. velutina var. velutina and D. velutina var. nardifolia (Griseb.) Rúgolo, which were differentiated by the following characters: the upper leaf sheath indumentum (dense pubescence vs. glabrous), floret length [3-4.5 (-5.5) vs. 5-5.5 mm], and leaf blade shape (flexuous, unusually straight vs. curved), respectively. The varieties grow in humid meadows but differ in their altitudinal occurrence along the Andean distribution range. Deyeuxia velutina var. nardifolia grows along central Andean mountains of Peru, Bolivia, Chile, and northwestern Argentina at 3200-5000 m a.s.l. (Rúgolo de Agrasar, 2006; 2009; Rúgolo, 2012; Villavicencio & Rúgolo, 2014). On the other hand, D. velutina var. velutina extends along the Andean region of Chile and Argentina at 2300-4200 m a.s.l. (Rúgolo de Agrasar, 2006; 2009; Rúgolo, 2012). Nevertheless, infraspecific identification of recent collections of D. velutina varieties along the overlapping distributional range have been problematical due to ambiguous diagnostic characters.
Wiens (2011) considered that biogeographic patterns arise from dispersal limits, and species do not extend beyond their environmental niche, including abiotic and biotic factors. Therefore, differences in environmental variables (precipitation, temperature, seasonality, altitude, among others) could drive speciation since climatic-niche divergence may conduct incipient species to occupy different habitats (Hua & Wiens, 2013). These variables combined with geographical occurrences are commonly used to delimit cryptic species or species complexes through species distribution modeling analysis (e.g. Wiens & Graham, 2005; Raxworthy et al., 2007; Rissler & Apodaca, 2007; Fernández et al., 2017). This analysis is useful for studying geographic distribution and habitat preferences, aimed at species conservation and taxa delimitation (Sites & Marshall, 2003; Wiens, 2007; Schlick-Steiner et al., 2010). The combination of morphological analyses with ecological approaches may explain the influence of abiotic factors on the geographical distribution as well as their effect on specific changes in phenotypic variation patterns (Vogel Ely et al., 2018).
Traditionally, taxonomy has relied on morphological characters for the circumscription and identification of taxa. Morphometrics based on similarities (or dissimilarities) among specimens have successfully helped in delimiting groups or taxa (e.g. Sassone et al., 2013; Nicola et al., 2014; Fernández et al., 2017; Viera Barreto et al., 2018; Moroni et al., 2019). This approach and comparative studies, when considering multiple methods and independent datasets, are increasingly relevant to understand the processes that guide speciation (Sites & Marshall, 2003; Wiens, 2007; Schlick-Steiner et al., 2010).
In this contribution, we explore the morphological variability of Deyeuxia velutina varieties and its association with their geographical distribution. Our work aims to improve their delimitation through multivariate analyses and to select reliable characters from univariate analyses to provide an updated key for taxa identification. A study of species distribution model is used to map and predict the potential distribution area of the two varieties.
MATERIALS AND METHODS
Sampling
Preliminary taxa identification and character selection followed the criteria of Rúgolo de Agrasar (2006). A total of 124 morphological characters (i.e., 14 quantitative and 110 qualitative; Appendix A) were defined, coded, and added to the data matrix for analyses. Morphological data were obtained from 52 specimens (i.e., 32 of D. velutina var. velutina and 20 of D. velutina var. nardifolia; Appendix B), housed at BAA, CONC, LIL and SI herbaria (acronyms following Thiers,
2020, continuously updated).
Micromorphological studies
Florets, as well as small fragments of leaf sheaths (middle zone), were placed in glass tubes with xylene and exposed to ultrasound for approximately one hour to eliminate superficial wax and impurities. The material was air-dried, mounted, and coated with a gold-palladium (40%-60%) alloy by a Thermo VGScientific (West Sussex, UK). Photomicrographs were obtained using a Phillips XL 30 (Einhoven, The Netherlands) Scanning Electron Microscope (SEM) at the Museo Bernardino Rivadavia (Buenos Aires, Argentina).
Multivariate analyses
An ordination method was performed on standardized distance matrices (Legendre & Legendre, 1998; Oja & Paal, 2004); Principal Coordinate Analysis (PCoA) was used to analyze the morphological dissimilarity among the specimens using the Manhattan distance coefficient. A double centered matrix was employed as input to calculate eigenvectors in NTSYSpc v.2.2 (Rohlf, 1998) which were used to project the specimens in a bidimensional space. We used the packages “corrplot” (Wei & Simko, 2017) and “ggplot2” (Ginestet, 2011) in the R environment (R Core Team, 2019) for graphical representations. Finally, using a set of quantitative and non-correlated characters, a Linear Discriminant Analysis (LDA) was used to test whether morphological variation is enough to discriminate between varieties. In order to select the uncorrelated quantitative characters, we performed a Pearson’s correlation coefficient analysis using variables with r<0.5.
The discrimination accuracy was evaluated using a leave-one-out cross-validation procedure to show the percentage of specimens correctly classified on the a priori groups. The LDA was performed in the R environment using the package “MASS” (Venables & Ripley, 2002).
Univariate analyses
The analysis of variance (ANOVA) and a non- parametric Kruskal-Wallis test were performed to evaluate significant differences of quantitative and qualitative variables, respectively, and to select diagnostic characters between groups as identified from the taxonomical relationships among OTUs as represented in the PCoA (Tables 1 and 2). Moreover, the Kruskal-Wallis test was also used for non-normally distributed and/or heteroscedastic characters (Nelder & Wedderburn, 1972). All statistical analyses were implemented using the package “stats” in the R environment.
Species Distribution Modeling (SDM)
Georeferences were recovered from herbaria labels of 104 specimens. After removing duplicates, the SDMs were performed with 68 specimens of D. velutina var. velutina and 26 of D. velutina var. nardifolia (Appendix B). We used the geographic coordinates [latitude (S) and longitude (W)] of the Argentinian and Chilean specimens which were obtained by the collectors using a handheld Geographical Positioning System (GPS) unit. The remaining point localities obtained from herbarium specimens were georeferenced using Google Earth 9.124.0.1 (https://earth.google.com/web). In order to avoid taxonomical and geographical errors, we carefully checked the identity of the specimens and their geographical data. We retrieved data for the 19 bioclimatic variables and altitude from WorldClim 1.4 (http://www.worldclim.com/version1), with 2.5 minutes (~5 km2) spatial resolution (Hijmans et al., 2005). We defined a rectangular area from 57° to 20° S and 72° to 52° W to calibrate our models explicitly excluding Bolivia and Peru due to scarce collections in these countries. The species distribution models were constructed using MaxEnt v.3.4.0 (Phillips et al., 2004) in the R environment through the package “dismo” (Hijmans et al., 2016). Maxent analyses were performed setting maximum iterations to 1000, and all other options were left as default (logistic output, convergence threshold of 0.00001, 10,000 background points, regularization multiplier of 1, default prevalence of 0.5 and autofeatures). We selected 70% data for training and the remaining 30% for testing. The initial model was run using the 19 bioclimatic variables and altitude. To avoid overfitting of the model we built Pearson’s correlation coefficient matrix and excluded variables (r>0.7). To decide which of the highly correlated variables should be left in the model, we used the Jackknife test of variable importance. This test provides information on the performance of each variable in the model (i.e., the importance of each variable on species distribution and the unique information it provides; Baldwin, 2009). The correlation matrix was constructed with the package “raster” (Hijmans, 2020). As a result, the final model included 6 variables for D. velutina var. velutina and 7 for D. velutina var. nardifolia, respectively (Table 3). The Area Under the Receiving Operator Characteristic curve (AUC) was used to test model’s goodness-of-fit. After calibration, we projected the results on South America.
Linear Regression analysis
A linear regression model was applied to measure the relationship of Axes I and II of PCoA (axes that explained most of the variability of the ordination analysis) with the geographical data (latitude and longitude). Consequently, the relationships between morphological variation and environmental changes were determined in both varieties of D. velutina. All georeferenced specimens considered in the morphometric analyses were included in the regression. The linear regression model analysis was implemented using the package “stats” in R environment (R Core Team, 2019).
RESULTS
Multivariate analyses
Sixty-one variable characters out of 124 were retained to calculate the Manhattan distance to perform Principal Coordinates (PCoA). Based on PCoA, we recognized two groups of operational taxonomic units (OTUs), corresponding to Group I and Group II (mostly in concordance with D. velutina var. nardifolia and D. velutina var. velutina, respectively). A few individuals first identified as D. velutina var. velutina (Group II) were mixed with D. velutina var. nardifolia (Group I). Both groups were mainly separated along the principal coordinate axis II (Fig. 1). Axes I and II explained 36% of the total variation.
Based on groups I and II, a Discriminant Analysis was performed using uncorrelated quantitative characters. Differences between a priori groups resulted significant (p < 0.001) on the discriminant function; the lemma length and floret width contributed most to the negative end of the linear discriminant associated with D. velutina var. nardifolia, while the inflorescence length contributed most to the positive end associated with D. velutina var. velutina. A total of 90% of the individuals were correctly classified as D. velutina var. nardifolia, while 91% were correctly classified as D. velutina var. velutina.
Univariate analyses
Groups as determined by multivariate analyses (PCoA and DA analyses) were tested for significant differences using ANOVA and Kruskal-Wallis tests. Univariate analyses showed differences in the height of culms, length of inflorescences, ligules and lemma, presence of hairy indumentum on the adaxial leaf blade and ligular zone, length ratio between glumes I and II, and blade apex shape (Tables 1 and 2).
An updated key to the identification of varieties is presented herein based on additional statistically significant characters; floret and leaf sheath indumentum are illustrated in Fig. 2. Quantitative characters that present overlapping ranges of variation between varieties have not been considered in the key.
1. Leaf sheath densely villous; ligular zone sericeous; ligule margin entire; leaf blade flexuous, exceptionally erect. Glumes shorter than the floret, upper glume exceptionally as long as the floret; rachilla extension ½ the length of the floret .............. D. velutina var. velutina
1. Leaf sheath glabrous; ligular zone glabrous; ligule margin denticulate, sericeous; leaf blade curved. Glumes as long as the floret; rachilla extension ¼-⅓ the length of the floret ................................ D. velutina var. nardifolia
Geographical distribution of varieties
Deyeuxia velutina inhabits high-elevation habitats along the Andes, with an altitudinal interval between 1300 to 5400 m a.s.l. The species have been reported to occur from Peru to the center of Argentina, reaching Neuquén, and our potential distribution models predict suitable environmental conditions between 8°-35° S and 68°-77° W (Fig. 3).
The species distribution models mostly reflected the known distribution of each variety. The average AUC test values of the SDM models were 0.96 and 0.97 for D. velutina var. velutina and D. velutina var. nardifolia, respectively. Although both potential distributions overlap, D. velutina var. nardifolia present a higher habitat suitability in the northern area (Fig. 3A) and D. velutina var. velutina in the southern zone of the distribution of the species (Fig. 3B). Even though the distribution of D. velutina var. nardifolia was modeled in Argentina and Chile, this variety showed a high habitat suitability in Bolivia and southern Peru. During a recent field trip, specimens of D. velutina var. nardifolia [Ferrero & Iacobucci 19 (BAA)] were also recorded for San Juan (Argentina), constituting the first report for this province and the southernmost limit of its geographical distribution. The models showed that the Altitude (alt) contributed to the potential distribution of both varieties. In addition, the Precipitation of Driest Quarter (BIO17) and the Precipitation of the Warmest Quarter (BIO18) contributed to the potential distribution of D. velutina var. nardifolia and D. velutina var. velutina, respectively.
Linear Regression analysis
Linear regression between PCoA axes and geographic coordinates was analyzed to test whether the geographic distribution is associated with morphological Groups I and II. Significant values for the regression model were found for Axis I vs Latitude (r2= 0.13, p≤0.001), and Longitude (r2= 0.12, p≤0.001), while Axis II showed a higher adjustment to the model only with Longitude (r2= 0.21, p≤0.001).
DISCUSSION
While molecular phylogenetic studies are valuable tools for understanding the evolutionary relationships between organisms, multivariate analysis has proved to be useful for the selection of taxonomic groups based on morphological variation. Moreover, univariate analyses contribute to the selection of diagnostic characters especially when the study involves several variables. Hence, integrative taxonomy is a multisource approach for combining information and disciplines for species delimitation (Schlick-Steiner et al., 2010). In this work, we used traditional morphometry and species distribution models to improve the accuracy of taxa identification based on an integrative taxonomy approach (Will et al., 2005). Our integrative approach has been effective in shedding light on the infraspecific classification of D. velutina varieties supporting the identity of both taxa as different varieties, like in other similar studies (e.g. Sassone et al., 2013; Nicola et al., 2014; Fernández et al., 2017; Viera Barreto et al., 2018; Moroni et al., 2019).
Morphometric studies were helpful in selecting reliable characters to separate these varieties. However, it is of particular interest that individuals with intermediate morphological characteristics have been found in areas where the varieties are sympatric (San Juan and Mendoza in Argentina and the Metropolitan Region of Chile; Fig. S1), indicating these areas as a possible hybridization zone. Taxonomic groups associated with both varieties, Groups I and II (Fig. 1), were differentiated in the PCoA and corroborated by LDA. Moreover, based on the univariate analyses, some of the characters considered by Rúgolo de Agrasar (2006) and a set of additional characters were considered for discriminating D. velutina var. velutina from D. velutina var. nardifolia: ratio of glumes to spikelet, rachilla length, ligular zone and ligular margin. These characters shed light on morphological differentiation and overlapping of known characters. Afterwards, we used the new extended key to assign misclassified specimens of D. velutina var. velutina as D. velutina var. nardifolia (Fig. 1).
Species Distribution Modeling
Deyeuxia velutina inhabits in areas with annual cumulative precipitation less than 400 mm and an average annual temperature ranging from 0° to 32°C. Both varieties partially overlap their potential distribution in areas with low habitat suitability. Notwithstanding, they are mainly segregated by the geographic latitude: D. velutina var. nardifolia occurs in the northern distribution area, while D. velutina var. velutina inhabits within the southern limit of the species’ range. The differences in the distribution of both varieties can indicate a divergence mechanism linked to new environmental conditions.
Several studies have indicated that the inclusion of elevation as a predictor variable of SDMs improves the quality of the prediction for high-elevation plant species (Austin, 2002; Körner, 2004; 2007; Oke & Thompson, 2015; Moroni et al., 2019). Oke & Thompson (2015) proposed that combining elevation, geographic and climatic data might predict more accurate species distribution models for high-elevation plants than those produced by climate data alone. In our results, altitude presents a high influence on both varieties potential distributions. “Precipitation of Driest Quarter” (BIO 17) and “Precipitation of the Warmest Quarter” (BIO 18) resulted the next most contributing variables in the potential distributions of D. velutina var. nardifolia and D. velutina var. velutina, respectively. The BIO 17 and BIO18 might be related to the presence of these varieties in humid areas and in the “Altoandina” phytogeographic province where the most common type of precipitation are snow and hail (Cabrera, 1978). On the other hand, the formation of circular and semicircular clumps, which have been described for high-Andean region as a consequence of ice and snow during the extreme winter (Hauman, 1918; Ruiz Leal, 1959; Rúgolo de Agrasar, 2006) would be associated to these bioclimatic variables. Moreover, Ruiz Leal (1959) assigned these structures as a response to altitude, humidity, low temperature, and snowfall. Finally, the Altitude contributed to the potential distribution of both varieties in agreement with previous observations of Rúgolo de Agrasar (2006) who reported D. velutina var. nardifolia in higher altitudes (3800-4700 m a.s.l.) than D. velutina var. velutina (2300-4200 m a.s.l.). Both varieties occur in the “the biome Montane Grasslands and Scrublands” defined by Olson et al. (2001, Fig. S1) and the “Altoandina” phytogeographic province (Cabrera & Willink, 1980). Alternatively, according to the main geographic units of the Andes defined by Weigend (2002), D. velutina var. nardifolia is distributed in the central Andes where the rainfall regime changes from summer-rainfall to winter-rainfall (Luebert & Pliscoff, 2006). Likewise, D. velutina var. velutina is distributed in the southern Andes where the elevation decreases from North to South (Pankhurst & Herve, 2007). In addition, our linear regression analyses supported that morphological differentiation that distinguishes the varieties of D. velutina also differentiate along the north-south and east-west gradients through the distribution and altitude of the Andes.