Study area

The Cupatitzio river basin forms part of the hydrological region of the Balsas River (RH18), within the Tepalcatepec-Infiernillo (I) hydrological sub-region (^{Inegi, 1985}; ^{DOF, 2011}). It is located in the western part of Michoacán State; it is made up of the of Uruapan, Nuevo Parangaricutiro, Gabriel Zamora, Francisco Javier Múgica and Parácuaro municipalities (^{DOF, 2011}). It is located between 18º49'58" and 19º36'11" N and 101º59'30" and 102º13'16" W, with an approximate area of 2 659.03 km² (^{Conagua, 2009}) (Figure 1).

Cuenca del Cupatitzio = Cupatitzio river basin; Elevación = Elevation; Clasificación de arroyos = Stream classification; Orden = Order, Zonas urbanas = Urban sites; Zonas de estudio = Study areas; Alta = High; Media = Medium; Baja = Low.

Figure 1 Geographic location of Pinus populations collected in the Cupatitzio river basin in Michoacán. 

The dominant vegetation is the coniferous forest with three types: pine forest, pine-oak forest and mountain mesophilic forest (^{Bello and Madrigal, 1996}). Soils are of the Andosol type in the upper part of the basin, Luvisol and Acrisol in the middle part, and Regosol and Vertisol in the lower part (^{Inegi, 1985}). The climate corresponds to the temperate climate, humid semi-warm type with abundant rainfall in summer (A)C(m)(w) (^{Inegi, 2010}).

The study area was divided into three geographical zones due to the different pine species: the upper part or recharge area of the basin, with three populations; Cherán (CH1), represented by Pinus montezumae Lamb, Quinceo (Q2) and San Lorenzo (SL3), represented by Pinus pseudostrobus Lindl. The middle part, which includes the Barranca del Cupatitzio National Park Reserve, with three populations: high reserve (RA4), medium reserve (RM5) and low reserve (RB6), represented by Pinus douglasiana Martínez. The third zone, the lower part of the basin, includes La Pinera (LP7) and Tzaráracua (TZ8) sites, represented by Pinus oocarpa Schiede, with eight sampled populations (Table 1 and Figure 1). It is worth mentioning that the town of Cherán is not within the Cupatitzio basin, however, as it is part of the upper recharge zone of aquifers, it was important to include it in the actual study.

Foliar material was collected for each species from five to 19 trees per population, based on the density of the trees, with a total of 102 samples, which were stored in plastic bags and placed in an ultra-freezer at -80 °C, until their use.

Molecular data

The extraction of DNA was done with the CTAB mini-prep method for pines (^{Delgado et al., 2013}). Five core microsatellites (SSRn) isolated from Pinus taeda L. (^{Karhu, 2001}) were tested, of which four were selected that were polymorphic and reproducible (PtTX3013, PtTX3025, PtTX2146, PtTX2123). The amplification of microsatellites was in accordance with ^{Elsik et al. (2000)}, with changes in the concentration of magnesium chloride (4mM) (^{Delgado et al., 2013}). The amplified fragments were separated by electrophoresis in polyacrylamide gels (7 % 8M urea, 0.5 % TBE buffer), run in manual sequencing chambers at 60 Volts, for 2-3 h from the size of the fragments (^{Karhu, 2001}). To observe them, the gels were revealed with silver nitrate (^{Delgado et al., 2013}), and their size was determined by means of a 10 bp marker (Invitrogen^TM), as a reference.

Data analysis

The estimators of the genetic variation that were obtained were the following: percentage of polymorphic loci (% P), total number of alleles (A), average of alleles per locus (N _A ), observed heterozygosis (H _o ) and expected (H _E ) (^{Nei, 1987}) with the formula:

The index of fixation or inbreeding (F_IS) that was obtained as:

The deviation of the Hardy-Weinberg equilibrium was estimated with the nonparametric method of Marco and Monte Carlo chains (MCMC) with 100 000 permutations (^{Guo and Thomson, 1992}). The analyzes were obtained with the program Arlequin v. 3.5.1.2 (^{Excoffier and Lischer, 2010}). Exact linkage disequilibrium tests were performed for each pair of loci for the populations of each species and a global one with Fisher's exact method, obtained with the GENEPOP v. Program. 4.4 (^{Rousset, 2008}). For the four SSRn used, the frequency of null alleles was determined with the FreeNa program (^{Chapuis and Estoup, 2007}).

The effective size (N _e ), that is, the number of individuals that contribute to reproduction (^{Hartl and Clark, 1997}), was estimated indirectly according to the formula:

modified for microsatellites (^{Excoffier and Lischer, 2010}). Since the parameter θ = 4Neμ (^{Kimura, 1968}), and assuming a mutation rate (μ) of 10^-3 per generation, the N _e was estimated as θ / 4X10_-3 (^{Boys et al., 2005}; ^{Delgado et al., 2011}; ²⁰¹³; ^{Rebolledo et al., 2018}).

The genetic structure was determined with an analysis of molecular variance (AMOVA), between and within the populations of each species, with the use of the models of infinite allele mutation (AMI) (^{Weir, 1996}) and step by step (SMM) (^{Slatkin, 1995}). The statistical significance was obtained with 1 000 nonparametric permutations (^{Excoffier and Lischer, 2010}). The parameter M, which is the migration rate (M = Nm) between pairs of populations, was obtained with the formula M = (1 / F _ST -1) / 4 (^{Slatkin, 1993}). All estimates were made with the Harlequin program 3.5.1.2 version (^{Excoffier and Lischer, 2010}).

To analyze the genetic relationships of the populations, a Neighbor-joining dendrogram was constructed, based on standardized genetic distances (Da) (^{Nei et al., 1983}). The robustness of the topology was evaluated with a bootstrap value of 10 000 (^{Takezaki and Nei, 1996}), with the use of the Poptree2 program (^{Takezaki et al., 2010}). The association between the geographic and genetic structure of the populations was obtained through a Bayesian cluster analysis with the BAPS program (^{Corander et al., 2008}). This method uses information from the spatial distribution of populations to define the most probable genetic structure. They were assumed from 1-8 groups (K), with 10 replicates per K with 10 000 iterations, preceded by 10 000 initial iterations.

Finally, an isolation analysis was made by distance between pairs of populations and their geographical distribution, by means of a Mantel test with 1 000 permutations (IBD) (^{Bohonank, 2002}). The paired values of the F_ST statistic (^{Slatkin, 1995}) and the absolute geographical distances in kilometers obtained through the Mercator transformation with the Quantum GIS Geographic Information System (v. 2.16 Nødebo) program (^{QGIS Development Team, 2016}).

Diversity and genetic structure

A total of 37 alleles were recorded in a range of 13 (PtTX2146) to seven (PtTX3013) per locus. P. montezumae and P. pseudostrobus had the highest number of alleles (18 and 17 respectively), followed by P. douglasiana (14) and P. oocarpa (13). The total number of alleles between the populations of each species was different: those of LS3 of P. pseudostrobus with 23 and RB6 of P. douglasiana with 19 had the highest number; Q2 and RA4 of the same species, the smallest number, with 11 and nine alleles respectively (Figure 2, Table 2). They were detected from one to four unique alleles (present in a single population), with frequencies between 0.033-0.150. The RB6 population of P. douglasiana had three unique alleles, CH1 of P. montezumae, Q2 and SL3 of P. pseudostrobus and TZ8 of P. oocarpa had two alleles, while the rest of the populations did not identify any.

Frecuencias alélicas = Allele frequencies; Tamaño de los fragmentos = Size of fragments; ID de las poblaciones = ID of the populations.

Figure 2 Allele frequencies for the four loci analyzed in the Pinus genus populations of the Cupatitzio river basin, Michoacán, Mexico. 

When comparing the results obtained with those of other studies conducted with SSRn (^{Delgado et al., 2011}, 2013; ^{Karhu et al., 2006}), the number of alleles was higher than those of a population of P. pseudostrobus (19) and another of P. montezumae (18) distributed in Michoacán (^{Delgado et al., 2013}), the only two populations of P. caribae var. hondurensis (Sénécl.) W. H. Barret & Golfari for Mexico (17) (^{Delgado et al., 2011}) and five populations of P. radiata D. Don in California, United States of America and Baja California, Mexico (8) (^{Karhu et al., 2006}). P. oocarpa was the only species with an average value lower than those recorded, except for the study of P. radiata (^{Karhu et al., 2006}). These figures indicate that although the populations of the Cupatitzio basin are fragmented, they still preserve allelic variants, so it could be considered that the gene drift has not acted in an important way in the fixation or loss of this form of diversity in the populations.

The H _o on average was lower than the H _E for the four species analyzed. The average of H _o in P. pseudostrobus and P. douglasiana were the highest (H _o = 0.378 and 0.368, respectively), while the lowest values were calculated for P. montezumae (H _o = 0.338) and P. oocarpa (H _o = 0.275). In contrast, H_E values differed between species; P. pseudostrobus recorded the highest value, followed by P. montezumae, P. oocarpa and P. douglasiana (Table 2). The above shows that the values deviate from the Hardy-Weinberg equilibrium, and were significant for the four species (p <0.05) (Table 2).

