Genetic diversity in micro-endemic plants from highland grasslands in southern Brazil

Abstract

Population genetic structure results from the interaction between historical events, current ecological conditions and life traits. The genetic structure and gene flow between populations are important to species dynamics, mainly for rare and endangered species that are more vulnerable to landscape changes and fragmentation. Here we evaluated the genetic diversity, population structure and gene exchange in Petunia bonjardinensis, P. reitzii and P. saxicola, three rare species endemic to subtropical highland grasslands in southern South America. We analysed the genetic diversity and structure considering historical events, such as founder effect and climate changes, and biological traits of each species. We also estimated the conservation status for these three species. We collected samples from all adult individuals and occurrence sites that could be found at the same flowering season and genotyped them for 13 nuclear microsatellite markers. Our results indicate that rarity is probably historical for these species, given that we found no genetic evidence for recent bottlenecks. Petunia bonjardinensis, with the largest occurrence area and population sizes, displayed the higher diversity indices. The other two showed lower genetic diversity and are geographically most restricted. Gene exchange among these species was low, although they share some ancestral genetic polymorphism. Historical migration, founder effects and Pleistocene climate cycles ae the main factors explaining genetic diversity, and this was also influenced by reproductive biology and recent habitat loss, whereas the landscape influences the structure. Based on IUCN criteria, the three species are endangered, and the main risk for their survival is probably anthropic activity in the occurrence area. We recommend an urgent programme for the preservation of these species in situ and ex situ.

INTRODUCTION

On a global scale, plant endemism is more common than previously thought (Enquist et al., 2019). Most ecologically stable areas (e.g. the Andes and Brazilian Atlantic Forest) harbour high concentration of species with small range sizes (Morueta-Holme et al., 2013). Rare plants are especially susceptible to global climate changes and anthropogenic pressures (Enquist et al., 2019; Vincent et al., 2020), due to the natural selection or demographic phenomena can be more pronounced in these populations as they usually have reduced levels of genetic diversity (Ellstrand & Elam, 1993) and great risk of local extinction (Blanco-Pastor, Fernández-Mazuecos & Vargas, 2013).

A common characteristic among endemic and narrowly distributed species is low levels of genetic variation that could be the cause or a consequence of their rarity (e.g. Hamrick & Godt, 1989; Gitzendanner & Soltis, 2000; Peakall et al., 2003). A small effective population size reduces genetic diversity through allele loss due to genetic drift and inbreeding (Ellstrand & Elam, 1993). Landscape fragmentation may further increase genetic erosion (Aguilar et al., 2008). Conversely, there are some examples of high genetic diversity in narrowly distributed plants (e.g. Song & Mitchell-Olds, 2007; Turchetto et al., 2016; Rodrigues et al., 2019), probably due to their evolutionary histories and reproductive strategies.

The open field areas in the subtropical highlands of southern South America (Fig. 1A), as part of the Brazilian Atlantic Forest ecoregion, hold high endemism and plant diversity (Iganci et al., 2011). Remarkably, the southern Brazilian subtropical highland grasslands (SHG) harbour high genetic diversity associated with niche stability (Barros et al., 2015; 2020). In SHG, mosaics between grasslands and araucaria forests dominate the natural vegetation, a landscape in which the Brazilian pine, Araucaria angustifolia (Bertol.) Kuntze, is the dominant tree species (Safford, 1999; Behling, 2002; Andrade et al., 2016). The area is an escarpment that rises c. 1000 m from the Atlantic coastal plain and slowly decreases in elevation to the floodplains in the west (Paraná and Uruguay basins). Hills, ridges, plateaus and valleys characterize the landscape, where the climate is humid with mean annual precipitation is 1500–2000 mm, mean annual temperatures between 15 and 20 °C and frequent snowfalls during the winter. The soils have a basaltic origin that cover the effusive rocks, which originated c. 120–135 Mya (Plá et al., 2020).

Many pollen-based records have elucidated grasslands and forest dynamics during the Pleistocene (Behling & Pillar, 2007) and revealed the role of past climate changes on species diversification in some plant groups (e.g. Lorenz-Lemke et al., 2010; Fregonezi et al., 2013; John et al., 2019). Glacial and interglacial cycles promoted grassland and forest expansion–contraction, that should have fragmented the SHG and facilitated several allopatric speciation events (Boucher, Zimmermann & Conti, 2016). Currently, the SHG have suffered severe habitat loss due to the human activities, being replaced by large mechanized plants and forestry (Plá et al., 2020), cattle ranching and natural forage fire management (Andrade et al., 2016) that are responsible by the landscape fragmentation (Supporting Information, Fig. S1). Only a few conservation areas can be found in SHG (Supporting Information, Fig. S2).

