Plant material

Samples of foliage and mature female strobili (cones) from the upper part of the crown were collected from 153 trees from various localities (Figure 1). We sampled 66 individuals of the typical variety P. pseudostrobus var. pseudostrobus (Chiapas, Mexico, Oaxaca, Puebla, Tlaxcala, and Veracruz), 60 of oaxacana (Chiapas, Oaxaca, Puebla, and Veracruz), 17 of coatepecensis (Veracruz) and 10 of apulcensis (Hidalgo).

Figure 1 Geographical location of trees sampled for each variety of Pinus pseudostrobus. 

Morphological and anatomical traits evaluated

Five healthy and developed fascicles (i. e., formed the previous year) and five mature, well-developed cones without apparent defects or malformations were selected from each tree. Information on 20 needle and 12 cone traits was collected using the methodology described by ^{Stead (1983}).

In foliage, the number of needles per fascicle (NN), sheath length (SL, cm), length (NL, cm) and thickness (NT, mm) of the longest needle in the fascicle, needle shape coefficient (NSc = NT/NL, mm∙cm^-1) and sheath/needle ratio (SR = SL/NL) were measured. Also, number of 'teeth' (NT), lines with stomata (LS), number of stomata per line (NS), and total stomata (TS = LS x NS) were counted in the mid-section of the longest needle, 5 mm from the dorsal side. In cross sections of the middle section of the same needle, the total number of resin canals (RC_T) with external (RC_E), internal (RC_I), middle (RC_M) and septal (RC_S) position was counted, and the proportion of resin canals with middle position was calculated (PRC_M = RC_M/RC_T). Similarly, the number of intrusions of the hypodermis into the endoderm, number of hypodermal cells touching the endoderm (CH_E), number of endodermal cells on the dorsal surface of the vascular bundle (EC) and the RC_H ratio (RC_H=CH_E/EC) were counted.

In closed cones, the length of the peduncle (PL, mm), length (CL_C, mm) and maximum diameter (CD_C, mm) of the cone were measured and their shape coefficient was calculated (CS_C = CL_C/CD_C); they were then dried at room temperature to open the scales and extract the seeds. Once the cones were opened, the maximum diameter (CD_O) was measured, and their degree of openness was calculated (OD_C = CD_O/CD_C). Three fertile scales were separated from the central part of each cone and the width (W_A, mm), depth (D_A, mm) and height (H_A, mm) of the apophysis were measured and the ratios between them were calculated (R_WD = W_A/D_A, R_WH = W_A/H_A and R_DH = D_A/H_A).

Statistical analysis

Univariate and multivariate analysis of variance was performed to determine the magnitude of trait variation between and within the OTUs, and to evaluate their ability to discriminate the putative varieties of P. pseudostrobus. In the univariate analysis, the MIXED procedure of SAS® version 9.4 (Statistical Analysis Software Inc., 2015) was used to estimate mean values for each putative variety and, when there were significant differences (α = 0.05), comparisons between varieties were made with Tukey's test. The discrete variables were transformed with the square root function to meet the criteria of normality and homoscedasticity; however, the results were similar to those of the original variables, so it was decided to present the analyses of the non-transformed values. In the multivariate analysis, carried out with the GLM procedure of SAS, the statistical significance of the variety effect was determined and the oaxacana and apulcensis varieties were compared by a contrast, using needle and cone traits separately.

To evaluate the effect of considering the oaxacana and apulcensis varieties as different UTOs, variance components were estimated for the variables that were significant in the univariate analysis. The VARCOMP procedure of SAS® (2015) was used with two alternative models: the first one considered the four putative OTU, while the second one grouped the samples of oaxacana and apulcensis in the same OTU.

The discriminatory ability of needle and cone traits was determined using a canonical discriminant analysis (CDA) with each group of variables separately and together. The 'candisc' function (^{Friendly, Fox, & Friendly, 2021}) of the R software v. 3.6.1 (^{R Core Team, 2020}) was used to estimate the canonical variables (CV) and the graphical representation of the canonical structure. Linear functions (discrimination criterion) and their cross-validation were calculated with the 'lda' and 'predict' functions of the MASS package with the same “a priori” probability (25 %) of membership (^{R Core Team, 2020}; ^{Ripley et al., 2021}). Prior to CDA, the STEPDISC procedure (^{SAS, 2015}) was used for the selection of variables from each group with significant contribution (P ≤ 0.05) to the discrimination model.

Morphological differentiation for needles and cones

