Pinus L. (Pinaceae) is a monophyletic and conspicuous genus comprising ca. 119 species of large (up to 80 m in height), long-lived, monoecious, and perennial trees distributed almost exclusively in the Northern Hemisphere (Mirov 1967; Richardson and Rundel 1998; Farjon 2010). Pines possess distinctive secondary needle-like leaves arranged singly or in fascicles (Farjon 2005, 2010). Leaf and seed cone characters are usually used to recognize pine species because they are more variable than other structures (Farjon and Styles 1997). Despite often being of great economic and ecological importance, uncertainty still exists on the geographic limits of natural ranges for some pine species, especially in Mexico (Shaw 1909; Critchfield and Little 1966; Perry 1991; Eckenwalder 2009; Farjon 2010; Debreczy and Rácz 2011). Clearer interspecific limits would substantially improve the use, management, and conservation of pines (Mirov 1967; Challenger and Soberón 2008).

Pinus subsection Ponderosae Loudon includes between approximately 14 and 16 species distributed from western Canada and the U. S. A. to Mexico and Central America (Little and Critchfield 1969; Gernandt et al. 2005; Eckenwalder 2009; Farjon 2010; Hernández-León et al. 2013). Divergent opinions surround the circumscription of species in this group (e.g. Price et al. 1998; Eckenwalder 2009; Farjon 2010; Debreczy and Rácz 2011). Morphometric studies have been carried out to better define the specific limits of variable taxa such as P. ponderosa P. Lawson & C. Lawson (Callaham 2013), P. hartwegii Lindl. (Matos 1995), and P. douglasiana Martínez and P. maximinoi H. E. Moore (Stead 1983a). Studies that integrate more diverse sources of data are needed for these and other pines, as are objective criteria for deciding whether or not to recognize taxa.

Pinus tenuifolia Benth. was described from material collected near Guatemala City, Guatemala (Bentham 1842); the name is a later homonym of P. tenuifolia Salisb. (1796), which in turn is considered a synonym of P. strobus L. (Farjon and Styles 1997). Shaw (1909) treated this taxon as a variety of P. pseudostrobus Lindl. (P. pseudostrobus var. tenuifolia Shaw), and described its distribution as from northwestern and central Mexico to Nicaragua. Martínez (1943) segregated P. douglasiana from P. maximinoi based on populations occurring in central and western Mexico with thicker, stiffer needles and thicker seed cone scale apophyses. Leaf and cone characters of P. maximinoi and P. douglasiana were included by Stead (1983a) in subsequent principal components and canonical discriminant analyses that illustrated the differences in leaf and cone size between these taxa. He concluded that they were legitimate species (Stead 1983a, 1983b). According to Farjon and Styles (1997), no other consistent differences have been identified between P. douglasiana and P. maximinoi. Although Martínez (1948) and Mittak and Perry (1979) stated that P. douglasiana and P. maximinoi can be distinguished by rough versus smooth branchlet surfaces, others such as Perry (1991) and Farjon and Styles (1997) have concluded that this character is inconsistent. Notwithstanding the morphological differences in size, Silba (1990) demoted P. maximinoi to a variety of P. douglasiana (P. douglasiana var. maximinoi [H. E. Moore] Silba).

Some pine species may be in a formative stage (Farjon and Styles 1997; Perry et al. 1998). Molecular studies suggest that Mexican and Central American taxa of Pinus subsection Ponderosae diversified relatively recently. A plastid DNA study of the group found only minor differences between P. douglasiana and P. maximinoi, with some P. maximinoi individuals also sharing haplotypes with P. pseudostrobus Lindl., demonstrating the lack of genealogical monophyly in P. maximinoi (Gernandt et al. 2009). Plastid sequences typical of P. pseudostrobus are more closely related to two other Mexican and Central American species of the ‘Montezumae Group’, P. montezumae Lamb. and P. hartwegii (Gernandt et al. 2009). Hybridization, introgression, or incomplete lineage sorting could explain the sharing of plastid haplotypes among pine species (Delgado et al. 2007; Syring et al. 2007; Willyard et al. 2009). A subsequent molecular clock-based estimate based on plastid DNA indicated that P. douglasiana and P. maximinoi shared a common ancestor in the Pleistocene (Hernández-León et al. 2013).

Pinus douglasiana and P. maximinoi are also difficult to distinguish in the field and herbarium from two closely related species in subsection Ponderosae, P. pseudostrobus and P. yecorensis Debreczy and Rácz. These four taxa differ in leaf and cone characters (Table 1). Farjon and Styles (1997) recognized two varieties of P. pseudostrobus (var. pseudostrobus and var. apulcensis). Pinus pseudostrobus var. pseudostrobus is the most similar in appearance to P. douglasiana and P. maximinoi (Table 1). Pinus yecorensis was described more recently (Debreczy and Rácz 1995). It was treated as a doubtful name, possibly a variety of P. pseudostrobus by Farjon and Styles (1997), but recognized as a legitimate species in a recent treatment of the trees of Sonora (Felger et al. 2001). Individuals from the type locality have sequences that belong to the same haplogroup as P. douglasiana and P. maximinoi, but with unique substitutions, which may support the species status of P. yecorensis (Gernandt et al. 2009).

Pinus douglasiana is confined to Mexico while P. maximinoi has a natural distribution that extends from northwest Mexico to Guatemala, Honduras, El Salvador, and Nicaragua (Martínez 1948; Mittak and Perry 1979; Stead and Styles 1984; Perry 1991; Farjon and Styles 1997; Gapare et al. 2001). Pinus douglasiana, P. maximinoi, P. pseudostrobus, and P. yecorensis sometimes occur in sympatry (Stead and Styles 1984; Perry 1991, Farjon and Styles 1997; Felger et al. 2001). However, disagreements exist whether certain populations in northwestern and central Mexico correspond to P. douglasiana, P. maximinoi, or P. pseudostrobus (Stead and Styles 1984; Perry 1991; Farjon and Styles 1997).