Petunia Juss. (Solanaceae) is a charismatic genus comprising 14 currently recognized species and the commercial hybrid Petunia × hybrida (Hook.) Vilm. (Stehmann et al., 2009). The genus is closely related to Calibrachoa Cerv. (also tribe Petunieae) from which it diverged c. 8.5 Mya (Särkinen et al., 2013). A molecular phylogenetic and biogeographical study (Reck-Kortmann et al., 2014) identified two main clades, strongly supported by the corolla tube length, and indicated that the genus originated in lowland grasslands, specifically in the Pampa region, with later migration to grasslands in higher elevations. In the short corolla tube clade, species are bee-pollinated, usually with bright purple corollas and violaceous pollen. Seven of these species displaying short corolla tube have diversified in highland grasslands, under the influence of late Pleistocene climate changes that promoted expansion and contraction of the distribution of their ancestor leading to diversification and speciation in allopatry (Lorenz-Lemke et al., 2010).

Among the Petunia spp. that grow in the SHG, P. bonjardinensis T.Ando & Hashim., P. reitzii L.B.Sm. & Downs and P. saxicola L.B.Sm. & Downs (Fig. 1B) are endemic and narrowly distributed species, for which only a few sites and individuals are known (Fig. 1A), and all of them are included in the short corolla tube clade. Although individuals of these species occur in different elevational ranges, they share several morphological traits (Stehmann et al., 2009) and have shallow genetic diversity as estimated based on plastid sequences (Lorenz-Lemke et al., 2010). Like other Petunia spp., these diploid, annual, herbaceous plants form small patches of just a few individuals that blossom synchronously from late October to early January. Concerning their reproductive biology, P. bonjardinensis displays floral characteristics usually associated with mellitophilly, the flowers displaying an exerted stigma positioned above the anthers (Ando & Hashimoto, 1993). Based on the bright corollas and other morphological floral traits, it has been suggested that P. reitzii and P. saxicola should be considered bird-pollinated, probably by hummingbirds that have been observed visiting and hovering in front of their flowers (Ando et al., 1999). A controlled experiment reported that these three species are self-incompatible (Robertson, Goldberg & Igić, 2011).

Recently, several studies have been conducted evaluating the genetic diversity in Petunia spp. based on nuclear microsatellites (Segatto et al., 2017; Turchetto et al., 2015a, b, 2016, 2019b; Silva-Arias et al., 2017; Backes et al., 2019; Giudicelli et al., 2019; Rodrigues et al., 2019; Schnitzler et al., 2020; Turchetto, Schnitzler & Freitas, 2019a), and many of these studies compared the genetic diversity based on nuclear polymorphism with the evolutionary history described from plastid analyses. These studies revealed different levels of diversity for species with restricted distribution and those with wide distribution that mostly agreed with phylogeographical results (Segatto et al., 2014; Ramos-Fregonezi et al., 2015; Turchetto et al., 2015a, b; 2016; Backes et al., 2019). The narrowly distributed species from the SHG, despite their phylogeographical history having been previously described (Lorenz-Lemke et al., 2010), have still not been characterized for nuclear diversity. Thus, the present study aimed to: (1) evaluate the genetic diversity and population structure of three species from SHG, P. bonjardinensis, P. reitzii and P. saxicola; (2) identify the conservation status of these species based on some International Union for the Conservation of Nature (IUCN) criteria and (3) evaluate the influence of reproductive systems, historical climate changes and distribution range on these species structure and diversity. We used data from nuclear microsatellites genotyped in all adult individuals found across the entire distribution range of the three species during a reproductive season.

MATERIAL AND METHODS

Sampling

