Morphometric analysis provides evidence for two traditionally defined species of the Tillandsia erubescens complex (Bromeliaceae)

1Doctorado en Ciencias en Biodiversidad y Conservación, Área Académica de Biología, Instituto de Ciencias Básicas e Ingeniería, Universidad Autónoma del Estado de Hidalgo, Mineral de la Reforma, México 2Centro de Investigaciones Biológicas, Instituto de Ciencias Básicas e Ingeniería, Universidad Autónoma del Estado de Hidalgo, Mineral de la Reforma, México 3Departamento de Botánica, Instituto de Biología, Universidad Nacional Autónoma de México, Mexico City, México *Corresponding author: hleoni@uaeh.edu.mx RESEARCH ARTICLE


INTRODUCTION
Species are the fundamental unit of biological classification; thus, carefully determining their identification, delimitation, and description is crucial in studies of systematic, evolutionary biology, biogeography, ecology, population genetics, and conservation biology (de Queiroz 2007;Freudenstein et al. 2017;Pinheiro et al. 2018).
The Tillandsia erubescens Schltdl. complex (Tillandsioideae, Bromeliaceae) currently includes eight morphologically similar species, some of them difficult to identify (Granados Mendoza et al. 2016). These species are endemic to Mexico and are distributed in temperate high mountain forests in a range of 1150-3200 m a.s.l., from the north to the south of Mexico (Granados Mendoza 2008).
Within the green corollas group, T. quaquaflorifera can easily be distinguished from the other two species by its deep violet and broadly ovate leaf sheaths, glaucous-greenish leaf blades, bright red scape bracts, and horizontally curved inflorescence. This species is restricted to the humid oak forest in the northern portion of the Sierra Madre del Sur (SMS), at elevations ranging from 700 to 3100 m (INEGI 2021).
However, morphological criteria used to separate T. erubescens from T. arroyoensis have been rather ambiguous (Weber 1983), presumably inconsistent (Granados Mendoza 2008), or not explicitly stated (Espejo-Serna et al. 2004). Tillandsia erubescens was first described by Schlechtendal in 1845. After this, two varieties were proposed by Weber (1983); namely, T. erubescens var. arroyoensis W.Weber & Ehlers and T. erubescens var. patentibracteata W.Weber & Ehlers. According to these authors, T. erubescens var. arroyoensis differs from the type variety in that all its structures are narrower, while T. erubescens var. patentibracteata shows a longer inflorescence with divergent bracts. Later, Espejo-Serna et al. (2004) changed the taxonomic rank of T. erubescens var. arroyoensis to the species rank (T. arroyoensis), and synomized T. erubescens var. patentibracteata with T. erubescens, albeit without stating the specific criteria supporting their decision.
According to Gardner (1986) and Granados Mendoza (2008), the three taxa share the presence of leaf blades densely covered by cinereous winged trichomes, simple and pendulous inflorescences, pink floral bracts, green petals, and green filaments and styles (exposed portion). Granados Mendoza (2008) suggested that T. erubescens var. erubescens differs from T. erubescens var. arroyoensis in its solitary growth or in colonies of up to four individuals (vs colonies of 25 to 40 individuals), threefold larger rosettes, ecarinate (vs carinate) and acute (vs rounded) sepals. However, a revision of herbarium and recently collected living specimens showed that the morphological features proposed by Granados Mendoza (2008) to differentiate T. erubescens var. erubescens from T. erubescens var. arroyoensis could be insufficient, since several individuals present overlapping characters (A. Martínez-García pers. obs.).
It has been reported that T. erubescens var. erubescens has a wide distribution, being located in pine oak forests and xerophilous scrub in the main mountain ranges of Mexico. That is, the Sierra Madre Occidental (SMOc), the Trans-Mexican Volcanic Belt (TMVB), the Sierra Madre del Sur (SMS), and the southern part of the Sierra Madre Oriental (SMOr) (Granados Mendoza 2008), at elevations from 1100 to 2900 m, at an average temperature of 17°C in a temperate subhumid climate (INEGI 2021). On another hand, T. erubescens var. patentibracteata is only documented in the central part of the SMOc (Weber 1983) in oak forests with climatic conditions similar to those recorded for T. erubescens var. erubescens. Tillandsia erubescens var. arroyoensis has been reported in pine oak forests only in the northern part of the SMOr (Granados Mendoza 2008) at elevations ranging from 1500 to 2900 m and an average temperature of 17°C in a temperate subhumid climate (INEGI 2021). Since T. erubescens var. erubescens and T. erubescens var. arroyoensis have not been documented in sympatry, their distribution could be determined by the presence of biogeographic barriers; however, this aspect has not been formally studied.
The aims of the present study were: 1) to test if the three taxa proposed by Schlechtendal (1845) and Weber (1983) (i.e. T. erubescens var. erubescens, T. erubescens var. patentibractetata, and T. erubescens var. arroyoensis) belong to separate taxonomic entities through linear morphometric analyses of both vegetative and reproductive characters; 2) to evaluate which taxonomic rank better reflects the morphological distance of the resulting taxonomic entities; 3) to address the geographic distribution patterns of the taxonomic entities recognized herein; and 4) to analyze the relationship among geographic distance, climatic variables, and elevation with morphological variables.