Studies of the climatic factors associated with the distribution of pines could help to understand their taxonomy and variation. The species considered here occur across the complex topography and diverse vegetation types of Mexico, Guatemala, Honduras, El Salvador, and Nicaragua. Pinus douglasiana grows in warm to temperate climatic zones at elevations between 1,400 and 2,500 m (Stead and Styles 1984; Perry 1991, Farjon and Styles 1997). Pinus maximinoi occurs principally in climates that favor cloud forests and wet subtropical forest at lower elevations, primarily between 900 and 1,800 m (Mittak and Perry 1979; Farjon and Styles 1997; Dvorak et al. 2000). Pinus pseudostrobus thrives in warm-temperate to cold-temperate climates, commonly at elevations between 1,900 and 3,000 m (Farjon and Styles 1997). Pinus yecorensis occurs on dry, hot, and rocky slopes between approximately 1,300 and 1,600 m (Debreczy and Rácz 1995, 2011). Modeling potential distributions can help determine whether differences in climatic requirements exist between taxa. Depicted potential distribution areas show zones where suitable climatic conditions prevail in order to determine where species occur (Raxworthy et al. 2007; Soberón and Nakamura, 2009).

Pines, like many trees, have high intraspecific variation and low mutation rates per unit time (Petit and Hampe 2006). Pines are further characterized by high outcrossing rates, huge effective population sizes, and weak barriers to gene flow, all of which obscure their interspecific limits (Ledig 1998; Delgado et al. 2007; Willyard et al. 2007). Reproductive isolation, monophyly, diagnosability, and ecological differences are contingent attributes, not an unavoidable result of speciation processes; however, the presence of these attributes often helps differentiate species (de Queiroz 1998; Sites and Marshall 2004).

Our main objective is to integrate three different lines of evidence that could help understand the process of speciation and divergence in recently formed sister taxa; morphology, molecules, and climate. Synthesis of this information is needed to understand complex morphological variation in Pinus. Given the deficiency of fixed qualitative differences for this group, we quantify morphological variation and evaluate it using statistical techniques.

Table 1.

Main diagnosable characters for Pinus douglasiana, P. maximinoi, and closely related taxa. Information based on Stead and Styles (1984), Farjon and Styles (1997), Debreczy and Rácz (1995), and Gernandt et al. (2009). ^a In this work only hypodermal intrusions that come in contact with the endodermis are referred to as such. Pinus maximinoi can lack hypodermal intrusions at its southernmost distribution. ^b Persistent peduncles are those that remain on the branch after the mature cone falls from the tree.

Materials and Methods

Plant Material—Field work was conducted from 2010–2013. To capture the greatest possible morphological variation, field collections were made from representative sites throughout the natural range of P. douglasiana and P. maximinoi in Mexico. Branches and seed cones of the study species were collected from the bottom of the crown. Material was pressed, dried, and deposited in the National Herbarium of Mexico (MEXU). Herbarium specimens from Guatemala, Honduras, and Nicaragua was examined, but no fieldwork was conducted in these countries. A total of 287 individuals were included (Appendix 1). Our sampling for P. yecorensis was limited, but we included it in comparisons whenever possible given its presumed close relationship with P. douglasiana and P. maximinoi. Locality data for 16 collections of P. yecorensis and three collections of P. pseudostrobus were taken from online databases (Appendix 1).

For the molecular study two or three fresh leaves were conserved in a freezer at -20°C. To obtain anatomical data, ∼1.0 cm segments from the medial part of fresh leaf samples were preserved in a FAA solution. For herbarium material, ∼1.0 cm segments from the medial part of a leaf were rehydrated in boiling water and then emersed ∼24 hrs in FAA. Transverse sections were made manually. The sections were cleared in 50% sodium hypochlorite for 2 min, rinsed several times with distilled water, dehydrated in a series of alcohol rinses (50, 75, 96, and 100%), stained with fast green in 100% alcohol, and mounted on a slide with resin. Because the data set was incomplete, not all collections were included in all analyses. The number and details of the individuals included in each analysis are specified below.

Preliminary Taxonomic Identification—A preliminary identification was made of each individual based upon field observations and microscopic examination of the leaves (Table 1). The persistence of the cone peduncle and presence of hypodermal intrusions in the leaf were the principal characters used to distinguish P. pseudostrobus from P. douglasiana or P. maximinoi (Martínez 1948; Mittak and Perry 1979; Stead and Styles 1984; Perry 1991; Farjon and Styles 1997). Length and thickness of the cone peduncle in P. pseudostrobus are variable. Forty-five new collections from northwestern Mexico, principally the western extreme of the Trans-Mexican Volcanic Belt (TMVB; Jalisco and Michoacán) and the Sierra Madre Occidental (SMOC; Sinaloa, Sonora, and Chihuahua) have heterogeneous seed cone morphology, deciduous cone peduncles typical of P. douglasiana and P. maximinoi, and predominantly lack leaf hypodermal intrusions. We refer to them here as Pinus aff. douglasiana. Some individuals identified as P. pseudostrobus lack leaf hypodermal intrusions, but have both persistent and deciduous peduncles and have flat cone apophyses that are similar to P. maximinoi. Collections of P. yecorensis lack hypodermal intrusions; this taxon has subtle differences from P. douglasiana such as stout branches and short and more spherical cones with stout peduncles (Gernandt et al. 2009).