At the population level, the average of the Ho was also lower than the H _E . The highest values were found in the SL3 populations (H _o = 0.449, H H _E = 0.674) of P. pseudostrobus and RB6 (H _o = to 0.580, H _E = 0.615) of P. douglasiana, while the lowest were obtained for the population of RA4 (H _o = 0.189; H _E = 0.240) of P. douglasiana and P7 of P. oocarpa with a much lower H _o (0.250) than the H _E (0.516) (Table 2). Most values deviate significantly from the Hardy-Weinberg equilibrium (p <0.05), suggesting different intensities of inbreeding.

In general, the results reflect a heterogeneity in the amount of genetic variation within the different populations and zones. However, the values are within those reported for other pine species, such as P. oocarpa (H _E = 0.642) (^{Dvorak et al., 2009}), P. patula Schiede. exSchl. & Cham and P. tecunumanii Eguiluz et Perry (H _E = 0.567 and 0.646 respectively) (^{Dvorak et al., 2009}), P. pseudostrobus and P. montezumae (H _E = 0.277 and 0.274 respectively) (^{Delgado et al., 2013}), P. mugo Turra (H _E = 0.440) (^{Zukowska and Wachowiak, 2017}), P. sylvestris L. (H _E = 0.772) (^{Bilgen and Kaya, 2017}) and P. caribaea var. hondurensis (H _E = 0.575) (^{Rebolledo et al., 2018}), indicating an intermediate and sufficient genetic diversity for the maintenance and possible adaptation of the species to environmental changes.

Endogamy and effective size

The average of the inbreeding index (F _IS ) was high and significant (F _IS = 0.286), so it deviates from the Hardy-Weinberg equilibrium (Table 2). The populations with the highest levels were LP7 (F _IS = 0.393 **, P = 0.000) and TZ8 (F _IS = 0.544 **, P = 0.000), both of P. oocarpa. These results suggest that populations of the Cupatitzio basin are not in genetic equilibrium, probably due to the effect of reproduction among close relatives. The values are relatively higher than those obtained with the use of SSRn in other pine species (P. patula, F _IS = 0.148, P. tecunumanii, F _IS = 0.109) (^{Dvorak et al., 2009}); P. caribaea var. hondurensis, F _IS = 0.163 (^{Delgado et al., 2011}). This shows that although populations contain important levels of genetic variation (such as higher number of alleles or high levels of H _E ), there are processes such as inbreeding that cause the frequencies of genotypes and alleles to deviate from genetic balance (^{Hartl and Clark, 1997}; ^{Holsinger, 2000}).

The high values of inbreeding, and therefore, the deficit of heterozygotes observed in the population, can be explained by the presence of null alleles (alleles not expressed) that give rise to erroneous values that deviate populations from the panmixia (^{Shinde et al., 2003}). In the analysis for the detection of null alleles, only the PtTX3025 locus in the TZ8 population of P. oocarpa showed a frequency higher than 0.200 (0.295). However, for this species, significant inbreeding values were also detected (F _IS = 0.150, P <0.05) with the use of SRNS (including the PtTX3025 locus) for 51 populations distributed in Mexico, Guatemala, Honduras and Nicaragua (^{Dvorak et al., 2009}). The TZ8 population had a significant inbreeding coefficient (F _IS = 0.204, P <0.05), so it is pertinent to consider that the deviation of panmixia from this locus is determined by the effect of inbreeding and not by the frequency of null alleles. Likewise, no evidence of linkage disequilibrium was recognized for any of the pairs of loci analyzed in each species or with the global test (P> 0.05), which indicates that those used are segregated independently, and, therefore, they do not overestimate the F _IS values obtained.

The effective size (N _e ) was again higher in P. pseudostrobus, followed by P. montezumae, P. oocarpa and P. douglasiana. However, the values at the level of the populations of each species were contrasting; the SL3 population of P. pseudostrobus presented the highest N _e (2 104), while the Q2 population had one of the lowest values (N _e = 360). In P. douglasiana the RB6 population had a high value (N_e = 1 439), and in the RA4 population the lowest value was obtained (N _e = 183). In the two populations of P. oocarpa, with the highest inbreeding levels, the estimates of N _e were intermediate (Table 2). Therefore, it could be considered that the resulting N _e estimates are enough to preserve the genetic variability of the species in the Cupatitzio basin.