Plant material and data collection
A total of 110 specimens were studied, including new collections from the field, and type and regular herbarium specimens housed at the herbaria CHAP, CHAPA, ENCB, HAL, HGOM, JES, MEXU, UAMIZ, and WU (herbarium acronyms following Thiers continuously updated), which were previously named in herbaria as one of the following a priori taxa: T. arroyoensis, T. erubescens var. arroyoensis, T. erubescens, and T. erubescens var. patentibracteata, and were preliminary identified as such. Of the 110 specimens studied, only 44 were used in the multivariate analyses because they had complete and appropriate structures that could be measured. In some cases, floral characters (such as sepals, petals, and filaments) were not present and could therefore not be measured, and in many cases the samples only had fruits. In other cases, it was not possible to obtain permission to remove floral structures for rehydration and measurement. In the end, measurements were taken from 21 specimens corresponding to T. erubescens var. arroyoensis, 21 to T. erubescens var. erubescens, and two belonging to T. erubescens var. patentibracteata (which are the only existing records for this taxon) (supplementary file 1). Each specimen was treated as an individual taxonomic unit (OTU), conserving the previous identifications (as variety or species) as a reference.
Measurements were taken directly from dry specimens; large structures such as leaves and inflorescences were measured with a ruler, and smaller structures such as floral bracts were measured with a precision digital calliper (0.1 mm, Weston). Additional measurements were recorded from photographs taken by the first author with a CANON Powershot G10 camera and including a scale. In the case of the nomenclatural types, data were obtained from high resolution scans that included a scale available from the herbaria websites. All images were processed with TpsDig2 software (Rohlf 2018).

Linear morphometric analysis
Sixteen quantitative morphological characters were measured, six vegetative and six reproductive, as well as four proportions between two of the characters. An explanation of how these measurements were taken is provided in table 1. Due to the wide ranges of values, the data were standardized with the function (log + 1). Subsequently, a factor analysis (FA) was performed to reduce the dimensionality of the data. In order to find potential patterns of similarity between the OTUs and the consequent formation of groups (McGarigal et al. 2000), a cluster analysis (CA) was conducted based on a Euclidean distance matrix constructed with the selected variables from the FA. If the cophenetic value is close to 1, it is assumed that the clustering results give a good representation of the original distances.
To determine if there were statistically significant differences between the retrieved groups (taxa), a discriminant function analysis (DFA) was implemented. The morphological characters were used as discriminating variables and the grouping variable was the species. To measure the similarities or differences between the groups (species), Wilks' lambda (λ) was used, which is easy to interpret because its values vary from 1 to 0, where 1 means total similarity and 0 total difference. All the aforementioned analyses were performed using the software STATISTICA v.7.0 (StatSoft Inc. 2004).