Analysis of Quantitative Variation—A total of 219 individual collections were included in the morphological analyses. Only collections with leaves that were mature in size and appearance were measured; similarly, only collections with open seed cones were measured. Data were obtained for herbarium specimens (mainly from Guatemala, Honduras, and Nicaragua) although it was not always possible to obtain all measurements. The herbarium specimens were included in the analyses given the desirability of evaluating their variation.

Morphological Data Set—Fourteen continuous or meristic characters of the leaf and cone identified as informative by Stead (1983a) and an additional three characters proposed here (characters three, five, and sixteen) were evaluated (Fig. 1): 1) LON: length of needle (cm), based on the mean of five leaves, except for herbarium material which was based on two or three; 2) LOS: Length of sheath (cm), based on the mean of five fascicles, except for herbarium material which was based on two or three; 3) TOS: Thickness of sheath (mm), based on the mean of five sheaths, except for herbarium material which was based on two or three; 4) WON: Width of needle (μm), based on one observation per individual; 5) TON: Thickness of needle (μm), based on one observation per individual; 6) Ncan: Number of resin canals, based on one observation per individual; 7) Int: Number of hypodermal intrusions, based on one observation per individual; 8) NEC: Number of endodermal cells in the leaf, based on one observation per individual; 9) NSLD: Number of stomatal lines on the leaf dorsal surface, based on one observation per individual; 10) NSLV: Number of stomatal lines on the leaf ventral surface, based on one observation per individual; 11) LOC: Length of the (seed) cone (mm), based on one observation per individual; 12) WOC: Width of (seed) cone at the widest point (mm), based on one observation per individual from the largest cone available; 13) WOA: Width of (seed scale) apophysis (mm), based on the mean of five cone scales evenly spaced around the widest point of a single cone, except for herbarium material which was based on two to five scales; 14) HOA: Height of apophysis (mm), based on the mean of five cone scales evenly spaced around the widest point of a single cone, except for herbarium material, which was based on two to five scales; 15) TOAU: Thickness of (seed cone scale) apophysis and umbo together (mm), based on the mean of five scales evenly spaced around the widest point of a single cone, except for herbarium material which was based on two to five scales; 16) LOCS: Length of cone scale (mm), based on the mean of five scales evenly spaced around the widest point of a single cone; 17) PPED: Persistence of the peduncle. This character was coded as absence (0) or presence (1). The morphological data are available from the Dryad Digital Repository:  http://dx.doi.org/10.5061/dryad.931h7.

After preliminary examinations, we concluded that the smoothness of the branchlet surface was not consistent and decided not to consider the character. Measurements of leaves and cones were obtained from dried and pressed material using a ruler calibrated to 0.1 cm and a Vernier caliper (Mitutoyo digimatic plastic caliper) with a resolution of 0.1 mm and an instrumental error of ± 0.2 mm. For the WON and TON variables, measurements were made with a light microscope using a 10× (occasionally 5×) objective lens and an ocular with a micrometric ruler.

One-way ANOVA—All statistical analyses were carried out with R 2.14.2 (R Development Core Team 2011). Individuals with leaf hypodermal intrusions and deciduous cone peduncles from Michoacán and the State of Mexico in the central TMVB were intermediate between P. douglasiana and P. maximinoi with respect to leaf width and cone scale apophysis thickness. To clarify their identity three groups were compared using the WON and TOAU variables: individuals clearly identified as P. douglasiana (n = 55) from Jalisco and Sinaloa, individuals identified as P. maximinoi (n = 69) from Chiapas, Guerrero, and Oaxaca, and the undetermined individuals from the TMVB (n = 21). An analysis of variance (ANOVA) was performed for this comparison. Levene's test (Fox and Weisberg 2011) and Bartlett's test were used to corroborate the supposition of homogeneity of variance (homoscedasticity). A Tukey test was performed when significant differences were found. Non-parametric tests were used when the data did not fulfill the assumptions of the model. In such cases the Kruskal-Wallis test was used to test for differences and the Wilcoxon rank sum test was used for multiple comparisons. Based on this comparison (see results) these individuals were tentatively classified as P. douglasiana.

Multivariate Analyses—We used principal components analysis (PCA) to explore the formation of groups that indicates the existence of taxa (Stead 1983a; Manly 1994). Only individuals with complete data for leaves and mature cones were included. Sixteen quantitative variables were scaled to standardize them. The PPED was excluded because it is a categorical variable. The results from two analyses are reported. The first (PCA-1; n = 145) included P. douglasiana (n = 76) and P. maximinoi (n = 69), including seven specimens of P. maximinoi from Central America that lacked hypodermal intrusions. The second (PCA-2; n = 214) included P. douglasiana (n = 76), P. maximinoi (n = 69), P. aff. douglasiana (individuals from Michoacán and Jalisco that were morphologically intermediate between P. douglasiana and P. maximinoi, but more similar to P. douglasiana; n = 41), P. pseudostrobus (n = 24), and P. yecorensis (n = 4). A correlation matrix was used to perform the analyses. For P. douglasiana and P. maximinoi, we statistically evaluated the existence of gaps in multivariate space to determine whether discrete groups could be identified (Zapata and Jiménez 2012). The analysis was performed using a γ of 0.95 and β of 0.90 for the overlap in the ellipsoidal region of tolerance. The same number of collections was used as in PCA-1.

