Introduction
Genetic data on threatened plant populations can facilitate the development of adequate conservation strategies to reduce extinction risk (Hedrick, Reference Hedrick2001). This is particularly important for species affected by habitat fragmentation, which reduces the number of individuals per population and leads to isolation between populations, raising the probability of local extinction (Heinken & Weber, Reference Heinken and Weber2013).
The genus Magnolia is represented in Cuba by seven endemic taxa (Palmarola et al., Reference Palmarola, Romanov, Ch. Bobrov and González-Torres2016). The Critically Endangered subspecies Magnolia cubensis subsp. acunae Imkhan. is endemic to the Guamuhaya mountain range (González-Torres et al., Reference González-Torres, Palmarola, Bécquer, Berazaín, Barrios and Gómez2013, Reference González-Torres, Palmarola, González-Oliva, Bécquer, Testé and Barrios2016) and is threatened by deforestation and land conversion for cattle farming and coffee production. Recent studies on the distribution and conservation status of the subspecies’ populations (Palmarola et al., Reference Palmarola, González-Torres, Barrios, Albelo and León2012; Granado, Reference Granado2015) reported that the two main subpopulations have 416 and 70 individuals, respectively. However, there is no information on genetic diversity, the degree of habitat fragmentation or the interaction between these factors. Our study aimed to evaluate the effect of habitat fragmentation on the structure and genetic diversity of the M. cubensis subsp. acunae subpopulations in the Guamuhaya mountain range, to support effective conservation management of this subspecies.
Study area
The study was carried out in the 201.35 km2 Protected Natural Landscape Topes de Collantes and the 60.91 km2 Ecological Reserve Lomas de Banao, in the Guamuhaya mountain range in central Cuba (Fig. 1; CNAP, 2014). Each of these protected areas contains one of the two main subpopulations of M. cubensis subsp. acunae (Granado, Reference Granado2015).
Fig. 1 The two study areas in Guamuhaya, Cuba.
The climate of Guamuhaya corresponds to the Western Caribbean subregion and is classified as humid tropical, although at high elevations it could be considered mild warm (Dominguez & Acosta, Reference Domínguez, Acosta, Domínguez, Torres-Martínez and Puerta2012). The variety of ecosystems in the region contain high levels of biodiversity, which, particularly with respect to the flora, makes the Guamuhaya one of the most biodiverse and endemics-rich localities in Cuba (Ruiz et al., Reference Ruiz, Naranjo, Albelos, Rodríguez, Cruz and Duardo2011).
Methods
Fragmentation analyses
We used a Landsat 8 satellite image of Guamuhaya mountain range, taken in February 2014, with 30 m spatial resolution and 11 spectral bands (Roy et al., Reference Roy, Wulder, Loveland, Woodcock, Allen and Anderson2014), and cropped the areas of interest using a map of protected areas (CNAP, 2014). We extracted georeferenced points of the habitat types in which M. cubensis subsp. acunae occurs from the vegetation map of Estrada et al. (Reference Estrada, Martín, Martínez, Rodríguez, Capote, Reyes, Galano, Cabrera, Martínez, Mateo, Guerra, Batte and Coya de la Fuente2012). The points were used for supervised categorization of satellite images using the maximum likelihood method. In this process, pixels with a known land-cover type, located within the training areas, are used to categorize pixels of unknown land-cover type. We used four land-cover categories: submontane rainforest, montane rainforest, water and matrix, the latter defined as non-forested areas, including other vegetation, and agricultural and urban areas.
We converted the raster land-cover map to vector format and calculated seven fragmentation indices (Table 1). The variety of habitat patch shapes was classified according to Henao (Reference Henao1988) for the mean shape index and Hargis et al. (Reference Hargis, Bissonette and David1998) for the mean patch fractal dimension index. We used IDRISI Selva 12.1 (Eastman, Reference Eastman2012) for the supervised categorization and the Patch Analyst extension of ArcGis 9.3 (Esri, Redlands, USA) for the fragmentation analyses.
Table 1 Fragmentation indices used for the characterization of submontane and montane rainforest in the Protected Natural Landscape Topes de Collantes and the Ecological Reserve Lomas de Banao, two protected areas in the Guamuhaya mountain range, Cuba (Fig. 1).
Population genetic analyses
Sampling
We collected samples from 67 M. cubensis subsp. acunae individuals and stored them in self-sealed bags with silica gel. Of these, 58 were from Topes de Collantes (39 leaf samples from adult plants, 19 leaf samples from juveniles obtained from seeds) and 9 from Lomas de Banao (all leaf samples from adult plants), representing 10% and 13% of all individuals in these subpopulations, respectively. We collected seeds from fruits randomly selected from the Topes de Collantes subpopulation, using only one fruit per tree and one seed per fruit. The seeds were planted in nurseries in December 2014 and after two months, when the seedlings had six leaves, one leaf was collected for genetic analysis. Here we define population as the total number of individuals of the subspecies and subpopulation as the different population nuclei within the distribution area, according to IUCN (2001).
DNA extraction and genotypification
We extracted DNA from dried leaf tissue using a modified cetyltrimethylammonium bromide (CTAB) extraction protocol (Doyle & Doyle, Reference Doyle and Doyle1990), with MagAttract Suspension G solution mediated cleaning (Xin & Chen, Reference Xin and Chen2012). We genotyped individuals with 11 microsatellite markers or simple sequence repeats (Table 3) developed on four Neotropical Magnolia species: M. lacandonica (MA39), M. mayae (MA40), M. dealbata (MA41) and M. cubensis subsp. acunae (MA42) (Veltjen et al., Reference Veltjen, Asselman, Hernández Rodríguez, Palmarola Bejerano, Testé Lozano and González Torres2018), using the three-primer PCR multiplex method. PCR conditions and primer labelling followed Veltjen et al. (Reference Veltjen, Asselman, Hernández Rodríguez, Palmarola Bejerano, Testé Lozano and González Torres2018). Fragment analyses were executed by Macrogen Inc. (Seoul, South Korea) and we analysed the results in Geneious 8.0.5 (Kearse et al., Reference Kearse, Moir, Wilson, Stones-Havas, Cheung and Sturrock2012), using the microsatellite plugin.
Simple sequence repeats marker testing
We calculated the deviation from Hardy–Weinberg equilibrium, linkage disequilibrium and the inbreeding coefficient per locus (F IS) cf. Weir & Cockerham (Reference Weir and Cockerham1984) for each locus and for each subpopulation, using 10,000 dememorization steps, 100 batches and 5,000 iterations per batch in Genepop 4.3 (Rousset, Reference Rousset2008). We calculated deviations of both the uncorrected (Waples, Reference Waples2015) and (sequential Bonferroni) corrected p-values to the nominal level of α = 0.05 for both analyses.
Population structure and genetic diversity
We applied the classification method according to Manel et al. (Reference Manel, Gaggiotti and Waples2005) and performed runs in Structure 2.3 (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000) under the following conditions: 10,000 Markov chain Monte Carlo replicates after an initial burn-in of 10,000, using correlated allelic frequencies and assuming the admixture model. To obtain probability values of allocation of individuals to each genomic group, we used five repetitions for each set K (number of groups), with K set to run from 1 to 8. We determined the most probable number of groups from the value of ΔK obtained according to the method of Evanno et al. (Reference Evanno, Regnaut and Goudet2005). The results were visualized in Structure Harvester Web 0.6.94 (Earl & von Holdt, Reference Earl and von Holdt2012). An individual was considered to be a member of a genetic group when its probability of belonging to that group was > 0.5.
We quantified genetic diversity within each subpopulation and within maturity classes (i.e. adult or juvenile individuals) using the number of alleles per locus (A), number of private alleles (A p), allelic richness (A R), observed heterozygosity (H o), expected heterozygosity (H e), percentage of polymorphic loci (P) and inbreeding coefficient (F IS). We estimated genetic differentiation between subpopulations and across generations through pairwise comparisons of F ST values (Weir & Cockerham, Reference Weir and Cockerham1984) and interpreted these following the criteria of Hartl & Clark (Reference Hartl and Clark1997). Measures of genetic diversity were calculated with GenAlEx 6.1 (Peakall & Smouse, Reference Peakall and Smouse2012). We calculated F ST values, allelic richness and significant deviations from zero for F IS values with FSTAT 2.9.3.2 (Goudet, Reference Goudet2001).
We estimated genetic distances between individuals and subpopulations and carried out a principal coordinates analysis from the matrix obtained. To identify possible patterns of isolation by distance in the genetic differentiation of the studied subpopulations, we performed Mantel correlation tests with 10,000 permutations between genetic distances and geographical distances for all pairs of adult individuals of both subpopulations. We used GenAlEx for both analyses.
To investigate the occurrence of any bottlenecks in the sampled subpopulations, we characterized allele frequency distribution by locus, and evaluated deficit or excess of heterozygotes for the infinite alleles model, the two phase model and the stepwise mutation model, using the Wilcoxon test in Bottleneck 1.2.02 (Cornuet & Luikart, Reference Cornuet and Luikart1996).
Results
Fragmentation analyses
The supervised categorization projected on the distribution map of the habitat of M. cubensis subsp. acunae in Topes de Collantes and Lomas de Banao showed submontane rainforest covered a larger area than montane rainforest (Fig. 2). Lomas de Banao had greater landscape homogeneity and a smaller area categorized as matrix (other vegetation, and agricultural and urban areas). In Topes de Collantes the matrix land-cover type was primarily in the south and to a lesser extent at higher altitude, in the central part of the protected area. In Lomas de Banao, matrix areas were mainly in the submontane rainforest regions, whereas the montane rainforest was less fragmented.
Fig. 2 Land-cover types in (a) Protected Natural Landscape Topes de Collantes, (b) Ecological Reserve Lomas de Banao, Cuba, resulting from satellite Landsat image classification (Roy et al., Reference Roy, Wulder, Loveland, Woodcock, Allen and Anderson2014). Matrix refers to other vegetation, and agricultural and urban areas.
The number of patches in the four Subpopulation × Class categories was 107–169 (Table 2); in Topes de Collantes patches of both forest types were smaller and less variable in their size compared to Lomas de Banao.
Table 2 Fragmentation indices of submontane and montane rainforest in the Protected Natural Landscape Topes de Collantes and the Ecological Reserve Lomas de Banao in the Guamuhaya mountain range, Cuba.
Table 3 Locus names, repeat motifs and measures of genetic diversity of the 11 simple sequence repeat markers designed for the genetic characterization of Magnolia cubensis subsp. acunae in the Topes de Collantes and Lomas de Banao subpopulations.
1A, average number of alleles per locus; A R, allelic richness; H o, observed heterozygosity; H e, expected heterozygosity; P, exact probability of Hardy–Weinberg equilibrium test; F IS (W&C), F IS values of Weir & Cockerham (Reference Weir and Cockerham1984) per locus; *p < 0.005 Bonferroni corrected probability were considered statistically significant; NI, non-informative comparison because it is a monomorphic locus or presents low values of H o and H e.
The mean shape index indicated that for the submontane rainforest, the most common patch form was rectangular-oblong, whereas for the montane rainforest, most patches could be considered amorphous. The mean path fractal dimension for both forest types showed complex forms analogous to fractal objects (Table 2).
Population genetic analyses
Simple sequence repeats marker
We found no significant deviations from Hardy–Weinberg equilibrium for 10 out of the 11 tested loci. Only MA42_279 in Topes de Collantes significantly deviated from Hardy–Weinberg proportions. The linkage analysis showed that the loci were not in allelic association, with the exception of MA40_045 and MA42_279; MA42_279 was discarded from all subsequent analyses.
Of the 11 microsatellite loci, 10 were polymorphic for both subpopulations; MA41_076 was monomorphic in Lomas de Banao. We found 68 alleles among all loci (1–11 per locus), with a mean of 6.18 alleles per locus. Of the 24 rare alleles found, 22 (91.67%) were private alleles of Topes de Collantes. Detailed information on the genetic diversity estimates per locus is given in Supplementary Table 1.
Population structure and genetic diversity
The optimal value for K was 2 (Fig. 3). Of the adult trees, 75% were assigned to a genetic group that aligned with their sampling location. All individuals of Lomas de Banao clustered in genetic group 2, together with 30.77% of the adult individuals and 100% of the juveniles of Topes de Collantes. This means 12 of 39 adult individuals from Topes de Collantes had a probability of > 50% to belong to genetic group 2, which aligns with the individuals of Lomas de Banao. This analysis did also not distinguish between juveniles from Topes de Collantes and adults from Lomas de Banao within the second genetic group (Fig. 4).
Fig. 3 ΔK values obtained with K (number of groups) from K = 2 to K = 8, and five simulations per analysis, visualized with Structure Harvester Web 0.6.94 (Earl & von Holdt, Reference Earl and von Holdt2012), to determine the most probable number of groups of Magnolia cubensis subsp. acunae individuals from the Topes de Collantes and Lomas de Banao subpopulations, according to the method of Evanno et al. (Reference Evanno, Regnaut and Goudet2005).
Fig. 4 Probability of genetic group allocation of M. cubensis subsp. acunae individuals of Topes de Collantes and Lomas de Banao subpopulations, inferred at K = 2 based on allelic frequencies on simple sequence repeats data.
Values for all diversity measures were lower in Lomas de Banao than Topes de Collantes. Similarly, juveniles of Topes de Collantes were less diverse than the adult population of this locality. The inbreeding coefficient (F IS) was higher, and statistically significant, in Lomas de Banao compared to Topes de Collantes (Table 4).
Table 4 Measures of genetic diversity in adult and juvenile individuals of M. cubensis subsp. acunae in the Topes de Collantes and Lomas de Banao subpopulations.
1A, mean number of alleles; A R, allelic richness; H o, mean observed heterozygosity; H e, mean expected heterozygosity; HW, number of loci deviated from the proportions of Hardy–Weinberg; F IS, population inbreeding coefficient, significant deviations from zero are indicated with * (p = 0.05); P, percentage of polymorphic loci.
No genetic differentiation was apparent between the two subpopulations using the F ST fixation index (Table 5). The principal coordinates analysis showed no differential clustering for individuals from different subpopulations or maturity classes (Fig. 5). The Mantel correlation test had a low correlation coefficient (r = 0.034), considered statistically non-significant (p = 0.32), indicating no effect of distance on genetic variation.
Fig. 5 Principal coordinates analysis of genetic distances matrix with 10 microsatellite markers in M. cubensis subsp. acunae in Guamuhaya, Cuba.
Table 5 Fixation index F ST (above diagonal) and F IT (below diagonal) values, determined from 10 microsatellite loci, for all paired comparisons, of sampled subpopulations of M. cubensis subsp. acunae.
The Wilcoxon test for Topes de Collantes showed differences between the estimated heterozygosity value at equilibrium and that obtained by simulating diversity under the infinite alleles model. For Lomas de Banao, the analysis also revealed differences for the two stage mutation model and stepwise mutation model (Table 6). The distribution of allelic frequencies displays an L-shaped distribution (with a high proportion of low frequency alleles and a few frequency alleles near 1) for Topes de Collantes, which is not the case for Lomas de Banao (Fig. 6).
Fig. 6 Allele percentage distribution by frequency ranges for 48 adult individuals of subpopulations of M. cubensis subsp. acunae in Guamuhaya, Cuba.
Table 6 Probability of deficiency and excess heterozygosity (H e) observed at equilibrium, compared to that estimated under the infinite alleles model, two phase model and stepwise mutation model for the subpopulations of M. cubensis subsp. acunae.
*p < 0.01 was considered statistically significant.
Discussion
Fragmentation analyses
Habitat fragmentation was most severe in the montane rainforest of Topes de Collantes, with smaller, irregularly shaped patches and greater patch perimeter to area ratio. Edge effects and the quality of the surrounding area also influence the characteristics of remnant habitat patches (Heinken, & Weber, Reference Heinken and Weber2013). A matrix with a similar structure to a remnant forest fragment will have less influence than one with a different structure (Fischer & Lindenmayer, Reference Fischer and Lindenmayer2007). In the montane rainforest, where many fragments are surrounded by submontane rainforest, the matrix effect is thus smaller, whereas the submontane rainforest is mainly surrounded by matrix areas.
The high degree of forest fragmentation was expected, given the history of land-use changes in the area. Coffee has been produced in Topes de Collantes since the 19th century, resulting in selective logging across large areas and the replacement of natural vegetation by a monoculture. Furthermore, the wood of native tree species such as M. cubensis subsp. acunae is also used to build coffee plantation infrastructure (Domínguez et al., Reference Domínguez, Torres Martínez and Puerta2012).
Population genetic analyses
Simple sequence repeats markers
Deviations from the Hardy–Weinberg equilibrium found in the MA42_279 locus in both adults and juveniles of Topes de Collantes can be explained by the presence of inbreeding. F IS was significantly larger than zero, which indicates a deficit of heterozygotes. Although we cannot discard this phenomenon, the presence of the so-called Wahlund effect (i.e. a reduction of heterozygosity in a population caused by subpopulation structure) should affect the loci equally (Waples, Reference Waples2015). It therefore seems most plausible that apparent inbreeding at one of 11 loci is a chance effect, rather than a result of actual inbreeding, which should be evident across most loci. Another possible interpretation is the presence of null alleles or errors in genotyping (Pompanon et al., Reference Pompanon, Bonin, Bellemain and Taberlet2005). However, in our study the markers that presented such errors were eliminated from the final selection (Veltjen et al., Reference Veltjen, Asselman, Hernández Rodríguez, Palmarola Bejerano, Testé Lozano and González Torres2018).
Population structure and diversity
The analysis of population structure showed two genetic groups that partly corresponded with the two adult demographic groups defined prior to the analysis, and diversity analyses showed absence of differentiation between subpopulations. The subpopulation in Topes de Collantes was more diverse, whereas that in Lomas de Banao showed lower values for all diversity indices. This could be a result of the different subpopulation sample sizes. This is supported by the fact that 50% of the loci used in our study showed higher values for genetic diversity indices with samples from 39 adult plants from Topes de Collantes, compared to the analyses carried out by Veltjen et al. (Reference Veltjen, Asselman, Hernández Rodríguez, Palmarola Bejerano, Testé Lozano and González Torres2018) on 20 individuals. The Topes de Collantes subpopulation has the highest number of individuals, with 75% of the global population (Granado, Reference Granado2015), whereas Lomas de Banao includes 15% of the global population. The sample size was thus larger for Topes de Collantes compared to Lomas de Banao, but both sample sizes are considered comprehensive and cover a similar proportion of the known individuals per subpopulation. However, allelic richness by locus for a fixed sample size (Allendorf et al., Reference Allendorf, Luikart and Aitken2013), was also higher for Topes de Collantes in most loci evaluated.
The diminution of the diversity across generations in Topes de Collantes may reflect a decrease of pollen and seed dispersal. Pollen dispersal distances for trees visited by small insects such as beetles in closed-canopy forests often do not exceed 300 m (Dick et al., Reference Dick, Hardy, Jones and Petit2008), and it is possible that the loss of one genetic group in the juveniles of Topes de Collantes compared to the adult trees of this subpopulation is a result of reduced pollen movement. However, given the mix of c. two-thirds genetic group 1 and one-third genetic group 2 in the adult generation, a similar representation was expected in the juveniles. There are three potential explanations for this:
(1) Gene flow between subpopulations: the geographical proximity of the subpopulations increases the gene flow mediated by fruit dispersal and tends to homogenize the present genetic variation. In a fragmented habitat, the dispersion is affected by both the isolation of patches and the movement of the dispersers (Laaka-Lindberg et al., Reference Laaka-Lindberg, Korpelainen and Pohjamo2003). The potential dispersers of this subspecies, the Cuban trogon Priotelus temnurus, fieldfare Turdus pilaris and western spindalis Spindalis zena, are permanent or seasonal residents in Cuba (Garrido & Kirkconnell, Reference Garrido and Kirkconnell2011) and could move between the two subpopulations, which are 33 km apart. However, detailed studies of the ecology, behaviour and foraging strategies of potential seed dispersers in this mountainous region are required to elucidate the flow of diaspores between the subpopulations. It is unlikely that pollen moves between the subpopulation because there are only few reports of pollen travelling distances > 10 km for insect-pollinated species (Petit & Hampe, Reference Petit and Hampe2006).
(2) Evolutionary history: the simplest interpretation of small genetic distances between subpopulations is that they share a recent common ancestor. The homogeneity in Lomas de Banao could be an indication of the founder effect (Slatkin, Reference Slatkin2004); i.e. a few, genetically similar individuals from Topes de Collantes could have founded this subpopulation. However, it is more likely that both subpopulations once belonged to a single, larger population.
(3) The species’ crossing system: cross-pollination is one of the mating systems of Magnolia (Thien, Reference Thien1974), guaranteeing that a few migrants per generation are sufficient to counter genetic differentiation (Holsinger & Weir, Reference Holsinger and Weir2009). These species are thus generally characterized by low genetic differentiation between populations. This could explain the genetic similarity between the two subpopulations of M. cubensis subsp. acunae.
Topes de Collantes tested positive for a past bottleneck, assuming the infinite alleles model. However, for microsatellite markers, the stepwise mutation model is considered more appropriate (Putman & Carbone, Reference Putman and Carbone2014), which would indicate no trace of a past bottleneck for this subpopulation. Furthermore, the distribution of allelic frequencies in Topes de Collantes is typical for a population that has not suffered a recent bottleneck (Allendorf et al., Reference Allendorf, Luikart and Aitken2013). In Lomas de Banao the existence of a bottleneck is possible considering the low number of individuals, the low genetic diversity and the high values of inbreeding. However, the results should be interpreted with caution because the Wilcoxon test requires at least 10 individuals (Cornuet & Luikart, Reference Cornuet and Luikart1996) and could be compromised by the small sample size from Lomas de Banao (9 individuals).
Implications for conservation
The principal consequences of ongoing habitat fragmentation are progressive reduction of subpopulation sizes and increased distance between habitat fragments, which can affect a population's genetic variation, heterozygosity, inbreeding, gene flow and genetic divergence between subpopulations (Heinken & Weber, Reference Heinken and Weber2013). These potential genetic consequences of fragmented habitats are reflected in our results, with high levels of inbreeding and low genetic diversity in Lomas de Banao, reduction of the genetic diversity across generations in Topes de Collantes and low differentiation between subpopulations. Our findings merit a new approach to the conservation and management of M. cubensis subsp. acunae.
The Lomas de Banao subpopulation should be considered a priority for the population reinforcement actions currently being carried out in the Guamuhaya region. The results of the structure analyses, principal coordinates analysis plots and F ST values predict there is no marked genetic structure in the subpopulations. Our findings suggest the two subpopulations can be considered a single evolutionary unit and conservation entity (Moritz, Reference Moritz1994). Granado (Reference Granado2015) proposed nine areas for Magnolia population reinforcement in Guamuhaya, taking into account the species’ modelled (potential) distribution, the area's protection status and landscape use. From these analyses, and considering the similarity between the subpopulations, we suggest using individuals from both subpopulations for such reinforcement actions.
Despite past habitat fragmentation and loss of natural vegetation cover in these protected areas, the improvement of the legal framework on biological diversity in Cuba in recent years, which prohibit agricultural, forestry and fruit production in protected areas and set goals to reduce the impact of agroforestry, will support the recovery of the studied subpopulations and the ecosystems to which they belong.
Historically, there was little interest in conserving M. cubensis subsp. acunae in Topes de Collantes, but after 12 years of conservation projects focussing on this subspecies, awareness has increased. People are now interested in protecting this subspecies and its habitat because it represents a symbol for the community and its cultural identity. Many farmers cultivate coffee in the shade of M. cubensis subsp. acunae and claim that coffee quality is superior when shade is provided by native species such as Magnolia. Restoration of the montane and submontane rainforest, population reinforcement and support for these matters amongst local communities are essential for the long-term improvement of the conservation status of this endemic Cuban subspecies.
Acknowledgements
We thank the PlantLife Conservation Society, Whitley Fund for Nature, Fauna & Flora International, Fondation Franklinia, Instituto de Ecología, A.C., Magnolia Society International, Ghent University, National Botanical Garden of Cuba and Cuban National Center of Protected Areas for their support. We are grateful to all who helped during field work.
Author contributions
Study design and fieldwork: all authors; data analysis: MH, PA, EV; writing: MH; revisions: all authors.
Conflicts of interest
None.
Ethical standards
This research abided by the Oryx guidelines on ethical standards. It was conducted in Cuba, the country of residence of the first author and most collaborators, and all legal requirements regarding permits were followed (permit No. 15/16 CICA-CITMA). Field researchers adopted the precautions to avoid the accidental introduction of invasive and pathogenic organisms and no specimens were killed for the purposes of this study.