Table 1 shows that the four P. pseudostrobus OTUs under study were differentiated by needle morphology, as differences were found in 13 of the 20 traits evaluated (P ≤ 0.05); variables related to resin canals and hypodermal cells were not significant (P > 0.1617). Although the coefficient of needle shape, sheath length, number of needles, number of intrusions of the hypodermis into the endoderm and the proportion of resin canals with middle position were significant (P ≤ 0.05), the mean comparison test showed no differences in the average values of varieties, probably due to the imbalance in the sample size of the varieties compared. In the case of stomata per line and total stomata, the values were higher for apulcensis variety compared to oaxacana variety.

Table 1 indicates that cone morphology also showed differences between the variants of the species. Variety had a significant effect on all variables evaluated for cones (P ( 0.0110). Except for closed cone length, in all cases significant differences were found between oaxacana and apulcensis with contrasting mean values for peduncle length, open cone diameter, cone shape coefficient, depth and height of the apophysis. Multivariate analysis on needles and cones showed significant differences between varieties according to Wilks' Lambda, Pillai's trace and Hotelling's trace statistics (P < 0.0001), as well as the contrast between oaxacana and apulcensis, with a Wilks' Lambda of 0.0411 for needle traits and <0.0001 for cone traits.

Variance structure between and within varieties

The variance structure changed when differentiating var. oaxacana from apulcensis. According to Table 2, when considering the four varieties, the variance between OTU (σ² _var) increased between 0.1 and 14.7 % (average 6.8 %) with higher increases in the stomata per line and total needles variables, and in apophysis height, peduncle length and open cone diameter. The increase was most evident in traits with significant variation between varieties (Table 1). The increase in between-group variance and reduction in within-group variance (σ² _arb(var)) (i. e., greater within-variety homogeneity) when subdividing one or more groups is interpreted as evidence of needle and cone morphological differentiation between the two OTUs that were separated in the model. ^{Perry (1991}) discriminated the varieties because oaxacana has a longer, more prominent, and erect cone umbo than var. apulcensis; he also mentions that the needles are shorter and the cones longer in apulcensis, although no differences were found in this study (Table 1).

Important traits for variety discrimination

According to Table 3, the STEPDISC analysis considered seven needle traits and eight cone traits as the most important for discriminating putative OTU. Among the selected cone traits, apophysis height stands out, which is important in the taxonomic identification of the OTU, especially for var. oaxacana (^{Farjon, Pérez de la Rosa, & Styles, 1997}; ^{Perry, 1991}); however, peduncle length, cone shape coefficient and degree of opening, apophysis depth and depth/height ratio, length of closed cone and diameter of open cone were also important in differentiating varieties, which had not been reported in previous studies.

In the joint analysis, 11 variables were selected, all eight from cones and three from needles: total stomata, middle resin canals and endodermal cells, the latter was not selected in the separate analysis (Table 3). Although NL was not important in the joint analysis, this trait and those related to resin canals were useful in the study by ^{Delgado et al. (2007}) to differentiate between putative hybrids of P. pseudostrobus; also, needle number and length, and middle resin canals showed wide variation with a clinal pattern in P. patula (^{Dvorak, Jordan, Romero, Hodge, & Furman, 2001}). The variance explained by unselected needle traits in the joint analysis, especially needle length, may be contained in the cone variables by the correlations observed with the apophysis depth/height ratio (r = -0.42), apophysis height (r = 0.41) and open cone diameter (r = 0.37).

Canonical discrimination analysis

According to Table 4, the canonical correlation (C_r) was significant (P ≤ 0.0271 for all canonical variables (CV). In needle CV₁ explained 58 % of the variance between varieties, while in cone CV₁ explained 72 %, and 68 % in the joint analysis (Table 4; Figure 2). The c_r value between the canonical values of individuals and the first two canonical variables were higher (0.64 ≤ c_r ≤ 0.87) in the cone and joint analysis compared to the needle analysis (0.48 ≤ c_r ≤ 0.59). That is, needle traits showed lower ability to separate varieties (Figure 2a). In contrast, the first two CV of cone traits separate three of the varieties, but there is overlap between var. pseudostrobus and apulcensis (Figure 2b). Although the inclusion of needle traits does not add additional information to that of cones in discriminating the first three varieties, it does reduce the overlap between var. pseudostrobus and apulcensis (Figure 2c).

The canonical structure matrix shows the raw relationship between needle and cone traits with CV (Table 5). The coefficients represent the relative contribution of each trait to CV. In needles, the traits total stomata (TS) and lines with stomata (LS) are the most important associated with CV₁, and needle length (NL) and number of stomata (NS) with CV₂; in cones, the traits apophysis height (H_A) and depth (D_A), as well as their ratio (R_DH) are the most important associated with CV₁, while open cone diameter (CD_O) and cone shape coefficient (CS_C) are important in CV₁ and CV₂. In the joint analysis, no needle trait is relevant for CV₁, but TS and endodermal cells (EC) are important in CV₂ and contribute to the separation of var. apulcensis from the others (Figure 2c).

Figure 2 Structure of canonical variables (CV₁ and CV₂) of morphoanatomical traits of needles (a), cones (b) and as a whole (c) differentiating four putative varieties of Pinus pseudostrobus. Needles: LS = lines with stomata, TS = total stomata, SR = sheath /needle ratio, NS = number of stomata, NN = needles per fascicle, NL = needle length, EC = endodermal cells. Cones: CD_O = open cone diameter, H_A = apophysis height, D_A = apophysis depth, CL_C = closed cone length, CS_C = cone shape coefficient, R_DH = apophasis depth/height ratio, PL = peduncle length, OD_C = openness degree. 

The morphological similarity of cones between var. pseudostrobus and apulcensis may be due to the semi-sympatric relationship they share (^{Carvajal & McVaugh, 1992}) and the wide variability of the former, especially the apophysis height, which sometimes hinders its distinction (^{Farjon et al., 1997}); however, when foliage traits are added, differentiation capacity increases (Figure 1c). It appears that the number of stomata and endodermal cells, associated with the size of the vascular bundle, are useful for distinguishing varieties. On the one hand, stomata density is an important functional trait in the control of water balance and CO₂ capture per unit leaf area (^{Hetherington & Woodward, 2003}); moreover, the size of the vascular bundle (xylem and phloem) affects water efficiency and safety, which is a fundamental basis for the regulation of gas exchange (^{Domec, Palmroth, & Oren, 2016}). In this regard, the results support findings that var. apulcensis occupies part of the natural habitat of the typical variety but is absent in more xeric environments (^{Farjon & Styles, 1997}). Cone traits clearly separate var. oaxacana from apulcensis, so these two OTUs should not be considered synonyms.

Cross-validation results reflected the reduced contribution of needle traits to the discrimination between varieties. Table 6 shows that, when using both groups of traits, similar hit percentages were achieved to those obtained with cone traits alone, with average errors of 7 % and 9.2 %. Cone traits were more efficient in the identification and correct discrimination of the UTOs. The greater simplicity and lower cost of measuring cone traits, compared to foliage traits, compensates for the limitations and costs of cone collection, which is only possible at certain time of the year. Nevertheless, the usefulness of some functional traits of adaptive value in needles, such as total stomata and endodermal cells should not be excluded. This is the case of distinguishing P. patula var. longipedunculata in the Sierra Madre del Sur with morphological markers of adaptive value (^{Dvorak et al., 2001}).

In the cross-validation analysis, var. pseudostrobus presented the lowest percentage of hits with the two groups of traits (Table 6). The difficulty in its differentiation is due to the high morphological variability because of the wide geographic distribution and interaction with other processes such as diversification, followed by gene flow and low rates of evolution (^{Gernandt & Pérez-de la Rosa, 2014}), which cause greater phenotypic variation compared to the other varieties and complicate its taxonomic delimitation (Figure 1). The wide overlap observed in the morphological variation of this variety compared to the others may be due to individuals in intertaxa transition zones or in the process of speciation, because of introgressive hybridization with other species of the same phylogenetic group (^{Gernandt et al., 2009}). ^{Delgado et al. (2007}) were only able to classify 40 % of the putative hybrid individuals of P. pseudostrobus, due to the great variability of morphological traits caused by hybridization with other species

The results of the study confirm the wide intraspecific variation of P. pseudostrobus, coinciding with findings in other phenotypic traits such as seedling performance and survival in different growth environments (^{Villegas-Jiménez, Rodríguez-Ortíz, Chávez-Servia, Enríquez-del Valle, & Carrillo-Rodríguez, 2016}; ^{Viveros-Viveros et al., 2005}; ^{Viveros-Viveros et al., 2006}), frost tolerance (^{Viveros-Viveros, Sáenz-Romero, López-Upton, & Vargas-Hernández, 2007}) and isoenzyme variation (^{Viveros-Viveros, Tapia-Olivares, & Sáenz-Romero, 2014}). The identification of P. pseudostrobus OTUs and the recognition of the extent of intraspecific variation in traits of adaptive value are indispensable requirements for the proper management and sustainable use of the genetic resources of the species.