Clustering methods were also used to explore the formation of groups using the data from the 16 continuous variables and the categorical variable (PPED). The k-means method (MacQueen 1967) was used to form groups. For all these analysis the variable data were standardized and the algorithm of Hartigan and Wong (1979) was used. Two analysis that included all the taxa (P. aff. douglasiana, P. douglasiana, P. maximinoi, and P. pseudostrobus) were conducted. The same collections were used as in the PCA-2 analysis. These analyses were performed including and excluding the variable PPED. For these two analyses a table was constructed to compare the prior classification based on morphology, with the groups obtained with the k-means method.

Linear Model Analysis—Linear regression was used to evaluate the relationship between elevation and morphological variation as represented by the PC1 data from the PCA-1 analysis (P. douglasiana and P. maximinoi). To test the validity of the regression, the normality of the residuals was evaluated with the Shapiro-Wilk normality test; alternatively, homogeneity of variance was inspected visually using diagnostic graphics.

Molecular Study—Because the ycf1 gene (∼5,100 bp) represents one of the most variable regions of plastid DNA in pines (Gernandt et al. 2009; Parks et al. 2009), an ∼800 bp fragment of ycf1, corresponding to amplicon B of Gernandt et al. (2009), was amplified and sequenced. Sequences of this ycf1 fragment were identical for five of eight P. douglasiana and P. maximinoi individuals in a previous study (Gernandt et al. 2009). However, the fragment is variable enough to distinguish a P. douglasiana-P. maximinoi-P. yecorensis haplogroup from another corresponding to P. pseudostrobus. We also chose it for further characterization with the expectation that greater sampling would reveal additional variation in the study taxa.

DNA Isolation and Amplification—DNA was isolated from40–60 mg of ground, frozen needle tissue using the CTAB method (Doyle and Doyle 1987). The ycf1 fragment was amplified by the polymerase chain reaction (PCR) technique using primers reported by Gernandt et al. (2009), Pt98232F and Pt99121R. The following concentrations were used for PCR: 1 × buffer, 1.25mMMgCl2, 0.4 μMprimer, 0.2 μMof each dNTP, and 0.6–1 U Platinum Taq DNA polymerase (Invitrogen, Carlsbad, California). Approximately 5–50 ng genomic DNA was added to each reaction. The amplification conditions were as follows: initial denaturalization at 94°C (3 min); 30 cycles of denaturation (94°C for 1 min), annealing (50°C for 50 sec), and extension (72°C for 80 sec); and a final extension at 72°C (5 min).

DNA Sequencing, Edition, and Analysis—The PCR products were sent to the University of Washington High Throughput Genomics Center (Seattle, Washington) for purification and sequencing. Forward and reverse sequence reads were assembled and edited with Sequencher 4.8. (Gene Codes, Inc. Ann Arbor, Michigan). Assembled sequences were aligned with Clustal W (Thompson et al. 1994), as implemented in BioEdit 7.0.5 (Hall 1999) using default parameters and adjusted by eye. Ninety-five new sequences from P. maximinoi, P. douglasiana, P. pseudostrobus, and P. yecorensis were obtained and analyzed. These were aligned with 12 sequences downloaded from GenBank (Appendix 1). The alignments were trimmed at the extremes adjacent to the primers to avoid sites with missing or low quality base calls. Sequences were used to assemble a matrix 803 bp in length, without indels and without missing data. The haplotype of each individual is available from the Dryad Digital Repository:  http://dx.doi.org/10.5061/dryad.931h7. A haplotype network was built to estimate gene genealogies with the statistical parsimony method implemented in TCS 1.21 (Templeton et al. 1992; Clement et al. 2000).

Fig. 1.

Cone and leaf variables. A. Fascicle. B. Fascicle sheath. C. Cross sectional anatomy of a secondary leaf. D. Stomatal lines on the dorsal face of a leaf. E. Cone. F. Cone scale. G. Lateral view of a cone scale (See methods for an explanation of abbreviations).

Species Distribution Modeling—The potential distributions of P. douglasiana, P. maximinoi, P. pseudostrobus var. pseudostrobus, and P. yecorensis were modeled using MaxEnt 3.1 (Phillips et al. 2006) with the default settings. Climatic and elevational variables, with a spatial resolution of about ∼1 km², were used for the respective distribution modeling. Bioclimatic data for Mexico, Belize, Guatemala, Honduras, El Salvador, Nicaragua, Costa Rica, and Panama were obtained from WorldClim 1.4 (Hijmans et al. 2005). A predictive model based on presence-only data was generated using geographic coordinates of collections from the natural distribution of the study taxa. Only one occurence record per species per grid cell (∼1 km²) was included.

Records with locality descriptions but lacking coordinates were georeferenced using Google Earth 5.1.3509.4636 (beta) and the Mexican website of the Instituto Nacional de Geografía e Informática (INEGI 2014). For P. douglasiana and P. maximinoi records, 17 and 39 points, respectively (Appendix 1) correspond to individuals that we identified with multivariate analyses or herbarium collections. The type locality of P. douglasiana in Cerro Juanacata near the town of Jala in Nayarit, Mexico (Martínez 1943, 1948) was not included in the modeling of P. douglasiana because we failed to locate the species when we visited this area. In agreement with the conclusions of Stead and Styles (1984), we did not find the typical form of P. maximinoi in Sonora and Sinaloa. Points for P. maximinoi from cold and high altitude localities in the eastern TMVB (Vigas-Xalapa, Veracruz and Nevado de Toluca, State of Mexico) were not included because the collections that we examined from these sites lack leaf hypodermal intrusions and were thus redetermined as P. pseudostrobus. Three additional georeferenced points for P. pseudostrobus were obtained from the online botanical database, Tropicos (Missouri Botanical Garden 2015). Locality data were obtained for P. yecorensis from six MEXU specimens, 16 records from the website of the Consortium of Intermountain Herbaria (2014), and one locality from the original species description (Debreczy and Rácz 1995). According to our statistical analyses, the distribution points of P. maximinoi and P. douglasiana included for modeling do not include sites where the two species are in sympatry. To visualize zones of sympatry, the potential distribution of these species was displayed on the same map. To test each resulting model, a partition of the ocurrence localities was made setting the random test option in MaxEnt at 25%. Likewise, the scores of the area under the receiver operating curve (AUC) for the training data were considered as statistical tests to evaluate the performance of maximum entropy modeling. It has been reported that scores greater than 0.90 are indicative of a good performance (Phillips et al. 2006; Baldwin 2009). The inferred potential distribution was trimmed considering the logistic threshold value associated with maximum training sensitivity plus specificity value and the resulting models were displayed graphically in ArcMap 10.0 (Environmental Systems Research Institute, Inc., Redlands, California).

