Morphological analysis
Of the three leaf characters (leaf length, leaf width and ratio of leaf length to width), putative hybrid individuals consistently had two morphological characters, leaf length and leaf width, intermediate in value between B. alternifolia and B. crispa (Table 1). Ratios of leaf length to leaf width were significantly larger in B. alternifolia than the other two taxa (Table 1). Of the seven floral characters, corolla tube width and anther height of B. × wardii were intermediate between the values of the two assumed parental species, whereas herkogamy did not differ significantly between the three taxa (Table 1). The remaining four floral characters, corolla tube length (TL), corolla lobe length (CLL), corolla lobe width (CLW) and style length (SL) all showed a similar pattern, with the characters in B. alternifolia having significantly lower values than those measured in B. crispa or B. × wardii (Table 1).
The two putative parental species are morphologically clearly distinct. In the PCA of 10 morphological characteristics, the first and second principal components explained 52.17% and 12.98% of the total variation, respectively. The two-dimensional scatter diagram based on PC1 and PC2 showed clearly the separation of B. alternifolia and B. crispa. Individuals of B. × wardii fell between the two parent species, with a slight overlap with B. alternifolia and a large overlap with B. crispa. Apart from the character HE, there is little difference in the correlation coefficients between the other nine traits (0.29-0.38) (Figure 2 a).
Petal color reflectance
The reflectance spectrum of the petals showed some variation between the different species. In B. crispa and B. × wardii there was a clearly marked peak in the reflectance spectrum at 485 nm, with extremely low variation. However, there was no obvious peak in the reflectance spectrum in B. alternifolia (Figure 2 b).
Sequence analyses of the four nuclear genes
NrETS: The total length of the nrETS region alignment was 380 bp in all individuals, including 20 and 17 nucleotide substitutions in BH and TJ population, respectively (for variation sites, see Additional file 1: Table S1). These variable sites generated 28 haplotypes, among them five, six and thirteen haplotypes belong to B. alternifolia, B. crispa and B. × wardii in BH population, whereas six, three and seven haplotypes in TJ populations, respectively. Three haplotypes were shared by the two populations. Haplotype network analysis identified two major clusters by six nucleotide substitutions. One cluster comprised five haplotypes of B. alternifolia, one haplotype of B. crispa and seven haplotypes of B. × wardii from BH population, and six haplotypes of B. alternifolia, one haplotype of B. crispa and three haplotypes of B. × wardii from TJ population. The other cluster comprised five haplotypes of B. crispa and six haplotypes of B. × wardii from BH population and two haplotypes of B. crispa and four haplotypes of B. × wardii from TJ population (Figure 3 a). Haplotypes of B. crispa from allopatric population were nested in this cluster and forming a small group with three step mutations from the main haplotype (Additional file 9: Figure S2).
Two haplotypes from B. crispa nested with the B. alternifolia cluster, and was found in BHCR2, BHCR7 and TJCR13, which show different haplotypes from both clusters. Of the putative hybrid individuals, all but two individuals (BHWI9 and TJWI1) had two divergent haplotypes nested with each of the two divergent clusters. The remaining two individuals were homozygous for a B. crispa haplotype at this locus (Figure 3 a).
GapC2: The total length of the gapC2 region alignment was 606 bp for all individuals, including 21 nucleotide substitutions, one 2-bp insertion/deletion and one 1-bp insertion/deletion in BH population and 22 nucleotide substitutions and one 1-bp insertion/deletion in TJ population(for variation sites, see Additional file 1: Table S1). A total 21 haplotypes were observed from these loci with three haplotypes were shared between these two populations, among which one, nine and nine haplotypes belong to B. alternifolia, B. crispa and B. × wardii in BH population, whereas three, six and ten haplotypes in TJ populations, respectively. The haplotype network fall into two major clusters separated by 11 nucleotide substitutions. One cluster contained the only one haplotype of B. alternifolia, one haplotype of B. crispa and two haplotypes of B. × wardii from BH population, and two haplotypes of B. alternifolia and three haplotypes of B. × wardii from TJ population. The other cluster contained eight haplotypes of B. crispa and seven haplotypes of B. × wardii from BH population, and six haplotypes of B. crispa and seven haplotypes of B. × wardii from TJ population. Haplotypes of B. crispa from allopatric population were nested in this cluster and shared one haplotype with B. crispa and B. × wardii from BH and TJ (Additional file 9: Figure S2).
Two individuals of B. crispa, BHCR2 and BHCR7, had a haplotype found in the B. alternifolia cluster but all their other haplotypes were found in the B. crispa cluster. In B. alternifolia, there was an exception of two individuals, TJAL6 and TJAL12, which had a haplotype found in the B. crispa cluster, but all the other haplotypes were found in B. alternifolia cluster. All B. × wardii individuals but one (BHWI9) showed two divergent haplotypes originating from both clusters, and the individual BHWI9 was homozygous for a haplotype from B. crispa cluster (Figure 3 b).
PPR24: The total length of the PPR24 region alignment was 647 bp for all individuals, including 43 and 46 nucleotide substitutions in BH and TJ population, respectively (for variation sites, see Additional file 1: Table S1). This gene showed a high level of genetic variation, with 40 haplotypes formed and two haplotypes were shared between the two populations. B. alternifolia, B. crispa and B. × wardii in BH populations had four, ten and fourteen haplotypes, whereas there were one, eleven and fifteen haplotypes in TJ population, respectively. In the haplotypes network analysis, this region was divided into two major clusters separated by 21 nucleotide substitutions. One cluster contained four haplotypes of B. alternifolia, one haplotype of B. crispa and seven haplotypes of B. × wardii from BH population, and the only haplotype of B. alternifolia and eight haplotypes of B. × wardii from TJ population. The other cluster contained nine haplotypes of B. crispa and seven haplotypes of B. × wardii from BH population, and all ten haplotypes of B. crispa and seven haplotypes of B. × wardii from TJ population. One B. crispa individual (BHCR7) and all B. × wardii individuals showed two divergent haplotypes originating from both B. alternifolia and B. crispa clusters (Figure 3 c).
PPR123: After sequence alignment, the PPR123 region was 735 bp in length, and included 23 and 45 nucleotide substitutions in BH and TJ population, respectively (for variation sites, see Additional file 1: Table S1). These variable sites generated 26 haplotypes, among which two, four and eight belong to B. alternifolia, B. crispa and B. × wardii in BH population, whereas two, ten and ten haplotypes in TJ population, respectively. Haplotype network analysis could be divided into three major clades with seven, eight or nine mutational steps between them. One clade comprised two haplotypes of B. alternifolia, one haplotype of B. crispa and five haplotypes of B. × wardii from BH population, and two haplotypes of B. alternifolia and one haplotype of B. × wardii from TJ population. The second clade comprised three haplotypes of B. crispa and three haplotypes of B. × wardii from BH population, and five haplotypes of B. crispa and seven haplotypes of B. × wardii from TJ population. The third clade comprised five haplotypes of B. crispa and two haplotypes of B. × wardii from TJ. (Figure 3 d).
All haplotypes from B. alternifolia individuals fell into the first cluster. For the B. crispa individuals, only one from BH population nested with the first clade derived from BHCR7, ten individuals (TJCR3/5/6/7/9/10/11/12/13/15) had two divergent haplotypes, falling into the second and third clusters, four individuals (TJCR1/4/8/16) had haplotypes nested in the second clade and two individuals (TJCR2/14) had haplotypes nested in the third clade. Of the putative hybrid individuals, all but one individual (BHWI9) had two divergent haplotypes, one in the first cluster and the other in either the second or the third cluster. The remaining BHWI9 had two haplotypes found in the second clade (Figure 3 d).
Sequence analyses for the combined chloroplast regions
The combined length of the cpDNA fragment alignment (including rpl16, trnD-trnT, trnS-trnfM) from the three taxa was 2054 and 2065 bp, and contained ten and five nucleotide substitutions in BH and TJ population, respectively (for variation sites, see Additional file 1: Table S1). Eight haplotypes were inferred in total, and each taxon had three haplotypes in BH population. Haplotype network analysis indicated that most putative hybrid individuals in BH (75%) shared two haplotypes with B. crispa, while they did not share any haplotypes with B. alternifolia. The remaining four B. × wardii individuals (25%) had a unique haplotype with one mutational step away from the predominant haplotype of B. crispa. The haplotype network of cpDNA could not be divided into two clades. Even though B. alternifolia and B. crispa did not share any haplotypes, the differences between their haplotypes were only one or two mutational steps. All individuals in TJ population had only two haplotypes, of them one was derived from six B. alternifolia individuals (TJAL2/3/4/5/8/11), and the remaining 42 individuals shared another haplotype consistent with the haplotype only in individual BHAL13 (Figure 3 e).
NewHybrids analysis
Analysis of the four studied nuclear genes using NewHybrids among the three taxa showed that all individuals with B. alternifolia morphology were assigned to a pure parental species with high posterior probabilities (>0.988). Of the individuals morphologically identified as B. crispa, all but three (BHCR2, BHCR7 and TJCR13) could be identified as B. crispa with high posterior probabilities (>0.996). The individuals BHCR2 and TJCR13 were assigned to B. crispa with lower probability of 0.775 and 0.840, and the individual BHCR7 was assigned to F1 with probabilities of 0.785. Of the individuals morphologically identified as B. × wardii, all but two individuals were assigned to the F1 class with high posterior probabilities (>0.974). The remaining individual BHWI2 was assigned to the F1 class with probability of 0.778 whereas BHWI9 was assigned to be B. crispa with a high posterior probability (0.997) (Figure 4 a).
Population structure analysis
The Bayesian clustering-based structure analysis of all individuals based on all variation loci showed that the highest △K value was found with K=2, indicating that the individuals sampled in this study were divided into two genetic clusters. When K=2, the genetic component of almost all individuals morphologically identified as B. alternifolia mostly formed one cluster and almost all individuals identified as B. crispa formed another cluster. The remaining two individuals of B. alternifolia (TJAL6 and TJAL12) showed low genetic admixture from B. crispa, moreover, one of the remaining three individuals of B. crispa (BHCR7) as well as most B. × wardii individuals showed equal proportions from both clusters. The remaining B. × wardii (BHWI9) and B. crispa (BHCR2 and TJCR 13) individuals had large genetic component belonging to B. crispa (BHWI9: 0.840; BHCR2: 0.755; TJCR13: 0.830), indicating a backcrossed generation (Figure 4 b).