Currently there are few studies that have estimated the N _e with the use of SSRn in pines, as is the case of Pinus resinosa Ait, distributed in the United States of America, with fragmented populations (N _e between 62-222) (^{Boys et al., 2005} ), of P. caribaea var. hondurensis (N _e between 402-868) (^{Delgado et al., 2011}), with only two populations fragmented and restricted to the southeast of Mexico, as well as some populations of this variety and P. caribaea var. caribaea Morelet, and P. caribaea var. bahamensis (Griseb.) W. H. Barret, distributed in the Caribbean basin (N _e between 161-537) (^{Rebolledo et al., 2018}).

When comparing these results with those obtained in the study, it is observed that most of the values are relatively higher. It is therefore very likely that the fragmentation of some of the populations such as RB6 (N _e = 1 439) is very recent and that the effects of genetic drift or inbreeding are not yet reflected in the reduction of their genetic variability, and therefore in the estimates of their effective sizes.

Genetic structure and genetic flow

The estimation of the genetic differentiation between the populations of P. pseudostrobus and P. douglasiana was higher and significant with the F _ST statistic (0.152 and 0.131, respectively) than with R _ST , although for P. pseudostrobus both indexes were significant (P <0.05) (Table 3). It is likely that for these species migration makes a greater contribution to the genetic differentiation of populations compared to the mutation, which does not follow a strict model (^{Crow and Aoki, 1984}; ^{Balloux and Goudet, 2002}). For P. oocarpa the F _ST statistic was lower (0.094) to R _ST (0.295), although both were significant (P <0.01). In this case, the effect of the mutation is not negligible with respect to migration and could be the most probable cause of the differentiation observed (^{Balloux and Goudet, 2002}) (^{Hardy et al., 2003}).

At present, there is no consensus on the accuracy of the estimates with both mutation models, but it has been experimentally and empirically tested, that it depends, to a large extent, on the size of the populations, the sample size and the number of analyzed loci. (^{Gaggiotti et al., 1999}; ^{Balloux and Goudet, 2002}), with a greater sensitivity of the R _ST statistic to small sample sizes or to the use of few loci (<8), because their estimates have a higher variance than with F _ST (Golstein et al., 1995; ^{Slatkin 1995}; ^{Balloux and Goudet, 2002}; ^{Hardy et al., 2003}). For this last statistic, the number of populations does not have to be infinite, nor necessarily very large for the genetic structure, and, therefore, the estimates of Nm, to be more accurate (^{Niegel, 2002}).

For the estimation of the genetic structure of Mexican pine species through the use of SSRn in particular, the F _ST statistic was used, which resulted in values between 0.059 in P. patula and 0.131 in P. oocarpa (^{Dvorak et al., 2009}; ^{Alfonso et al., 2014}); in the present study higher values were obtained, except for P. oocarpa. Thus, the estimated Nm for the latter species is relatively higher (Nm = 2.41) compared to P. pseudostrobus and P. douglasiana (Nm = 1.40 and Nm = 1.66 respectively).

In works for Mexican pine species, the Nm values calculated with SSRn vary between 1.81 in P. patula and 13.2 in P. caribaea var. hondurensis (^{Dvorak et al., 2009}; ^{Delgado et al., 2011}). However, it is documented that units greater than 1.0 Nm are sufficient to homogenize the allelic frequencies among the populations of some plant species and prevent genetic differentiation (^{Hartl and Clark, 1997}).

Conifers, in general, have dispersion mechanisms by wind, both pollen and seeds, so the estimated genetic flow is generally greater than 1.5 units (^{Hamrick et al., 1995}). Although they also have the capacity to self-pollinate since most of them are mono- and self-compatible species (^{Williams, 2008}). Therefore, it can be considered that, although some of the populations show important levels of inbreeding, they have the possibility of exchanging allelic variants (by the values of Nm obtained), which prevents a marked genetic separation, besides that they contain alleles that are shared among the populations of each zone, as well as differential alleles among the populations of the species of interest.

Genetic relationships and isolation by distance

The Neighbor-Joining analysis formed two groups: the first one represented by the populations distributed in the lower area of the sub-basin (P7 and TZ8), of P. oocarpa; and the second, located in the basal part of the dendrogram, formed by the population of P. montezumae (CH1) and one of the two populations of P. pseudostrobus (SL3) from the upper area of the sub-basin. The most derived or recent populations are composed of those at the middle zone, represented by P. douglasiana and a population of P. pseudostrobus (Q2) (Figure 3).

Zona Baja = Low zone, Zona Alta= Upper zone; Zona Media = Middle zone.

At the base of the branches are indicated the robustness percentages of the topology obtained with 1 000 re-sampling replications.