To evaluate the relationship between climate and morphological variation we used k-means to compared the groups formed with the climatic variables, the morphological variables, and the preliminary clustering criteria. The same 19 climatic variables used for modeling and the 16 quantitative morphological variables were standardized for the analysis. One collection per site was included. Sites lacking morphological information were excluded. The morphological data correspond to individual collections identified as P. aff. douglasiana (n = 9), P. douglasiana (n = 12), P. maximinoi (n = 20), P. pseudostrobus (n = 7), and P. yecorensis (n = 2). The corresponding information for each of the climatic variables (layers with a resolution of 0.5 sec) and the sites included in the modeling were extracted using the raster package (Hijmans and van Etten 2012). A table was constructed to compare the resulting classifications.

Results

Comparison of Quantitative Character Means—The comparison of P. douglasiana and P. maximinoi in western and southern Mexico, respectively, with individuals from Michoacán and the State of Mexico in the TMVB indicated that this group is statistically different from typical P. douglasiana or P. maximinoi from other localities. Variances among groups in the TOAU character after transforming the data were statistically different (Levene test F_{2, 142} = 4.75; p = 0.010). A non-parametric test (Kruskal-Wallis test) showed significant differences among groups (χ² = 103.5887; df = 2; p < 0.01; see Fig. 2A). Although statistically different, the mean value of TOAU from the Michoacán and State of Mexico collections was closer to the value of P. douglasiana than to P. maximinoi from Chiapas, Guerrero, and Oaxaca (3.72, 4.45, and 2.29 mm, respectively). For WON, collections from Michoacán and the State of Mexico had needles that were wider on average than P. maximinoi from Chiapas, Guerrero, and Oaxaca but thinner than P. douglasiana from Jalisco and Sinaloa (Fig. 2B). The group variances were not statistically different (Levene test; F_{2, 142} = 2.82 ; p = 0.063) but the residuals were not distributed normally (Shapiro-Wilk normality test, W = 0.97, p = 0.0060). A non-parametric test (Kruskal-Wallis test) showed significant differences among groups (χ² = 83.3541; df = 2; p < 0.01), and a nonparametric pairwise comparisons (Wilcoxon rank sum test) showed that all groups were different (p < 0.01). In this comparison of the WON variable, unidentified individuals from Michoacán and the State of Mexico were more similar to P. maximinoi (Fig. 2B). However, based on the high similarity of TOAU with P. douglasiana we decided to treat tentatively these individuals as this species.

Multivariate Data—The PCA-1 and PCA-2 did not result in discrete groups in multivariate space among species. For the PCA-1 (inclusion of P. douglasiana and P. maximinoi only), the sum of variation explained by the first and second principal components (PC1 and PC2) was less than 66% (Table 2). According to the analysis to distinguish gaps in multivariate space, the overlap in the ellipsoidal region of tolerance at a γ of 0.95 and β of 0.90 is such that no gaps were detected and therefore no discrete groups were found (Fig. 2E). All coefficients of PC1 (PCA-1) had negative values (Table 2), indicating that differences between individuals of P. douglasiana and P. maximinoi explained by PC1 were mostly in size and not form. Most variables, but mainly NEC, TOAU, and TOS, contributed to this size difference (Table 2). In the dispersion graph (Fig. 2C), collections from Michoacán and the State of México grouped with individuals from Jalisco and Sinaloa, which supported their identity as P. douglasiana. For PCA-2 (Fig. 2D), the sum of variation explained by PC1 and PC2 was less than 61% (Table 3). With the exception of number of leaf intrusions, the positive values of the PC1 variables indicated that there were differences in size among taxa (Table 3). The variables HOA, NEC, and TOAU contributed most to this size difference. The negative and positive coefficient values for PC2 indicated that there were differences in form among individuals, principally between P. pseudostrobus and P. douglasiana and P. maximinoi (Table 3). These differences indicated an inverse relationship among the cone variables and leaf variables, mainly the number of leaf hypodermal intrusions.

Table 2.

Coefficients of the first and second component (PC1 and PC2) of the first principal component analysis (PCA-1).

Fig. 2.