Geographic distribution and climatic variables
Distribution records from 274 accessions corresponding to the three studied varieties were obtained from the Global Biodiversity Information Facility (GBIF 2020) and Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO 2020) databases, as well as from the collection locality named on the herbarium label. Those records without coordinates were manually georeferenced with the online version of Google Earth (Google Earth 2020). Once the result of the morphometric analysis was obtained, the information about the distribution documented in the online databases was compared to the distribution registered on the herbarium labels. Finally, the distribution projections were built from the taxonomic entities resulting from the morphometric analysis based on a comprehensive and well curated geographic database. The species distribution map was generated with QGIS v.3.4 (QGIS Development Team 2020), projecting the coordinates of each record on a layer of the biogeographic provinces of Mexico (Morrone et al. 2017). Additionally, we explored if, 1) climatic variables 2) elevation, and 3) geographic distance have an influence on morphological differences among taxonomic entities. For the above, we constructed a data matrix including the coordinates of the 44 localities, which corresponds to each of the measured specimens. Moreover, we used DIVA-GIS v.7.5.0 (DIVA-GIS Development Team 2021) to extract the values of the climatic variables of each site (temperature, precipitation) available in the WorldClim database (Fick & Hijmans 2017) at 30 arc-seconds resolution. Whenever available, elevation data were taken from the labels of the herbarium specimens and when not, it was taken by georeferencing the locality through Google Earth (Google Earth 2020). We performed a series of single Mantel tests using the R package vegan v.2.5-7 (Oksanen et al. 2020; R Core Team 2020) with significance determined using 9999 permutations, testing the correlation between various dissimilarity matrices: 1) climatic variables temperature and precipitation (Euclidean distance), 2) elevation (Euclidean distance), 3) geographic distance (Haversine distance), and 4) resultant scores from the morphological FA analysis (Euclidean distance).