We collected fresh and young leaves from 90 P. bonjardinensis individuals (from 12 collection sites), 23 P. reitzii individuals (from five collection sites) and 27 P. saxicola individuals from its only known location (Supporting Information, Table S1). No statistics that involved multiple population comparisons were estimated for P. saxicola. The samples were taken from all adult individuals found in the same flowering season (Fig. 1). The number of individuals reflects the current rarity of each species as we searched for the species in all recorded locations obtained from SpeciesLink (http://www.splik.org.br).

These species are not sympatric, inhabiting different micro-environments. Plants of P. bonjardinensis grow in open and disturbed areas, such as roadside slopes, at 1200–1500 m, individuals of P. reitzii occur at 800–1000 m and P. saxicola inhabits humid and shadowed rocks at 800 m in elevation. The main characteristic differentiating these species is the habit (decumbent in P. bonjardinensis, ascendent in P. reitzii and saxicolous in P. saxicola), and the species differ also in the stigma position (above the longest pair of anthers in P. bonjardinensis, below them in P. reitzii and slightly exserted in P. saxicola; Stehmann et al., 2009). Petunia reitzii and P. saxicola can be separated based on leaf shape (linear-lanceolate and elliptical, respectively; Ando & Hashimoto, 1993).

We took the geographical coordinates using the Global Positioning System (GPS) and one voucher per collection site. We deposited the vouchers at Universidade Federal de Minas Gerais herbarium—BHCB-UFMG, Belo Horizonte, Brazil (Supporting Information, Table S1).

DNA extraction and microsatellite genotyping

We extracted total genomic DNA from silica gel-dried leaves using a CTAB-based protocol (Roy et al., 1992). We estimated the DNA quality and quantity using a Nanodrop DN-1000 Spectrophotometer (Thermo Fisher Scientific Co., Waltham, USA) and a Qubit Fluorometer (Thermo Fisher), respectively. We screened 13 microsatellite loci previously described (Bossolini et al., 2011) for P. ×hybrida (Supporting Information, Table S2) following the protocol adapted for Petunia spp. (Turchetto et al., 2015b). We labelled the PCR products with the fluorescent dyes FAM, NED, PET or VIC (Thermo Fisher). We performed PCR amplifications for each SSR locus separately and multiplexed. They were denatured and size-fractionated the PCR products using capillary electrophoresis on an ABI 3100 genetic analyser (Thermo Fisher) with a LIZ (500) molecular size standard (Thermo Fisher). We used the Genemarker v.1.97 software (Softgenetics LLC, State College, PA, USA) to identify the alleles. Additionally, we visually verified and scored all alleles. All individuals exhibited a maximum of two alleles per locus, as expected for diploid species, and the sizes of the alleles were compatible with the repetition for each locus (Bossolini et al., 2011).

Genetic diversity

We used Micro-Checker (http://www.microchecker.hull.ac.uk/) to estimate genotyping errors due to stutter bands, allele dropout and null alleles. For each species, we estimated the number of alleles per locus (A), allele richness (AR), the number of private alleles (PA), Nei’s unbiased gene diversity (GD; Nei, 1987) and inbreeding coefficients (F_IS; Weir & Cockerham, 1984) with FSTAT v.2.9.4 (Goudet, 1995). We also computed the observed (H_O) and expected (H_E) heterozygosity under the Hardy–Weinberg equilibrium (HWE) after Bonferroni’s correction using Arlequin v.3.5.1.2 (Excoffier & Lischer, 2010,). The pairwise genetic differentiation between the collection sites for the three species based on F_ST was performed in Arlequin, and significance was tested with a permutation test. Differences in the distribution of genetic variation between collection sites were quantified through the pairwise R_ST estimated in Arlequin. We compared F_ST and R_ST to estimate the role of drift and mutation in the population differentiation based on Hardy et al. (2003). We compared the diversity indices among species through a Kruskal–Wallis test using the IBM SPSS Statistics v.24.0 (IBM Corp., Armonk, USA).

To verify the contribution of stepwise-like mutations to the genetic differentiation, we estimated Slatkin’s R_ST (Slatkin, 1995) obtained from an analysis of variance of allele size. Then, we tested the hypothesis that F_ST = R_ST based on Hardy et al. (2003), using the software SPAGeDI v.1.5 (Hardy & Vekemans, 2002). The comparison of F_ST and R_ST provides insights on the role of genetic drift and mutation on the population differentiation because R_ST is expected to be larger than F_ST under stepwise-like mutations and equal when solely drift causes the differentiation (Hardy et al., 2003). The comparisons were run for P. bonjardinensis and P. reitzii populations separately, considering all 18 populations and all individuals of each population forming a group.

We estimated the hierarchical analysis of molecular variance (Excoffier, Smouse & Quattro, 1992) in Arlequin. In the variance analysis, we used the length of the amplified product of SSR markers of each accession as the value of microsatellite alleles to evaluate the proportion of variation can be attributed to differences between species (F_SC), among populations per species (F_CT) and within populations (F_ST).

Genetic structure

We inferred the genetic structure within and between species using Structure v.2.3 (Pritchard, Stephens & Donnelly, 2000), with no previous information about taxonomy or geographical origin for each sample, and estimated the best K value with ∆ K (Evanno, Regnaut & Goudet, 2005) in Harvester v.0.6.93 (Earl & von Holdt, 2012) online version. We performed Structure runs using 10⁶ Markov chain Monte Carlo (MCMC) repetitions after 2.5 × 10⁵ burn-in period and ten iterations per K value, evaluating different numbers of clusters (1 to 7). The resulting bar plot from the summarized iterations for the best K was generated using Clumpak online (http://clumpak.tau.ac.il/contact.html). Genetic relationships among individuals of the three species were examined employing the discriminant analysis of principal components (DAPC) (Jombart, Devillard & Balloux, 2010, from Adegenet package; Jombart, 2008) on the SSR data set, with taxonomic priors.

We also analysed the spatial genetic structure (SGS) for P. bonjardinensis and P. reitzii using SPAGeDi (Hardy & Vekemans, 2002). We estimated the pairwise kinship coefficient (F_ij; Loiselle et al., 1995), regressed on the logarithm of the spatial distance between individuals in each population. We tested the significance through 10 000 permutations and the jackknife approach to obtain standard errors across loci. For graphical representation, the average pairwise F_ij values were computed against pairwise spatial distances. We defined the intervals of distance classes aiming to maximize the number of individual pairs in each class. The amount of SGS for each species was estimated using Sp statistics [$Sp=blog/(1F1)$⁠; Vekemans & Hardy, 2004].

Species relationships and demographic history

We estimated an unrooted neighbor-joining (Saitou & Nei, 1987) tree (NJ) from the D_C (Cavalli-Sforza & Edwards, 1967) pairwise genetic distance between individuals (within and between species) with Populations v.1.2.32 (Langella, 2021). The Dc distance is based on shared alleles assuming a stepwise mutation model (SMM) for microsatellite changes. Branch support values in the NJ were obtained through 1000 bootstraps. We visualized and edited the tree using FigTree v.1.4.2 (Rambaut, 2008). Dc genetic distances between species were estimated with the same method.

We estimated the evolutionary history of this species group using the population model selection in Migrate-N v.4.4.4 (Beerli, 2009) with the SSR data set. Seven evolutionary scenarios were tested (Fig. 2), including the three possible rooted trees between these species, with and without migration, and an island model, in which the three species exchange migrants without a shared history. Effective population sizes (N_e), effective migration rates (m) and divergence times between species were estimated with the Bayesian coalescent approach (Beerli & Felsenstein, 2001; Beerli, 2006). We used the Brownian motion model, with starting conditions based on F_ST and uniform prior distribution to estimate theta (θ) and M values. We showed the proposal and prior distribution for all parameters used in the analyses in Supporting Information, Table S3. We ran four long chains with three temperatures (1.0, 1.5 and 3.0), resulting in a total of 10 000 genealogies record every 50 steps, after the first 1000 have been discarded as burn-in. Estimates were computed per locus and summarized as weighted values across all loci per species. We used the constant mutation rate across all loci and uniform distribution to estimate population size and gene flow. To obtain N_es, we used a punctual SSR mean mutation rate (μ) of 10^–3, and the two extremes 10^–4 and 10^–2 (Zhang & Hewitt, 2003) with the formula N_e = θ/ 4 × μ. To estimate the number of migrants, we used M with the θ value of the recipient population (m = M × θ). The models were compared using Bayes factor tests based on Bezier-approximated marginal likelihoods (Beerli & Palczewski, 2010).

As landscape fragmentation is a reality throughout the geographical distribution of these species, we tested for evidence of population reduction in each species using Bottleneck v.1.2.02 (Piry, Luikart & Cornuet, 1999). Bottleneck assumes one population that has undergone a recent bottleneck has more heterozygotes than expected under the HWE. We tested the excess of heterozygotes using a Wilcoxon signed-rank test (Cornuet & Luikart, 1996) with a two-phase model (TPM) of microsatellite evolution (Di Rienzo et al., 1994), a variance of 12% (Piry et al., 1999), different percentages (70, 85 and 95%) of SMM (Di Rienzo et al., 1998) and 1000 replications.

We evaluated the recent gene flow between species employing a Bayesian multilocus genotyping procedure implemented in BayesAss v.3.0.4 (Wilson & Rannala, 2003). We ran the analysis with an MCMC chain length of 3 × 10⁷ iterations and 3 × 10⁶ iterations as burn-in. Runs were carried out in duplicate, and we assessed the convergence of the MCMC using Tracer v.1.7 (Rambaut et al., 2018) to validate data.

Conservation status

We adopted the criteria from the IUCN (2012; http://www.iucnredlist.org/) to identify the threat categories of P. bonjardinensis, P. reitzii and P. saxicola. Conservation status was based mainly on the distribution area, as estimated with GeoCAT (Geospatial Conservation Assessment Tool) software (Bachman et al., 2011). The geographical distribution area for each species was estimated from the value of two parameters (extent of occurrence and the area of occupancy). The conservation status followed the nomenclature of IUCN: Least Concern (LC), Not Threatened (NT), Vulnerable (VU), Endangered (EN) and Critically Endangered (CR). We used all records obtained from SpeciesLink and Global Biodiversity Information Facility (http://www.gbif.org). Records obtained from the databases were verified manually for incongruence, and all duplicates were removed. Only records with GPS coordinates and detailed localization were used. We estimated the extent of occurrence and the area of occupancy based on total records and records obtained from 2000 to the present due to the expressive increase in silvicultural plantations in SHG since then (Hermann et al., 2016).

RESULTS

Genetic diversity and differentiation

For the three species, all pairs of loci were in linkage equilibrium (P < 0.001, Bonferroni’s adjusted value at α = 0.05), and all loci, except PM195 and PM8, were polymorphic for the three species (Supporting Information, Table S4). Locus PM195 was monomorphic for P. reitzii and P. saxicola, and PM8 was monomorphic for P. saxicola only. The frequency of null alleles was low (< 0.5%) across all loci and species, and we did not found evidence of genotyping errors.

In terms of genetic parameters (Table 1), P. bonjardinensis showed the highest values, especially in absolute indices (AL, AR and PA), as expected, given the sample size and number of sites studied. At the same time, P. reitzii and P. saxicola were equally diverse, considering the same indicators. The genetic diversity (GD) values were similar for all species. The three species were in HWE after Bonferroni’s correction, and inbreeding coefficients (F_IS) were not significant. Similarly, population differentiation indices (overall F_ST and R_ST) were higher in P. bonjardinensis than in P. reitzii. In some cases, differentiation was higher between populations of the same species than compared with populations of other species. The comparison between species revealed F_ST significant values (P. bonjardinensis vs. P. reitzii and P. reitzii vs. P. saxicola: F_ST = 0.09; P. bonjardinensis vs. P. saxicola: F_ST = 0.05; P < 0.05), whereas only the comparison between P. reitzii and P. saxicola had significant values based on R_ST (P. bonjardinensis vs. P. reitzii and P. bonjardinensis vs. P. saxicola: R_ST = 0.03; P. reitzii vs. P. saxicola: R_ST = 0.12; P < 0.05). For pairwise F_ST and R_ST values between the collection sites for the three species, see Supporting Information, Table S5.

Slatkin’s R_ST indicated differentiation among populations in all comparisons. Considering the species individually, P. bonjardinensis populations produced R_ST = 0.068 (SE = 0.014, P < 0.001) that did not differ from F_ST (P = 0.477), and P. reitzii had R_ST = 0.015 (SE = 0.011, P = 0.554) that did also not differ from F_ST (P = 0.498). Similarly, in the comparison of the total sampling, all the 18 populations, R_ST = 0.052 (SE = 0.011, P < 0.001) that did not differ from the overall F_ST (P = 0.782), and considering the three groups as individuals, R_ST = 0.032 (SE = 0.006, P = 0.0045) that did not differ from F_ST (P = 0.982).

The hierarchical AMOVA showed that most variation is found within populations (87%), with F_ST = 0.87, F_SC = 0.09 and F_CT = 0.04 observed in the comparison among species, reflecting some polymorphism sharing or gene flow among these three Petunia spp.

Genetic structure

We did not observe any within species structure based on collection sites, with all individuals in each species sharing the same genetic component. Considering all samples, the uppermost level of genetic structure obtained with ∆ K method from Structure results were K = 2 and K = 3 with similar likelihood. At K = 2 (Fig. 3A), P. reitzii and P. saxicola shared the same genetic component, whereas P. bonjardinensis formed another genetic component. At K = 3 (Fig. 3B), P. saxicola individuals shared a single component that is also presented in P. reitzii individuals in a proportion of c. 0.5. Petunia bonjardinensis individuals presented a high proportion of an exclusive component and a minor proportion of another component shared with P. reitzii individuals. At K = 5 (Supporting Information, Fig. S3), P. reitzii individuals presented a single exclusive component, and P. bonjardinensis individuals presented multiple genetic components. The DAPC analysis separated individuals by their species of origin, with a single P. saxicola individual clustered with P. bonjardinensis individuals, whereas a few others were found in intermediate positions (Fig. 4A). If we consider the coefficients of membership (q) of the DAPC result ≥0.9 as purebred individuals, only one individual from P. reitzii and one from P. saxicola showed a mixed membership with a higher proportion of the genetic component attributed to P. bonjardinensis (Fig. 4B). Finally, one P. saxicola individual was assigned as purebred P. bonjardinensis. These were the same individuals noted by their positions superimposed with individuals of other species in the DAPC graph.

The regression of F_ij on the logarithm of geographic distances was significant only for the first distance class in P. bonjardinensis (Fig. 5). Both P. bonjardinensis and P. reitzii had small values for the Sp statistics (0.004 and 0.002, respectively), suggesting the absence of SGS.

Species evolutionary relationships

The average intraspecific pairwise genetic distances (Dc) between individuals were similar for the species and lower than the interspecific distances (Table 2). Petunia reitzii and P. saxicola are closer to each other than to P. bonjardinensis. The neighbor-joining tree (Supporting Information, Fig. S4) shown individuals from all species intermixing across the tree. However, several small homogeneous groups were also found, but in most of these cases individuals were from different collection sites. All branches had low support values (BS ≤ 50).

Demographic history

We detected a heterozygote deficit with Bottleneck for the three species, independently of the SSR evolutionary model were used (Supporting Information, Table S6). However, these results should be interpreted with caution as < 20 SSR loci were analysed.

Of the seven demographic models tested in Migrate-N (Fig. 2), model A, in which P. reitzii and P. saxicola are sister species, and the three extant species exchange genes presented a probability of 1.0 (Table 3; Supporting Information, Fig. S5). The posterior distributions of the mutation-scaled demographic parameters from model A were shown in Table 4. Population sizes and migration rates were well estimated, but divergence times were not well estimated in any scenario (Supporting Information, Fig. S5). The N_e for P. bonjardinensis was significantly larger than the population sizes of the other two species, which were similar. Using μ of 10^–3, the mode and 95% confidence interval for the N_es for each species were: P. bonjardinensis: 1013 (817–1317); P. reitzii: 588 (408–775) and P. saxicola: 496 (342–667). The number of immigrants per generation were c. 1.0 when we used the mode estimation of θ (Table 5).

The recent gene exchange values between species, estimated by BayesAss as the fraction of the individuals in a population that are migrants from other population, were all small and statistically similar, although the mean was higher from P. bonjardinensis to the other species (Supporting Information, Table S7).

Conservation

We estimated the extent of occurrence and area of occupancy for the three Petunia spp. considering the total records and only the records from the year 2000 to the present separately. We evaluated the risk status for each species according to these IUCN criteria (Table 6). We observed a sharp decrease in the distribution area of all three species based on the two analysed criteria, and this was more pronounced in P. saxicola. The three species were classified as EN (endangered) according to these IUCN criteria, with P. saxicola becoming CR (critically endangered) in the last two decades considering the extent of occurrence.

DISCUSSION

We report the GD and conservation status of three endemic and rare Petunia spp. found in SHG based on 13 microsatellite loci. Genetic diversity provides useful information on the evolutionary history and relationships between species and microsatellites are especially informative for revealing the GD of rare species (e.g. Kalia et al., 2011).

The three Petunia spp. from SHG are currently distributed in small areas, fragmented by human activities such as urbanization and agriculture. Several studies have reviewed the correlation between GD and species range and rarity in plants (e.g. Blanco-Pastor et al., 2013), with narrowly distributed species displaying low genetic variability (Binks, Millar & Byrne, 2015) and widespread species showing high GD indices (Gitzendanner & Soltis, 2000). Petunia reitzii and P. saxicola seem to follow the expected depauperate GD for rare species, whereas P. bonjardinensis reveals diversity indices comparable to more widely distributed Petunia spp. (e.g. Turchetto et al., 2015b; to compare some diversity indices with other Petunia spp. and Calibrachoa spp., see Supporting Information, Table S8).

Petunia reitzii and P. saxicola were found in just few sites encompassing a small number of individuals. Moreover, their distribution is reduced and more homogeneous compared with P. bonjardinensis, and for both we observed a more severe reduction in extent of occurrence and area of occupancy in the last two decades. Species narrowly distributed and, particularly those that experienced habitat loss, are expected to have reduced levels of GD and rare species usually have significantly lower levels of polymorphism than their widespread congeners (Gibson, Rice & Stucke, 2008).

We observed that P. bonjardinensis had the highest diversity indices among the three species studied here as well as a larger distribution and higher number of individuals. These indices were also higher than those estimated for other Petunia spp., such as P. mantiqueirensis T.Ando & Hashim. (Backes et al., 2019) that also occurs in a highland grassland formation but c. 1000 km to the north. As P. bonjardinensis, P. mantiqueirensis is known for just a few collection sites and individuals, grows in the same elevation range and is also bee-pollinated.

For some plant species, the reproduction system, seed and pollen dispersal and other life-history traits are considered the main factors controlling GD and population structure (e.g. Thiel-Egenter et al., 2009). Such factors seem did not explain the diversity found in P. bonjardinensis as species that share with P. bonjardinensis such characteristics showed lower GD. Geographically closer to P. mantiqueirensis, Calibrachoa elegans (Miers) Stehmann & Semir is a bee-pollinated species (Mäder & Freitas, 2019) for which the genetic variability is low, even lower than that was observed in P. mantiqueirensis (Backes et al., 2019). These three species display the same floral syndrome and seed and pollen dispersal mechanism. In terms of the demographic parameters, P. bonjardinensis presented a higher population effective size than P. mantiqueirensis and C. elegans. Simultaneously, P. mantiqueirensis and C. elegans have experienced a severe historical bottleneck related to loss of niche suitability during the Mid Holocene (Backes et al., 2019) that was not found for P. bonjardinensis.

Additional factors as historical demography and genetic stochasticity (e.g. Lowe & Allendorf, 2010), habitat specificity and landscape and/or climate changes (Loveless & Hamrick, 1984; Shirk et al., 2014; Ciéslak et al., 2015; Shao et al., 2015) have been implicated in GD in narrowly or widely distributed species. Four Calibrachoa spp. are sympatric to the P. bonjardinensis in SHG, and all displayed similar diversity values (John et al., 2019) compared with P. bonjardinensis, suggesting the influence of similar demographic and landscape features on all.

Petunia bonjardinensis showed GD similar to P. secreta Stehmann & Semir (Turchetto et al., 2016), a lowland species that is also bee-pollinated, apart from the observed heterozygosity that was even higher in P. bonjardinensis. Petunia bonjardinensis and P. secreta have a similar number of known collection sites. However, P. secreta has fewer individuals found during the same flowering season and a significantly lower extent of occurrence and area of occupancy than P. bonjardinensis. Despite its rarity, the high diversity found in P. secreta has been explained based on the seed bank formation and secondary dormancy (Rodrigues et al., 2019), attributes that we were unable to test in P. bonjardinensis.

Although the effective pollinators for SHG Petunia spp. were not yet formally identified, solitary bees (Ceratina, Lanthanomelissa, Hexantheda, Augochlora and Pseudagapostemon; Hymenoptera, Apidae s.l.) were observed visiting P. bonjardinensis and P. reitzii, actively collecting pollen or only drinking nectar (J.R. Stehmann, unpublished data). All these bee genera usually forage on flowers on the ground, such as Petunia spp., and more rarely in the forest canopy, exploring many different species in SHG (Mouga et al., 2016). The behaviour of Pseudagapostemon sp. on P. secreta (Rodrigues et al., 2018) and Pseudagapostemon fluminensis in P. mantiqueirensis (Araújo et al., 2019), collecting pollen and promoting pollination, suggests that the same could occur with P. bonjardinensis and P. reitzii for which an undetermined species of Pseudagapostemon was observed visiting the flowers. Such facts, combined with the proximity between the populations of each Petunia sp. from SHG would explain, at least in part, the low microsatellite genetic structure observed. Similarly, the Calibrachoa spp. endemic to SHG also showed low genetic structure (John et al., 2019). One of these Calibrachoa spp. shows morphological traits attributed to bee-pollinated species, whereas the other two display characteristics linked to the bird-pollination syndrome, as suggested for P. reitzii and P. saxicola (Ando et al., 1999). Thus, the pollination syndrome and reproductive system of these species could not totally explain their GD and the absence of population structure.

Petunia spp. show an autochorous seed dispersion system (Stehmann et al., 2009), in which seeds fall freely when the mature fruit opens (van der Pijl, 1982). This dispersion system reduces the distance between related plants, favouring pollen exchange between relatives and increasing homozygosity and genetic structuring due to the biparental inbreeding as observed in other Petunia spp. (Turchetto et al., 2015a; Rodrigues et al., 2019). For the three species studied here, we did not observe increased mean inbreeding coefficients or significant HWE deviation in the species as a whole, which could be due to the self-incompatibility of these species (Robertson et al., 2011). This may also be one reason why these species are rare, as many pollination events could be unsuccessful. Thus, geographical range or reproductive biology do not seem to explain alone the GD and population structure found in the species from SHG. It is most likely that these species were influenced by historical events as observed in other species from different areas (e.g. Meirmans et al., 2011).

The Pleistocene climatic cycles impacted the SHG formation (Behling, 2002) and probably deeply influenced the endemic species of, e.g. Petunia and Calibrachoa, from this area (Barros et al., 2015, 2020). Petunia is thought to have originated in the Pampas region at the sea level and, under the influence of the glacial and interglacial Pleistocene cycles, secondarily migrated to the highlands (Reck-Kortmann et al., 2014). The same seems to have happened with Calibrachoa (Reck-Kortmann et al., 2015; Mäder & Freitas, 2019). During the glacial periods, cold-adapted grassland species have expanded northwards, with species of Petunia and Calibrachoa reaching the Serra da Mantiqueira, Minas Gerais Brazilian State (Backes et al., 2019). With the interglacial warmer and wetter climate, the Araucaria forest, which during cold and dry weather was restricted to more humid areas along shadowy ravines (Behling, Bauermann & Neves, 2001; Ledru et al., 2005), expanded and fragmented the open fields. Thus, Petunia spp. and Calibrachoa spp. from SHG were isolated and could differentiate, occupying the different elevations where they are currently found. This scenario could explain most the GD found here for these three Petunia spp. that descended from a largely distributed ancestor and diverged in allopatry, confined in small and disconnected areas.

Petunia bonjardinensis, P. reitzii and P. saxicola share a common ancestor (Reck-Kortmann et al., 2014) but occupy areas with different elevations and environmental conditions (Barros et al., 2015). The previously proposed evolutionary history for the highland Petunia spp. based on plastid haplotypes suggested that they have diversified in allopatry under the Pleistocene influence (Lorenz-Lemke et al., 2010). In the present work, we obtained an evolutionary scenario that agrees with these ideas and revealed that P. reitzii and P. saxicola are sister species. The divergence times for all Petunia spp. (c. 1.5 Mya, as estimated by Lorenz-Lemke et al., 2010; c. 2.5 Mya according to Särkinen et al., 2013) was widely superimposed to the Pleistocene period.

Although the three species showed a high level of polymorphism sharing (Supporting Information Table S2; Lorenz-Lemke et al., 2010), this is more likely to be due to inheritance from a common ancestor than to interspecific hybridization. This hypothesis is compatible with the low immigration rates and the results from the genetic structure analyses. Besides, since they occupy different elevations and seem to have different effective pollinators (Stehmann et al., 2009), gene exchange should be more difficult. However, all Petunia spp. can cross, at least under controlled conditions (Gerats & Vadenbussche, 2005), and the three species analysed here share some plastid haplotypes with P. altiplana T.Ando & Hashim. (Lorenz-Lemke et al., 2010) and belong to the same clade in the phylogenetic tree for the genus (Reck-Kortmann et al., 2014). To explain the plastid and microsatellite sharing, we cannot rule out hybridization with another widely distributed Petunia spp. from SHG, P. altiplana, which occurs close to populations of these three congeneric species.

Despite recent habitat loss, rarity is probably historical for these three species, given that we found no genetic evidence for recent bottlenecks. Although P. bonjardinensis presented an effective population size and GD twice as high as the other two, it is compatible with the fact that species has larger geographical distribution and census size compared with the other two analysed here. The N_e and GD absolute values were low, as found also based on plastid sequences (Lorenz-Lemke et al., 2010), suggesting that this species is naturally rare and has persisted in SHG for a long time. The small effective population size and GD of P. reitzii is compatible with the few collection sites and individuals of this species. Although with an effective population size similar to that of P. reitzii, P. saxicola grows in only one site, where paper and lumber industry and intense livestock farming are the main economic activities (Plá et al., 2020). Its present extreme rarity may have been caused by loss of its habitat (see the proportion of loss in extent of occurrence and area of occupancy for this species, Supporting Information, Table S8).

Understanding why species are rare has vital importance to conservation as they can act as biological indicators of overall biodiversity (Wiegand et al., 2020). Grasslands are globally threatened by climate changes and unsustainable land use, which often cause transitions among alternative stable states, and even catastrophic transition to desertification, and it is particularly concerning to highland grasslands (Song et al., 2020). In the last two decades, SHG has suffered severe losses in its wild green coverage (Global Forest Watch, www.globalforestwatch.org, accessed in December 2020; Hansen et al., 2016), mainly due to agriculture and urbanization. Such landscape transformation impacted the three Petunia spp. studied here, reducing the area of occurrence and population sizes and, consequently, their GD, increasing the risk of extinction for these already rare species.

We would suggest that these Petunia spp. should be preserved in situ, but they do not occur in officially preserved lands. An ex situ intervention could be carried out, creating a seed or germplasm bank to re-introduce the plants to their natural environment. It is also essential to begin detailed studies evaluating GD of these species over several generations to create an efficient conservation plan. Demographic and genetic studies are necessary to monitor the current and future decline of these species as stochastic phenomena can weaken the survival ability of populations and initiate an extinction process. Thus, the future status of populations should be carefully followed. Conservation measures should save a maximum number of populations, with priority given to those that are genetically most diverse in these three species.

CONCLUSIONS

Petunia bonjardinensis, P. reitzii and P. saxicola from the Brazilian southern highland grasslands have restricted geographical distributions, and only a few individuals each can be found per reproductive season. Their GD is comparable to that observed in endemic and rare plant species, although P. bonjardinensis display higher indices than other Petunia spp. with similar rarity and even comparable to other species with large range. Based on IUCN criteria these three species are endangered and active strategies for their preservation should be implemented. The main driver for these species diversity and structure was historical and related to the Pleistocene climate changes, secondary migration and founder effect, but currently the landscape fragmentation and habitat loss are strongly influencing them and putting their survival in risk.

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article at the web-site:

Table S1. Sampling information.

Table S2. Microsatellite loci general information.

Table S3. Proposal and prior distributions for parameters for models tested in Migrate-N.

Table S4. Genetic diversity indices per species and microsatellite loci.

Table S5. Pairwise F_ST and R_ST estimates for three Petunia species.

Table S6. Bottleneck signatures for three Petunia species.

Table S7. Recent gene flow among three Petunia species.

Table S8. Diversity indices for Petunia and Calibrachoa species obtained based on microsatellite polymorphism.

Figure S1. Southern subtropical highland grasslands (SHG) landscape.

Figure S2. Conservation areas in southern subtropical highland grasslands related to three Petunia species.

Figure S3. Bar plot for the best number of genetic components as obtained in StructureK = 5. Vertical bars represent individuals and dotted white bars separate the species.

Figure S4. Neighbor-joining tree obtained for individuals of three Petunia species based on DC genetic distances and shared SSR alleles. Small dots indicate P. bonjardinensis (blue), P. reitzii (red) and P. saxicola (yellow) individuals. (BS ≤ 50 for all branches).

Figure S5. Bayesian analysis—posterior distribution over all SSR loci for three Petunia species obtained for the best evolutionary model using Migrate-N.

ACKNOWLEDGEMENTS

We thank G. Mäder, J. N. Fregonezi and J. R. Stehmann for assistance in the field collection and J. R. Stehmann for species identification.

FINANCIAL SUPPORT

This study was supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) and Programa de Pós Graduação em Genética e Biologia Molecular da Universidade Federal do Rio Grande do Sul (PPGBM-UFRGS).

CONFLICT OF INTERESTS

The authors declare that they have no conflict of interest.

REFERENCES