Statistical comparisons of morphometric variables. A. Comparison of TOAU for individuals from the TMVB (VB; n = 21; mean = 3.78; sd = 0.45), P. douglasiana (Pdo; n = 55; mean = 4.32; sd = 0.93), and P. maximinoi (Pmi; n = 69; mean = 2.3; sd = 0.55) from other localities; the three groups were significantly different (p₂, 142 < 0.01). B. Comparison of WON among individuals from the TMVB (VB; n = 21; mean = 815.2; sd = 68.6), P. douglasiana (Pdo; n = 55; mean = 959.8; sd = 127.9), and P. maximinoi (Pmi; n = 69; mean = 715.4; sd = 103.4) from other localities; the three groups showed significant differences from each other (p₂ < 0.05). C. Scatter plot of PCA-1, P. douglasiana (black circles; n = 76); P maximinoi (open circles; n = 62); Specimens of P. maximinoi without intrusions (rhombs n = 7). D. Scatter plot of PCA2; P. aff. douglasiana (open triangles; n = 24), P. douglasiana (black circles; n = 76), P. maximinoi (open circles; n = 62), P. pseudostrobus var. pseudostrobus (open squares; n = 24), P. yecorensis (black triangles, n = 4). E. The proportion of individuals (β) of P. douglasiana (dashed line) and P. maximinoi (solid line) at different values of α (which define tolerance regions) and a constant confidence level (γ) of 0.95. The graph shows that at a given (α) there never exists the possibility of finding tolerance regions that include at least 90% of the individuals for each species. The dashed line shows the limit for β = 0.90. F. Linear regression showing the relationship between morphology and elevation P. douglasiana (black circles); P. maximinoi (open circles).

Table 3.

Coefficients of the first and second component (PC1 and PC2) of the second principal component analysis (PCA-2).

The results of the k-means clustering analysis (with k = 3) for the 16 quantitative morphological variables demonstrated a notable similarity between the individuals from the TMVB with those of P. douglasiana in Jalisco and Sinaloa (19 of 21 of these individuals grouped with P. douglasiana). The results for k = 4 (Table 4) considering all the individuals (P. aff. douglasiana, P. douglasiana, P. maximinoi, P. pseudostrobus, and P. yecorensis) and excluding the variable PPED demonstrated that all groups were similar, although there was a notable distinction for P. douglasiana and P. maximinoi. The most distinct taxon was P. maximinoi, with 41% of its individuals separated in a completely exclusive group. Including the variable PPED increased the congruence with the groups formed with the preliminary identifications. In general, there was a high similarity between the P. aff. douglasiana individuals and those of P. douglasiana, P. maximinoi, and P. pseudostrobus. Unfortunately the results from analyses that include PPED must be treated with caution because it was coded as a binary (presence or absence) variable.

Elevation and needle and cone size were positively correlated in the regression analysis (R² = 0.5275; F_{1, 141} = 159.52; p < 0.01). Residuals showed a normal distribution according to the Shapiro-Wilk normality test (W = 0.9942, p = 0.84). Pinus douglasiana, with needles and cone scale apophysis that are thicker than P. maximinoi, is distributed at higher elevations (Fig. 2F).

Table 4.

Comparison between clustering criteria. Clustering with the k-means criterion with k = 4, with and without the variable describing the persistence of the cone peduncle (PPED), and the initial classification. Based on 16 continuous leaf and cone variables and PPED. * Individuals lacking leaf hypodermal intrusions.

Fig. 3.

Haplotype network of ycf1 fragment for 107 individuals. The size of the pie graphs is proportional to the number of individuals.

Molecular Variation—Eight haplotypes were found. In the haplotype network (Fig. 3), 88% of the individuals identified as P. douglasiana had the same haplotype, which was also shared by 72% of the individuals identified as P. maximinoi. All individuals classified as P. aff. douglasiana had the typical haplotype of P. maximinoi and P. douglasiana. No P. pseudostrobus individuals had the haplotype typical of P. douglasiana or P. maximinoi. Exclusive haplotypes of P. yecorensis connected to the typical haplotype of P. douglasiana and P. maximinoi. One individual identified as P. maximinoi from Oaxaca and another from Guatemala (Guatemala), had haplotypes typical of P. pseudostrobus. The latter individual from Guatemala lacked leaf hypodermal intrusions.

Potential Distribution—Scores of the area under the receiver operating curve (AUC) for the training data were greater than 0.90 in all models. The AUC score was 0.990 for P. douglasiana, 0.980 for P. maximinoi, 0.977 for P. pseudostrobus var. pseudostrobus, and 0.989 for P. yecorensis. These high scores indicated that model performance was good (Phillips et al. 2006). The training omission rate and test omission associated with the logistic threshold had a value of 0.0 for all models.

The potential distribution for P. douglasiana was mainly confined to central and western Mexico (Fig. 4A). The results indicated that the major discontinuities in its distribution corresponded to areas of disruption between the SMOC and mountain ranges of the TMVB, but in general, its distribution was fragmented. Its potential distribution also included areas of sympatry in northern Mexico with P. yecorensis and western Mexico with P. maximinoi. The results also indicated some small areas with favorable conditions for P. douglasiana scattered in the mountain ranges of southern Mexico, namely the Sierra Madre del Sur, Sierra Norte de Oaxaca, and Sierra Madre de Chiapas. Throughout these areas it could occur in sympatry with P. maximinoi or P. pseudostrobus (Fig. 4A). Variables that contributed most to distribution modeling of P. douglasiana were elevation (24.6%) and temperature seasonality (23%). Considered alone, the mean temperature of the warmest quarter (19.0°C on average) was the most useful for explaining the distribution.

Fig. 4.

Real and potential distributions for the study group. A. Pinus douglasiana (n = 24), P. maximinoi (n = 41), and P. yecorensis (n = 20). B. Pinus pseudostrobus var. pseudostrobus (n = 16).

The potential distribution of P. maximinoi was located principally along the coastal slopes that extend to the Pacific Ocean in Mexico (Fig. 4A). In Central America, areas of distribution of P. maximinoi were predicted in small mountain ranges. Areas of potential distribution of P. maximinoi were fragmented and disrupted mainly among the mountain ranges of western and southern Mexico (Fig. 4A). Elevation contributed most to the distribution model (31.5%), followed by annual precipitation (22.9%). Considered alone, temperature seasonality was the most useful for explaining the distribution, contributing 11.9% to the model.