Groups recovered by linear morphometrics and taxonomic rank
Of the 16 characters analyzed in the FA, there were 13 characters that contributed the most to the differentiation of the groups: total size of the plant (TS), width of the rosette (WR), all characters of the leaf sheath and blade (WS, LS,  2) also reveals the existence of two groups; the first one includes specimens a priori determined as T. erubescens var. arroyoensis, and the second includes specimens identified a priori as T. erubescens var. erubescens and T. erubescens var. patentibracteata. Although the two specimens of T. erubescens var. patentibracteata are closely positioned in the morphospace, they are imbedded within the T. erubescens var. erubescens specimen morphospace. The dendrogram derived from CA confirms the two main groups as well ( fig. 3), with T. erubescens var. patentibracteata and T. erubescens var. erubescens in the same clade, and T. erubescens var. arroyoensis in a separate clade. The high value of the cophenetic correlation of this analysis (0.820) shows that the resolution of the CA (i.e. dendrogram) faithfully represents the structure of the original dataset.
Values close to 1 (as in this cluster analysis) means a good representation of the original distances. The DFA showed that the two groups (hereafter assumed to be species) are significantly different (F 13,30 = 32.128, p < 0.0000, Wilks' lambda = 0.06701, eigenvalue = 13.92204, Canonical R = 0.965911, χ 2 = 95.95079). In 100% of the   cases, both species were correctly classified based on values of the statistic and supported by morphometry.

Geographic distribution and climate variables
Non-overlapping distribution areas were found for the two resulting taxonomic entities. Tillandsia erubescens (including T. erubescens var. erubescens and T. erubescens var. patentibracteata) exhibits a very wide geographical distribution throughout the SMOc province in Sonora, Chihuahua, Sinaloa, Durango, Zacatecas, Aguascalientes, Nayarit, and Jalisco; TMVB in Colima, Michoacán, Puebla, Tlaxcala, Estado de México, and Mexico City; the Balsas Basin (BB) in Guerrero and Oaxaca; the SMS in Jalisco; and the southern part of the SMOr in Hidalgo, Querétaro, northern Veracruz, Puebla, and southern San Luis Potosí ( fig.  4). In contrast, T. arroyoensis is restricted to the northern portion of the SMOr province in Coahuila, Nuevo León, Tamaulipas, and northern San Luis Potosí ( fig. 4).
According to the results of the Mantel test (table 3, fig.  5), precipitation (r = 0.2546, p = 0.0001) and geographic distance (r = 0.2256, p = 0.0024) are correlated with the morphological differences between the two species. On the other hand, temperature (r = 0.04744, p = 0.2208) and elevation (r = 0.03197, p = 0.2621) have no correlation with the morphology of the species.

Morphometric delimitation, taxonomic rank, and definition of the three Tillandsia erubescens varieties
The statistical analyses allowed two groups to be clearly distinguished. One of them includes all specimens a priori identified as T. erubescens var. arroyoensis and the other contains the specimens identified as T. erubescens var. erubescens and T. erubescens var. patentibracteata. The characters (in order of contribution) that support this grouping were number of flowers (NF), plant total size (TS), those related to the width and length of the leaf sheath and leaf blade (WS, LS, WLB, LLB), and those related to the inflorescence (LPII, WI, LI). These characters are smaller and narrower in the first group (table 4), which agrees with Weber's (1983) original description of T. erubescens var. arroyoensis, mentioning that all its structures are narrower than in T. erubescens var. erubescens. We found that the number of flowers per inflorescence is fundamental for differentiating between the two varieties and is herein proposed as a diagnostic character. According to our observations, T. erubescens var. arroyoensis does not exceed two flowers per inflorescence, whereas T. erubescens var. erubescens ranges from five to 15. However, the characters associated with the leaf sheath and blade (WS, LS, WLB, and  (table 2).
The morphometric analysis did not support Weber's (1983) proposal of recognizing T. erubescens var.
patenitibracteata as a separate taxonomic entity, since in the statistical analysis their measurements always remain within the range reported for T. erubescens var. erubescens (table 4). Therefore, from a morphometric point of view, there is no basis for recognizing the varieties of the second group as different taxa, however, according to the presented data no specimens are separated from the variety type. Nevertheless, it is important to note that only two specimens identified a priori as T. erubescens var. patentibracteata (one of them the type), could be included in this study, which contrasts with the large number of specimens that were assigned to T. erubescens var. erubescens. Consequently, no morphological character was found to support the separation of these varieties as independent taxonomic entities, as  Table 3 -Comparative values of the characters of T. erubescens var. arroyoensis, T. erubescens var. erubescens, and T. erubescens var. patentibracteata. All measurements are given in centimetres.

Morphology vs geographic distance (km)
Figure 5 -Scatterplot of the Mantel test for morphology vs climatic variables (precipitation and temperature), elevation, and geographic distance.
these two specimens are not separated from T. erubescens var. erubescens by the analysis. Therefore, we consider that T. erubescens var. patentibracteata could be considered a synonym of the type variety. However, this proposal should be further investigated by including more specimens, if they are found, in future studies. The relationship between morphology and floral visitors is well-known (Benzing 1998), and although there are no specific data on the pollination syndrome of these three taxa, previous studies report that species with exserted stamens are mostly pollinated by hummingbirds (Carranza-Quiceno & Estévez-Varón 2008;Barfuss et al. 2016). Observations of hummingbirds visiting the flowers of T. erubescens var. erubescens were made, both in the field and in cultivation (Claudia Hornung-Leoni pers. obs.).

Geographic distribution and climatic variables
The two taxonomic entities recovered herein (i.e. T. erubescens and T. arroyoensis) showed non-overlapping geographic ranges. T. erubescens is widely distributed in SMOc, TMVB, BB, SMS, and in the southern part of the SMOr, while T. arroyoensis is restricted to the northern part of the SMOr province. The wide distribution of T. erubescens could be due to its epiphytic and epilithic lifestyle, which could have allowed it to colonize a wider variety of niches available in different types of vegetation, such as pine oak and oak forests, as well as xerophilous scrub. The exclusively epiphytic lifestyle of T. arroyoensis, together with its presence only in pine oak and oak forests, could be associated with its narrower geographic distribution (Granados Mendoza 2008).  According to the Mantel test, the precipitation distance matrix has a strong correlation with the dissimilarity matrix of morphology. This is, as the samples become more different in terms of precipitation, they also become more different in terms of morphology. It has been reported that within the microenvironmental conditions inherent to epiphytes, including the genus Tillandsia, the availability of water (which is taken in pulses during rainfall events, including mist and dew), is the environmental factor that most influences its distribution (Benzing 1998). This can cause displacements or even the disappearance of species in certain areas (Benzing 1998;Cach-Pérez et al. 2014). This could also explain the wide distribution of T. erubescens, which is distributed along the sites with the highest levels of precipitation. The morphology dissimilarity matrix also had a strong significant correlation with the geographic separation of the samples, that is, the more the samples were physically separated, the more the morphology of the specimens differed. This may be due to the fact that each group, even when it is present in similar climatic conditions, is immersed in a different niche, which could eventually favour the morphological differentiation between both species. In contrast, the species dissimilarity matrix did not have a significant correlation with the dissimilarity matrix of temperature or elevation. This suggests that, at least within the range in which these species are recorded, as temperature or elevation changes, the morphology does not necessarily change as well. This is consistent with the values obtained for temperature and elevation, which are similar in the species' distribution zones (17°C and from 1100 to 2900 m).
The remarkable environmental heterogeneity of the mountainous tropical forests of North and Central America, together with the complex climate history of these regions during the late Miocene and Pliocene, which included an increment in seasonality and aridity induced by a global cooling trend (Antonelli & Sanmartín 2011), could have promoted high rates of diversification of the subgenus Tillandsia, resulting in the current high richness of species present in this region (Granados Mendoza et al. 2017).
Previous studies have recovered T. erubescens and T. arroyoensis as sister species (Granados Mendoza et al. 2016). It could be hypothesized that their common ancestor was widely distributed and that climate changes during the Pleistocene could produce geographic isolation as a result of expansion, fragmentation, and divergence between its populations. This could eventually give rise to what we currently recognize as T. erubescens and T. arroyoensis. These kinds of speciation events have been reported to have occurred in different types of vegetation currently distributed in Mexico and North America (Graham 1999;Metcalfe et al. 2000), such as hummingbirds and flying squirrels inhabiting pine oak forest and mountain cloud forest (Kerhoulas & Arbogast 2010;Cavender-Bares et al. 2011;Ornelas et al. 2013Ornelas et al. , 2015. It should, however, be noted that this speciation hypothesis needs to be formally evaluated using, among other evidence, a dated phylogeny. On the other hand, it has been proposed that the SMOr province is not a natural biogeographic unit and it could be divided into southern and northern portions. The southern part comprises crasicaule, microphyllous, and rosetophyllous scrubland vegetation types (Salinas-Rodríguez et al. 2018) and is represented by several endemic species of mammals, birds, and cacti (León-Paniagua et al. 2004;Navarro et al. 2004;Del Conde Juárez et al. 2009). The northern part contains the highest mountains, with temperate forests dominated by oak, pine, and chaparral communities (Salinas-Rodríguez et al. 2018) and it is characterized by several taxa of coleopterans, mammals, reptiles, amphibians, and gymnosperms (Luna-Vega et al. 1999;Márquez & Morrone 2004;Medina-Romero 2009). The presence of T. arroyoensis in the north of the SMOc province supports the idea of this biogeographic division.
Although it is thought that the Bromeliaceae as a whole is not yet exposed to anthropic threats, some species have been systematically depleted since they are used indiscriminately by the human communities that inhabit the areas in which they are distributed (Espejo-Serna & López-Ferrari 2018). The correct delimitation of the species will allow a better estimation of the diversity of this group, its distribution area, and the problems that it is facing, as well as will promote studies or programs that address aspects of its conservation and sustainable use. Despite the fact that the advance in the floristic knowledge of this important family has been constant in recent years, more studies are required in groups as Tillandsia in which complexes of species of difficult delimitation and recognition are presented (Espejo-Serna & López-Ferrari 2018), indicating a morphological, ecological diversity, genetics, among others, that have not been yet sufficiently explored.

CONCLUSION
This study provides morphometric evidence for the recognition of two traditionally defined species within the green corollas clade of the T. erubescens group based on a set of quantitative characters. Through the use of linear morphometry, it was possible to recognize two species: T. arroyoensis and T. erubescens, from the three varieties originally proposed. The characters that enable them to be differentiated are the total size and the sheath and leaf blade (AVF, LVF, ALF, LLF) and the number of flowers. T. erubescens var. erubescens and T. erubescens var. patentibracteata are considered synonymous because no significant differences were found that would allow them to be separated. We conclude that precipitation and geographic distance play an important role in the morphological divergence between species. The distribution of T. arroyoensis and T. erubescens shows a geographic delimitation that divides the species and suggests that they are separated by geographic barriers; however, more studies are needed to explain the cause of this separation.