Plant species, study site and sampling
The two study species P. poissonii and P. anisodora are perennial herbs, growing on alpine meadows at altitudes between 2500 and 3700 m [23]. While P. poissonii is relatively widespread, occurring in the eastern Himalayas (Tibet, Yunnan and Sichuan), P. anisodora is only known from three populations in Yunnan province (China): two populations (Langduand Baishuitai) in Shangri-La county and one population (Xiaoyanjing) in Nujiang county (Additional file 6: Figure S1). We only know of one area where both species grow in sympatry and seem to form hybrids (Baishuitai). This hybrid zone spans an area of about 15 km2 in Baishuitai, Shangri-La, NW Yunnan (27.6°N, 100.04°W, 3100 m). The population is situated in wet alpine meadow, with a stream running through it. In 2017, individual counts for P. anisodora and P. poissonii were 517 and 2349 flowering plants, respectively.
Parental species and hybrids can be easily distinguished by flower color and floral orientation: P. anisodora has dark-red, drooping flowers; P. poissonii has magenta-colored flowers that are oriented perpendicular to the stem; putative hybrids have intermediate color and floral orientation (Fig. 5a, b, c).
During the flowering period in 2014, a total of 168 individuals were sampled: 108 from the hybrid zone (P. anisodora: 36 individuals; P. poissonii: 33 individuals; putative hybrids: 39 individuals); 20 each from the two allopatric P. anisodora populations - Xiaoyanjing (26.58°N, 99.44°W), Langdu (27.95°N, 99.7°W); and 20 P. poissonii from an allopatric reference population near Shangri-La county (28.2°N, 100.04°W). Leaves were collected and preserved in silica gel in individually numbered containers. We stated here that all samples collected and used in this study as well as pollination experiments performed below in the field do not need any permission. We declare plant materials used in the current research complied with government regulations.
Reproductive barriers in the hybrid zone
The flowering times of the two species overlap completely, starting early in July, and ending in mid August. Therefore we predicted that pollinator-mediated isolation and intrinsic incompatibilities could play a role for isolating parental species and hybrids, and might also determine hybridization patterns in the sympatric area.
Pollinator-mediated reproductive isolation
Pollinator observations were carried out in two types of setting: (1) natural setting - from 6 July to 10 July in 2015 pollinators were observed in natural either P. poissonii- or P. anisodora-dominated parts of the population, to assess their pollinator assemblages; and (2) artificially controlled plot setting - three plots were established within the hybrid zone to evaluate reproductive barriers between parental species. The first of these plots (Plot 1) was set up in a P. poissonii-dominated part of the population and comprised both parental species (10 plants each) in an alternating setup (Fig. 5d); Plot 2 had an identical layout to plot 1, but was placed in a P. anisodora-dominated part of the population (Fig. 5f); Plot 3 was placed inside the hybrid zone, and comprised both parents and hybrids (Fig. 5e). Individual plants in the grid setups were placed approximately 50 cm apart, and species / hybrid identity was assigned according to morphological characters. From 6 July to 12 July in 2015, both Plot 1 and Plot 2 were observed for 7 days; each observation day lasting from 09:00 to 17:00. However, if weather conditions were unfavorable (rain) for pollinators, field observations were stopped during those hours. Plot 3 was observed on five sunny days from 7 to 11 July, 2015. We recorded a flower approach by a pollinator as a visit only if the pollinator had foraged or interacted with the flower by touching the reproductive structures of the plant. All types of flower visitors were photo recorded, identified, and preserved in the insect collections of Kunming Institute of Botany, CAS.
Following Natalis and Wesselingh [18], we calculated overall ethological reproductive isolation (RI) in P. poissonii- as well as P. anisodora-dominated parts of the population as: 1 - (No. cross-species foraging bouts / total number of foraging bouts) × (heterospecific transitions / all transitions). Herein one bout means the visitation of an insect to the plot, from approaching the plot up to its departure, irrespective of how many plants within the plot were visited, whereas transition refers to the type (conspecific vs heterospecifc) and number of visitations within one bout.
Hand pollination experiments
Pollination experiments were carried out in 2015 and 2016. The main focus between the years was different: in 2015 cross-compatibility between the species was tested, and in 2016 fertility between hybrids and parental species was investigated. Although it would have been desirable to carry out both of the experiments in either year, this was not possible, as the number of flowering plants of P. anisodora would not allow appropriate sample sizes for both experiments in the same year.
Cross-compatibility of the parental species
In July 2015, cross-pollination experiments were carried out for 16 possible combinations between the flower-morph types and species. Because both species are heterostylous, four intra specific pollination treatments are possible (♀-♂): P-P; P-T; T-P; T-T. These same four treatments also apply to inter-specific cross-pollination, and hence 16 treatments are necessary to cover all possible combinations.
For each treatment a minimum of 30 plants were selected as mothers. Up to three flowers per plant were selected while in bud, and the rest removed mechanically; afterwards the plant was labeled and bagged with a nylon net. When the flowers had just opened, flowers of the respective pollen parent were collected, and recipient flowers were pollinated manually by brushing the anthers against the stigma. After application of pollen the nylon net-bag was replaced and left for a further 3 days. In September 2015, fruits were harvested, and seed numbers per fruit were counted. Several plants were damaged by grazing before they could be harvested, resulting in 15 to 35 flowers per treatment that could be assessed and which are given as sample size ( Additional file 1: Table S1).
Fertility between hybrids and parental species
In July 2016, we focused on barriers to backcrossing between hybrids and parental species. To this end we performed cross-pollinations between hybrids and the parental species. As before, the different plant morph types were treated separately, resulting in 16 treatment combinations (Additional file 1: Table S1). The experiment was carried out as in 2015, but for some treatments grazing damage was higher, resulting in assessable sample sizes between 12 and 33 pollinated flowers. Fruits were again harvested in September, 2016.
Data analysis
To detect differences of visitation percentages between parental species, a χ2 test was employed. For assessing the effect of factors below on fruit- and seed-set for the pollination treatment data we employed ANOVAs as implemented in SPSS 15.0 for Windows (Chicago, IL, USA).
For the 2015 data we included four factors: Mother Species (anisodora, poissonii); Flower Type of mother plant (pin, thrum); Pollen Source (anisodora, poissonii); and Cross Type (intra-morph, inter-morph). Additionally a three-way ANOVA analysis for intra-specific pollination treatments was done including the factors Mother Species, Flower Type and Pollen Source. In order to improve ANOVA assumptions of normality and homogeneity, fruit sets were arcsine square-root transformed whereas seed numbers per fruit and seed numbers per flower were both log (x + 1) transformed.
For the 2016 data 7 treatments did not produce any fruits. We first compared differences in fruit sets and seed numbers between cross-pollinations of P. poissonii x hybrids and P. anisodora x hybrids using a χ2 test (fruit sets for pollination treatments of P. poissonii x hybrids and P. anisodora x hybrids were calculated as the total number of fruits divided by the total number of flowers within each treatment) and a Mann-Whitney nonparametric test (mean seed numbers per treatment as replicates, thus both P. poissonii x hybrids and P. anisodora x hybrids cross-pollinations comprising 8 replicates).
Molecular analyses
DNA extraction and microsatellite genotyping
We extracted genomic DNA from the dried leaf tissue using a modified cetyl trimethyl ammonium (CTAB) protocol. Quantification of DNA was carried out with a SmartSpecTM Plus Spectrophotometer (Bio-Rad). To select candidate diagnostic microsatellites to identify parental species and hybrids, 219 SSR primers were designed based on sequences previously obtained on a MiSeq Benchtop Sequencer (Illumina, Inc., San Diego, CA, USA) for the close relative P. chungensis [24]. Detailed protocols for PCR amplification and condition regarding these primers followed Zhou et al. [24]. PCR products were directly analyzed on a 3730xl Sequence Analyzer (Applied Biosystems, Foster City, CA, USA), using a LIZ GeneScan-500 size standard. Resulting chromatograms were visualized and converted to diploid genotypes using automated allele-calling implemented in GENEMARKER v.4.0 (SoftGenetics LLC, State College, PA, USA). All automated genotyping was re-checked manually. All genotypes for each locus and individual were entered into an Excel file following the format of GenALEx 6.5 [25].
Population genetic analysis
All basic summary statistics for the microsatellite data, and an AMOVA were calculated using GenALEx 6.5. Deviation from Hardy–Weinberg equilibrium (HWE) for each locus, and linkage disequilibrium (LD) for all loci pairs, were assessed with the web version of GENEPOP, v4.0.10 [26]. For the HW test the heterozygote deficiency option was used. For the LD test the following Markov Chain parameters were set: 5000 dememorization runs; 1000 batches each with 5000 iterations.
Parent and hybrid assignment
To estimate parental ancestry proportions for each hybrid, microsatellite data was analyzed using the program STRUCTURE version 2.3.1 [27]. We adopted the admixture model with correlated allele frequencies [19]. No prior knowledge of the species was included in the analyzed data set. To determine the optimal number of groups (K), we ran STRUCTURE with K varying from 1 to 10, with five runs for each K value, and following Evanno et al. (2005) used ΔK to choose the most likely number of clusters. MCMC runs were then performed using a burn-in of 10,000 followed by 10,000 iterations.
Hybrid class assessment based on the microsatellite data was carried out with the program NewHybrids [28], using the following settings: burn-in of 10,000 followed by 100,000 MCMC iterations with the Uniform priorsmode. For the classification of individuals into each of the three groups, a relatively strict posterior probability of 0.9 was used, and only individuals reaching such a probability were taken as being classified; others were considered not properly classified by the program.
Hybrid swarm simulation
As the analyses were based on six loci only, we assessed their discriminatory power in NewHybrids by simulating several classes of hybrids based on parental individuals that had membership to either cluster in STRUCTURE close to 1 (‘pure’ parents). We simulated 5 hybrid categories (two parental species, F1and the first backcross to each parent) from 23 each of pure P. anisodora and P. poissonii (both > 90% probability in both structure and NewHybrids) using HYBRIDLAB version 1.0 [29]. For each category, 30 individuals were generated, hence a total of 150 individuals (30 of each parent, 30 F1s and 30 BC to each parent) were simulated for the hybrid swarm in Baishuitai, and these were subsequently analyzed in NewHybrids.