The potential distribution of P. pseudostrobus var. pseudostrobus was wide and occurred mainly in cold-temperate to warmtemperate zones (Fig. 4B). Its potential distribution occurred throughout the major mountain ranges, from northern Mexico to northern Nicaragua. The variables that most contributed to model performance were temperature seasonality (36%), elevation (27.1%), and mean temperature of the warmest quarter (17.7°C; 16.2%). Considered alone, the maximum temperature of the warmest month (25.6°C on average) was the most useful for explaining the potential distribution of this taxon.

The potential distribution of P. yecorensis was in the northeastern SMOC, where dry climates prevail (Fig. 4A). Precipitation seasonality contributed most to the distribution model for this species (31%), followed by precipitation of the warmest quarter (478.26 mm on average; 23.8%), and elevation (19.6%). Temperature seasonality and mean temperature of the driest quarter, each considered alone, were the most useful for explaining its potential distribution.

For the k-means analysis, taking into account the climatic variables we found a 100% correspondence between the groups formed with those based on the same method and morphology (without including the variable PPED; Table 5). Groups formed with the k-means analysis and morphology alone were not congruent with groups formed based on our a priori determinations.

Table 5.

Comparison of the classification taking into account the climatic variables and the k-means method versus the classification of groups that takes into account morphological variables, k-means, and the initial classification.

Discussion

Taxonomic Limits—Despite their morphological similarities, Pinus douglasiana and P. maximinoi have been considered distinct species in most recent taxonomic treatments (e.g. Price et al. 1998; Eckenwalder 2009; Farjon 2010; Debreczy and Rácz 2011). An exception is the demotion of P. maximinoi to P. douglasiana var. maximinoi by Silba (1990). The quantitative characters helped to differentiate partially P. douglasiana and P. maximinoi, but the two taxa share many characters and we were unable to identify any that were consistently diagnostic. The individuals from the TMVB with intermediate leaf width and cone scale apophysis thickness were more similar to P. douglasiana, but our decision to classify them as this taxon could be considered subjective due to the overlapping range of morphological variation with respect to typical individuals of P. maximinoi.

Overlapping variation in continuous characters has been observed previously in Pinus subsection Ponderosae (Matos 1995). Many species in the subsection may be of recent origin (Hernández-León et al. 2013), and many are partially sympatric (Farjon and Styles 1997), which may have delayed their divergence (Matos and Schaal 2000; Delgado et al. 2007; Willyard et al. 2009). Stead (1983a), sampling fewer populations in western Mexico, found groups that were better defined than the groups found in this study (Fig. 2D). Although our results are partially congruent with his, our PCA results are not completely comparable. In this study we did not sample as intensively within populations; we usually sampled fewer than 25 individuals at each locality, and did not always obtain mature cones. We also only considered continuous variables in our analyses. In contrast, Stead (1983a) sampled as many as 25 individuals per population, and coded the presence of a peduncle and the roughness of the branchlets as qualitative characters. There are other notable aspects that highlight the differences between our results and those of Stead (1983a, 1983b). These are related with the fortunate fact that current techniques allow the consideration of molecular characters. For example, we identified individuals as P. aff. douglasiana that aside from having an external morphology similar to P. douglasiana and occasionally to P. maximinoi, had the typical haplotype of these two species. Individuals of P. aff. douglasiana are similar to P. yecorensis in lacking leaf hypodermal intrusions and a persistent peduncle. Their morphological variation is heterogeneous, but at present we treat them as P. douglasiana. In this respect we differ from Stead (1983a; 1983b) in interpreting the presence of leaf hypodermal intrusions as an inconsistent character for diagnosing P. douglasiana and P. maximinoi. Furthermore, the individuals treated as P. aff. douglasiana may be hybrids formed between P. douglasiana and P. pseudostrobus. These were collected from the TMVB and SMOC in areas between 1,400 and 2,400 m elevation where P. douglasiana and P. pseudostrobus may be in contact. Although P. aff. douglasiana included forms that were morphologically intermediate, hybrids do not always take intermediate forms (Rieseberg 1995). Detailed studies are needed to explore the relationship between P. aff. douglasiana and other closely related species such as P. devoniana or P. pseudostrobus.

Whang and Pak (2001) reported that the stomatal apparatus of leaf cuticles are rectangular in P. maximinoi but elliptic in P. douglasiana, epidermal cells of P. douglasiana have vertical end walls while those of P. maximinoi are vertical and oblique, that the shape of the anticlinal walls in P. douglasiana are straight and sinuous while those of P. maximinoi are mainly straight, and that P. douglasiana has 7–9 epidermal cell rows between stomatal rows while P. maximinoi has 12–14. These characters and others merit further study.

Interspecific Gene Flow—Most individuals of P. douglasiana and P. maximinoi that we included have the same haplotype, and the same occurs with P. pseudostrobus. The typical haplotype of P. pseudostrobus is most closely related to haplotypes of P. montezumae and P. hartwegii, while the typical haplotype of P. douglasiana and P. maximinoi are more closely related to P. devoniana and other species in the Devoniana clade (Gernandt et al. 2009). Hybridization, introgression or incomplete lineage sorting could explain why other individuals of P. maximinoi and P. douglasiana share plastid haplotypes with P. pseudostrobus. Hybridization and introgression are well documented in closely related species of Pinus (reviewed by Ledig 1998; recent examples in North American pines include Matos and Schaal 2000; Delgado et al. 2007; Liston et al. 2007). A molecular study of multiple unlinked loci would be helpful for determining whether these patterns are due to hybridization and introgression or incomplete lineage sorting in these taxa.