Figure 3 Neighbor-Joining Dendrogram, based on the genetic distances (Da) of ^{Nei et al. (1983)} among the populations of Pinus genus of the Cupatitzio river basin.  

This result is very similar to that obtained with the Bayesian BAPS cluster analysis, which formed four groups (K4), with the highest marginal probability (LogML = -960.15, P = 1.000) (Figure 4). The first group represented by the CH6 population of P. montezumae; the second by the Q2 population of P. pseudostrobus; a third group formed by the SL3 population of P. pseudostrobus and the three populations of P. douglasiana; and the fourth group formed by P. oocarpa. They only differed in the second group formed by the population Q2, which in the dendrogram corresponds to the most derived population.

Figure 4 Bayesian BAPS grouping, obtained for the Pinus genus populations distributed in the Cupatitzio river basin. In color bars, the 4K groups selected with the most optimal partition of the data are represented (P = 1.000). 

Both types of analysis indicate a tendency to group populations with their geographical distribution, which are represented by P. montezumae and P. pseudostrobus in the upper part of the basin, by P. douglasiana in the middle part and P. oocarpa in the lower part.

This association was partially verified with distance isolation analysis (IBD), with a moderate but significant correlation (r² = 0.443, P <0.026), in which 44.3 % of the genetic differences observed are attributable to the geographical distance between populations.

It has been shown that IBD is more evident in some pine species that have a fragmented distribution, which have been obtained in several species analyzed with chloroplast microsatellites (SSRcp) (^{Delgado and Piñero, 2008}; ^{Moreno and Piñero, 2009}; ^{Jardón et al., 2011}).

Only three papers have been reported with the use of SSRn in pines, showing a low but significant association: Pinus pinaster Ait in the Mediterranean basin in Europe (r = 0.360) (^{Mariette et al., 2001}); P. resinosa that is located in the Northeast of the USA (r = 0.381) (^{Boys et al., 2005}); and in the Caribbean pine complex for the Caribbean basin (r = 0.263) (^{Rebolledo et al., 2018}). These results are similar to that obtained in the present study, so that a large part ((60 %) of the distribution of genetic variation could be linked to other factors, such as spatial fragmentation of populations, demarcated by both natural geomorphology of the basin, as by the change of land use, from forest to avocado or urbanization orchards (^{Bravo et al., 2009}). This relationship is also influenced by the intrinsic genetic differences of the four taxa studied, which is observed in the analysis of Bayesian clustering and in the dendrogram obtained.

Recommendations for conservation

First, it is suggested to carry out in situ conservation activities of the SL3 populations of P. pseudostrobus and RB6 of P. douglasiana, distributed in the upper and middle zones of the basin respectively, since they were the populations that presented the highest levels of H _E and of allelic variants, which could contain the representative gene pool of both species. The RB6 population of P. douglasiana is located within the Reserve area of the Barranca del Cupatitzio National Park, which has helped its conservation, while the SL3 population of P. pseudostrobus, which is distributed in the community of San Lorenzo, requires of mediate actions to delimit it in a formal way, which could be declared by the competent authorities as a Forest Management Unit (UMAFOR; ^{Conafor, 2009}), in such a way that both populations function as reservoirs of biodiversity and connectivity bridges with other populations or fragments (^{Aavick et al., 2014}).

Likewise, the collection of seeds from these two populations is recommended since they contain the greatest genetic variability. This type of seeds should be used for large-scale plant production for restoration or reforestation purposes; particularly for the case of the Purépecha Plateau, where the mentioned populations are immersed.

On the other hand, the populations RA4 (middle zone) of P. douglasiana, TZ8 (low zone) of P. oocarpa and Q2 (upper zone) of P. pseudostrobus, presented the lowest levels of genetic variation, and for two of them the highest levels of inbreeding were obtained (RA4 and TZ8). Demographically, for the population of P. douglasiana in the upper part of the RA4 Reserve, little regeneration was observed and the lowest number of adult trees (<15), probably as a consequence of the frequent fires that have occurred in the area. The TZ8 population of P. oocarpa is also represented by few adult trees (<30), with no regeneration; in this case, it is mainly due to urban settlements, since the site is immersed in the southeastern region of Uruapan city. It is necessary to carry out restoration and repopulation activities based on the introduction of plants obtained from trees of P. pseudostrobus, P. douglasiana and P. oocarpa native to the region, to increase their census size and reduce the effect of inbreeding in the medium term.

In general, conservation and restoration strategies should be implemented focused on increasing the size of the populations, so that genetic diversity is maintained and with it, the biological, evolutionary and ecological representation of the pine species of the temperate forests of the basin of the Cupatitzio river.