Geographic Distribution—The modeling of potential distribution of P. douglasiana, P. maximinoi, P. pseudostrobus, and P. yecorensis was satisfactory, but we should treat these results with caution given that the optimal conditions were not completely fulfilled in the performance of the modeling (Elith et al. 2011). Ecological niche models offer great promise for the development of hypotheses regarding the distribution of species, but they are subject to a series of limitations (Baldwin 2009; Soberón and Nakamura 2009). Models such as those used here only correlate known presence points with a subset of climate variables. Species distributions can also be determined by their biotic interactions with other species and the capacity that they have for dispersing to other regions (Soberón and Peterson 2005). In our case the sampling of collection sites included in the modeling cannot be considered a completely random sample given that, as usually happens, sites with difficult access are poorly represented. We conclude that the species considered here have a fragmented distribution and that this distribution is not well enough known to design a completely random sample. Another disadvantage in the number of sites used for modeling is that the sample size was small for all taxa except P. maximinoi. In this case we can feel somewhat confident because Maxent has been demonstrated to perform satisfactorily with small samples (Phillips et al. 2006; Baldwin 2009). A more accurate potential distribution may be established as we gain a clearer criterion for delimiting P. maximinoi and P. douglasiana. In this sense a completely reliable distribution can only be proposed when consistent diagnostic characters are found. We believe that the morphological differences between P. douglasiana and P. maximinoi, principally in the dimensions of the leaves and cones, are correlated with the environment where they prosper, and as a result we cannot reject the results of the potential distributions. In this case there is a notable similarity in the potential distribution of P. douglasiana and P. maximinoi with the areas of distribution of conifer forests, principally pines, in Mexico (Rzedowski 1990; INEGI 2014).

Disagreements exist regarding the distribution of P. douglasiana and P. maximinoi in central Mexico. Stead and Styles (1984) stated that P. douglasiana and P. maximinoi can occur in sympatry in southern Michoacán, but that only P. douglasiana occurs in central Mexico, which contradicts statements by Martínez (1948), Mittak and Perry (1979), and Perry (1991). This confusion can be due to the difficulty in distinguishing between P. douglasiana and P. maximinoi in Michoacán and the State of Mexico; the individuals from these states tend to be intermediate with respect to the typical individuals of P. douglasiana and P. maximinoi in the width of their needles (WON) and thickness of their cone scales (TOAU). The results of the potential distribution modeling showed that the predominant taxon in central Mexico is P. douglasiana, although P. maximinoi can occur in sympatry with P.douglasiana in southern Michoacán and dispersed throughout the State of Mexico (possibly at low elevations with hot and humid climates). Pinus maximinoi is susceptible to frost (Mitchell et al. 2013), which makes it unlikely that it could be distributed at higher elevations throughout the mountains of Central Mexico. Competition with species better adapted to tropical conditions may limit the southward distribution of P. douglasiana.

At present, the lack of qualitative diagnostic characters and the overlap of morphological quantitative characters between P. douglasiana and P. maximinoi prevent a clear taxonomic delimitation between the two. Differences in leaf and cone size could be the result of environmental influence on the morphology of these species. For example, these differences in size could be correlated with temperature. Temperature seasonality has an important influence on the distribution of the four species. It contributes substantially in the distributional modeling since the AUC values decrease when the model is re-evaluated with the permutated values on training presence and background data for this variable. Furthermore, in Central America P. maximinoi shows a clinal decrease in the number of hypodermal intrusions and a reduction of its genetic diversity and seed size from north to south (Stead 1983a; Dvorak et al. 2002).

Speciation and Selection in a Topographically Complex Landscape—Throughout the ranges of P. douglasiana, P. maximinoi, P. pseudostrobus, and P. yecorensis, geological events have given rise to a complex topography and complex climatic patterns that permit secondary contact, hybridization, and introgression (Ferrusquía-Villafranca 1993; Perry et al. 1998). Furthermore, pines in Mexico follow a distributional pattern determined by elevation and climate (Perry et al. 1998), and phenotypic plasticity in species could be an expression of plant tolerance to climatic changes. Leaf longevity could be strongly related with habitat water- and nutrient relation or stress (Richardson and Rundel 1998). Size differences between P. douglasiana and P. maximinoi could have originated by climatic changes correlated with differences in elevation and latitude. This same phenomenon could explain the emergence of P. yecorensis in northwestern Mexico. Pinus yecorensis has unique plastid haplotypes and although its distribution remains poorly understood, some ecological differences with respect to P. douglasiana and P. maximinoi have been identified, such as inhabiting a drier climate. Only P. yecorensis has both exclusive plastid haplotypes and a distinct combination of leaf and seed cone characters. It might occur in parapatry with P. douglasiana so it is possible to expect limited gene flow between the two. Detailed population studies of morphological and molecular variation in P. yecorensis are needed to determine its taxonomic status more conclusively.

If P. douglasiana, P. maximinoi, and P. pseudostrobus have achieved reproductive isolation, monophyly, or diagnosability, we have been unable to demonstrate it. Most individuals collected across a wide distribution range and identified as P. maximinoi or P. douglasiana share a haplotype (Fig. 3), and studies to date have concluded that variation in peduncle persistence, presence of leaf hypodermal intrusions, and the number and position of resin canals seem to unite more than separate P. douglasiana and P. maximinoi. We advocate treating P. maximinoi and P. douglasiana as separate species, as our data suggest that they are in the early stages of ecological and morphological divergence. Further study of these taxa using a wider diversity of morphological and molecular markers may help us to understand early patterns and processes of speciation in recently diverged conifers.