Diversifying selection of the anthocyanin biosynthetic downstream gene UFGT accelerates floral diversity of island Scutellaria species
- Bing-Hong Huang†1,
- Yi-Wen Chen†2,
- Chia-Lung Huang1,
- Jian Gao3 and
- Pei-Chun Liao1Email authorView ORCID ID profile
© The Author(s). 2016
Received: 7 June 2016
Accepted: 5 September 2016
Published: 17 September 2016
Adaptive divergence, which usually explains rapid diversification within island species, might involve the positive selection of genes. Anthocyanin biosynthetic pathway (ABP) genes are important for floral diversity, and are related to stress resistance and pollination, which could be responsible for species diversification. Previous studies have shown that upstream genes of ABP are subject to selective constraints and have a slow evolutionary rate, while the constraints on downstream genes are lower.
In this study, we confirmed these earlier observations of heterogeneous evolutionary rate in upstream gene CHS and the downstream gene UFGT, both of which were expressed in Scutellaria sp. inflorescence buds. We found a higher evolutionary rate and positive selection for UFGT. The codons under positive selection corresponded to the diversified regions, and the presence or absence of an α-helix might produce conformation changes in the Rossmann-like fold structure. The significantly high evolutionary rates for UFGT genes in Taiwanese lineages suggested rapid accumulation of amino acid mutations in island species. The results showed positive selection in closely related species and explained the high diversification of floral patterns in these recently diverged species. In contrast, non-synonymous mutation rate decreases in long-term divergent species could reduce mutational load and prevent the accumulation of deleterious mutations.
Together with the positive selection of transcription factors for ABP genes described in a previous study, these results confirmed that positive selection takes place at a translational level and suggested that the high divergence of ABP genes is related to the floral diversity in island Scutellaria species.
KeywordsAnthocyanin biosynthetic pathway CHS Floral diversity Island species Positive selection UFGT
Floral color patterns are recognized ecologically important factors that influence reproductive isolation and the response to environmental variation . Anthocyanins form a group of plant pigments that contribute to plant diversity in colorful flowers, leaves, stems, etc. Their production typically starts with chalcone synthesis by naringenin-chalcone synthase (CHS, EC 18.104.22.168), followed by a serial synthesis of flavanones, 3-OH-flavanones, flavan-3,4-diols, and anthocyanidins by the core-structural enzymes chalcone isomerase (CHI, EC 22.214.171.124), flavonone 3-hydroxylase (F3H, EC 126.96.36.199), dihydrokaempferol 4-reductase (DFR, EC 188.8.131.52), leucocyanidin oxygenase (LDOX, EC 184.108.40.206, synonym: anthocyanidin synthase, ANS), and UDP-glucose:flavonol 3-O-D-glucosyltransferase (UFGT, EC 220.127.116.11) (cf. ). Upstream and downstream genes of the same pathway might have different evolutionary scenarios. For example, the UFGT genes, which are downstream genes in the ABP (anthocyanin biosynthetic pathway) genes, have relatively higher non-synonymous replacement rates than the upstream genes . Different evolutionary mechanisms shape the genetic diversity of the genes that are upstream and downstream of ABP . Selective constraints and codon usage bias are prominent in CHS-D, while relaxed constraints and more insertion/deletion (indel) events have led to higher UFGT genetic diversity in Ipomoea sp. . Several studies have shown that upstream CHS has a low non-synonymous substitution rate (dN) and this suggests that it evolves with stronger selective constraints than downstream genes, which have an elevated dN [3, 4]. The anthocyanin-regulating transcription factors (TFs) have both a high dN and a synonymous substitution rate (dS), which suggests that there has been a relaxation of selective constraints, but not of positive selection due to ω (=dN / dS) < 1 [5–7]. It has also been suggested that floral diversification is mostly caused by downstream gene diversity [8, 9] and/or diversity in the corresponding TFs [10–12], and that regulatory gene variation is more important than structural mutation for floral anthocyanin diversity in rapidly evolving species . However, mutations of upstream genes (e.g., CHS) typically cause large physio-ecological changes, e.g., the elimination of nearly all flavonoids by the suppression of CHS [14, 15]. They have also been shown to affect pollination [15, 16].
Adaptation could proceed through changes in either TF or structural enzymatic genes, or a combination of the two types . Our previous study showed that the TF R2R3 type MYB11, which regulates the ABP genes PHENYLALANINE AMMONIA LYASE (PAL), CINNAMATE 4-HYDROXYLASE (C4H), CHS, CHI, and UFGT , is expressed in the inflorescence buds and is positively selected for in recently diversified Taiwanese Scutellaria species (Lamiaceae) . The signatures for positive selection pressures on anthocyanin-regulating TFs indicated that the diversification of Taiwanese Scutellaria species is affected by transcription-level evolutionary selection, but it is still unknown whether structural enzyme-level selective pressure (e.g., relaxation of selective constraints or positive selection) is involved. Taiwanese Scutellaria species have arisen many times in mainland China, particularly around c. 0.610 Mya, and have had a short species divergence time. This was followed by very recent within-island speciation events that occurred between 0.204 Mya and 0.015 Mya . Taiwan is a “semi-isolated” island that is separated from continental Asia by the Taiwan Strait, where the narrowest distance between China and Taiwan is 130 km. It was repeatedly connected to continental Asia during glaciation periods, but became isolated during interglacial periods. Repeated connection with the continent enriches species diversity . In addition, steep and rugged terrains create diversified environments and increases opportunities for adaptive radiation.
The phylogenetically close Taiwanese Scutellaria species are distinguished by different corolla colors and spots and their differential growth habitats . The varied floral color pigmentations are often determined by ABP genes and their regulators [11, 22]. These rapidly evolving species may show differential expression and cryptic genetic variation in response to different environmental conditions . The presence of cryptic genetic variation suggests that the plants exhibit low phenotypic variance under typical environmental conditions, but diverge in their response (i.e., high phenotypic divergence) to small genetic variations. For example, differentiation copigment composition and ABP regulator under stress may alter the accumulation or colour pattern of anthocyanins isolated from S. baicalensis [18, 24, 25]. Besides, norms of reactions from different genotype of ABP genes are often associated with environmental adaptation, e.g., light damage , autumnal senescence , and pollination syndrome . Moreover, arising new adaptive genotype with different norm of reactions may lead adaptive radiation . These adaptive genotypes could also be geographic structured due to natural selection and leaded to flower color diversity  and may cause variation of evolutionary rate in ABP gene in different geographic context. Therefore, we can expect that the adaptively evolving genes (e.g., ABP genes in Taiwanese Scutellaria species) will exhibit high ω (dN > dS) in rapidly evolving species, but not in other species. In contrast, if ω is similar between species, then they have similar patterns of molecular divergence.
In this study, we investigated whether the rapid divergence in the floral patterns of Scutellaria species is similar to the case study of Ipomoea  that the flora diversity is a result of the relaxation of selective constraints on genes that are downstream of ABP. We also asked a question that if the rapid divergence of Taiwanese Scutellaria species is associated with adaptive evolution, whether the signatures of positive selection could be detected on ABP genes. By assessing the genetic diversity of these island species, we offer a case of positive selection that accelerates the rate of gene evolution, which may explain the morphological diversification in a group of recently divergent species.
To test the neutrality of ABP genes in Taiwanese Scutellaria, all eight of the native Taiwanese Scutellaria species and an additional 23 species were sampled either from the field or from a seed bank such as B & T World Seed (Additional file 1: Table S1). All plants were grown in a greenhouse at the National Taiwan Normal University, Taipei, Taiwan, and the leaves were collected for DNA extraction.
To obtain the CHS and UFGT sequences, a polymerase chain reaction (PCR) analysis was performed in a MultiGene thermal cycler (Labnet International, Inc., Woodbridge, NJ, US) using 20 ng template DNA, 0.5–1 U Taq (Bernardo Scientific Corp., Taipei, Taiwan), 100 μM deoxyribonucleotide triphosphate, and 0.2 μM of each primer. Primers for CHS amplification were adapted from those used by Chiang et al. . The primers for UFGT amplification were designed from the transcriptomic assembly of inflorescences buds used by Huang et al. . The primer sequences for UFGT were ScUGT-F1: 5′-GTTCCAATGATAGCTCATGG-3′ and ScUGT-R1: 5′-GGAACATAGGCACTCAATTC-3′. The PCR program was set to 94 °C for 3 min for enzyme activation, followed by 35 cycles at 94 °C for 40 s, melting temperature for 40 s, and 72 °C for 80 s, with a 5-min final extension at 72 °C. The PCR products were sequenced directly in both directions using an ABI BigDye 3.1 Terminator Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA). All sequences were visually checked via chromatograms using an ABI PRISM®3730XL DNA Sequencer (Perkin-Elmer, Foster City, CA, USA).
Sequence alignments were conducted using the MUSCLE multiple sequence alignment software tool [31, 32]. Phylogenetic analyses were performed using the neighbor joining, maximum likelihood, and Bayesian method with the assistance of MEGA v.6 , PhyML v.3.0 , MrBayes v.3.2 . The best nucleotide substitution models for CHS and UFGT were the Tamura 3-parameter, with an assumption of a certain fraction of evolutionary invariable sites (T92 + G, gamma shape parameter 0.23), and the Kimura 2-parameter with a discrete gamma distribution (K2P + G, gamma shape parameter 0.42), respectively. Besides, T92 + G model was not available in maximum likelihood and Bayesian methods, so the second best model Hasegawa-Kishino-Yano with a discrete gamma distribution (HKY + G, gamma shape parameter 0.23) was adopted. These were determined by the Bayesian Information Criterion (BIC) method. Pairwise deletions for the treatment of gaps and 1000-times bootstrap replication were set for the phylogenetic reconstruction. To track the corolla colour evolution in skullcaps, we applied methods that originally designed for biogeographic history with consideration of ‘vicariant speciation’ to infer the trait evolution. Three methods were used: the dispersal-vicariance analysis (DIVA), dispersal-extinction-cladogensis analysis (DEC), and BAYAREA method incorporating with founder-event parameter j . In these analyses, ecological characters are analogues to the ancestral distributions and ecological events are analogous to speciation events . These analyses were conducted in BioGeoBEARS .
Detection signals for positive selection
To understand the degrees of genetic divergence within island-based evolving species (implying different selective pressure), the equality of means of evolutionary divergence over Taiwanese species (T/T) pairs dN/dS (ω) for CHS and UFGT was compared to the ω between non-Taiwanese and Taiwanese species (nT/T), and the non-Taiwanese species (nT/nT) Codon-based maximum likelihood models of ω were constructed by the codeml program in PAML v.4.2  and the HyPhy package . First, we estimated the averaged ω from the PAML Model M0, which describes a single ω for all sites. Distributions of pairwise ω values for CHS and UFGT, T/T, nT/T, and nT/nT were compared using Student’s t-test. We also drew a dN-vs.-dS plot and examined the rate saturation of non-synonymous and synonymous mutations. An ω-vs.-dS plot was constructed to illustrate the relative time series of selective pressures, where dS was the time index and ω was the indicator of selection pressure . Besides, ω values were also obtained from ratio of nonsynonymous and synonymous substitution rate calculated from relative rate test implemented in HyPhy packages . Sequences from Tinnea rhodesiana (Lamiaceae) were used as outgroup in relative rate test and ω values were compared using paired t-test. Finally, AMOVA analysis implemented in Arlequin  was conducted to confirm the geographic structure of both CHS and UFGT in Taiwanese and non-Taiwaneses species. The sliding-window analysis of the ω values between Taiwanese and non-Taiwanese species was performed using DnaSP v.5 , and the secondary structures were predicted using the Garnier-Osguthorpe-Robson method [44, 45]. The sliding window analyses and secondary structure predictions help improve understanding of the divergent regions and their corresponding secondary structure. To further confirm signature of positive selection in CHS and UFGT in Taiwanese sister species, McDonald and Kreitman (MK) test  were conducted. We further checked whether the polymorphisms of CHS and UFGT of Taiwanese species attribute to ancestral polymorphism, i.e., the balancing selection, polymorphisms of CHS and UFGT within and between Taiwanese and non-Taiwanese species were compared using Hudson, Kreitman and Aguade (HKA) test . Both MK and HKA tests were implemented in DnaSP v.5 .
The branch model (free-ratio) was tested against the M0 null modes (constant rate model) to detect different ω across lineages, and site models M2a (nearly neutral) and M8 (beta and ω) were tested against the M1a (nearly neutral) and M7 (beta) / M8a (beta and ω fixed to 1) respectively to detect variable selective pressures among sites. The likelihood ratio test (LRT) was used to evaluate the best evolutionary model. However, the PAML analysis does not consider recombination, so we also used the genetic algorithms for recombination detection (GARD), found in HyPhy, to check for recombination . The recombination parameter C  was used to evaluate the degrees of recombination using DnaSP v.5 . The fixed effects likelihood (FEL) and random effect likelihood (REL) models  in HyPhy were used to detect the signatures for the positive selection of specific codons, and the mixed-effects model of evolution (MEME) and branch-site REL (BSR) analysis were used to find the probable episodic selection events on subsets of lineages . Alterations to site-specific amino acid properties through the evolutionary process were diagnosed using the methods followed by Conant et al.  and Atchley et al. , and the property-informed models of evolution (PRIME) model. Each individual null model (no property change) was compared to the full model by LRT using the Chi-square distribution to assess significance.
The CHS and UFGT partial sequences from Scutellaria obtained in this study had alignment lengths of 756 bps (252 codons) and 741 bps (247 codons), respectively. All sequences belonged to an exon. These two genes were found to be analogs of TRANSPARENT TESTA 4 (TT4, E-value = 3 × 10−137) and UDP-glucosyl transferase 73 (UGT73, E-value = 2 × 10−11), respectively, after a tBLASTx search using the Arabidopsis thaliana genome. The number of CHS/UFGT variable sites, parsimonious informative sites, and singletons were 243/227, 159/192, and 84/35, respectively.
The averaged ω values for CHS and UFGT were 0.0269 and 0.5656, respectively. We compared the ω values for CHS and UFGT in Taiwanese species to non-Taiwanese species to ascertain whether the differential selective force act on Taiwanese (island) Scutellaria or just reflected a normal trend in anthocyanin evolution. The relative rate of ω values were calculated using outgroup Tinnea. Relative ω values were significantly different between CHS and UFGT comparisons (bootstrap hypothesis test for equality, P < 2.2e-16). Furthermore, the relative ω values for both genes are significantly different between Taiwanese and non-Taiwanese Scutellaria (two-tailed paired-t test, P = 0.013 for CHS and < 2.2e-16 for UFGT, Additional file 1: Table S3). The relative ω values of CHS were significantly higher in non-Taiwanese Scutellaria (one-tailed paired-t test, P = 0.006) while the relative ω values of UFGT were significantly higher in Taiwanese Scutellaria (one-tailed paired-t test, P < 2.2e-16). This suggested that there was a higher proportion of diversifying selection or relaxation of selective constraints in UFGT of island species, while there was a pervasive purifying selection in CHS of island species.
MK tests also reveal no fixed nonsynonymous substitution in CHS but slightly higher nonsynonymous fixed substitution in UFGT of Taiwanese sister species pairs, although variation is too limited to conduct statistical test (Additional file 1: Table S4). This suggested that UFGT of Taiwanese Scutellaria species suffered from higher divergent selective pressure from other species, but such divergent selection is none or less in CHS of Taiwanese Scutellaria species. Nonsignificant results of HKA test indicated no deviation from neutrality, suggesting no balancing or positive selection acting on both CHS and UFGT in either Taiwanese or non-Taiwanese Scutellaria species (Additional file 1: Table S5). However, such results could be biased due to homogeneous evolving in these two tested genes rather than rejection of positive selection .
Significant geographic structure can be detected in AMOVA in both CHS and UFGT (Additional file 1: Table S6, P < 0.0001 and 0.00684, respectfully), meaning that the Taiwan island plays an isolated environment for the independent evolution of both CHS and UFGT. In other words, the genetic diversification of CHS and UFGT in Taiwanese Scutellaria species is not or less affected by the gene flow from species of other places. Due to elimination of possibilities of ancestral polymorphisms and influence of immigrated genes, we hypothesized that the sources of high polymorphisms of ABP genes of Taiwanese species is driven by natural selection.
Summary of codon-based model analyses of PAML for CHS and UFGT of Scutellaria species
PS site (Pr ω > 1)
ω = 0.02694
ω0 = 0.01830 (p0 = 0.97111)
ω1 = 1.000 (p1 = 0.02889)
ω0 = 0.01830 (p0 = 0.97111)
ω1 = 1.000 (p1 = 0.02889)
ω2 = 35.71235 (p2 = 0)
α = 0.11018, β = 2.88094
α = 0.11218, β = 3.09870
p0 = 0.99761
p1 = 0.00239; ω1 = 1.000
α = 0.11218, β = 3.09870
p0 = 0.99761
p1 = 0.00239; ω1 = 1.000
ω = 0.56561
ω0 = 0.05827 (p0 = 0.62181)
ω1 = 1.000 (p1 = 0.37819)
ω0 = 0.09111 (p0 = 0.63671)
ω1 = 1.000 (p1 = 0.24187)
ω2 = 2.65082 (p2 = 0.12142)
47 N (0.965)
162 M (0.963)
α = 0.02259, β = 0.03044
α = 6.26619, β = 99
p0 = 0.62235
p1 = 0.37765; ω1 = 1.000
α = 0.28974, β = 0.72287
p0 = 0.84684
p1 = 0.15316; ω1 = 2.46858
40 L (0.958)
47 N (0.993)
131 L (0.964)
162 M (0.994)
171 K (0.957)
The inference of positive selection by PAML is usually an argument for not considering the effect of recombination. Therefore, we estimated the recombination rate by parameter C  and used GARD to find the breakpoints in recombination events. The KH topological test was used to check the topological incongruence of the trees on both sides of the breakpoint fragments. In CHS, C is 6.7 per gene and 0.0089 per site, and one breakpoint at codon 381 (ΔAICc = 84.3502) was found to cause different topologies between segments (P < 0.05 in the KH test) by GARD. In UFGT, C was smaller than for CHS (C = 1.8 per gene and 0.0024 per site), and no evidence of recombination was found after using the KH test in the GARD analysis. These results indicated that the CHS gene experienced more recombination events than UFGT.
To reduce recombination effects on the inference of positive selection, the FEL and REL models were used to re-estimate the ω of the codons. In CHS, 136 and 113 negatively selected codons (ω < 1) were detected by FEL at the 0.05 significance level and by REL at Bayes Factor = 50, respectively, but neither model detected sites that had been subjected to pervasive positive selection (ω > 1). In UFGT, three sites (codons 131, 132, and 163) and 12 sites (codons 27, 40, 47, 111, 112, 114, 131, 132, 147, 162, 163, and 171) were shown to have evolved under pervasive positive selection by the FEL and REL models, respectively. Among these 12 positively selected sites detected under REL, only three sites had posterior probabilities > 0.99 (codons 47, 162, and 163). Thirteen negatively selected codons were detected in FEL, but no codon was negatively selected by REL.
In addition, we used the MEME model to test whether specific codons were positively selected on specific lineages, which is defined as episodic positive selection. One in CHS and six in UFGT codon sites were detected under episodic positive selection at the 0.05 significance level (in CHS: codon 39, P = 0.005; in UFGT codon 8, P = 0.035; codon 38, P = 0.011; codon 132, P = 0.048; codon 163, P = 0.015; codon 218, P = 0.0003; and codon 236, and P = 0.039), respectively (Additional file 1: Figure S3, S4). With the exception of codon 8, the other five codons in UFGT that showed episodic positive selection were all found in S. zhongdianensis (Additional file 1: Figure S4). The episodic positive selection pressure on UFGT was also detected in the lineage of S. zhongdianensis under the BSR model, which suggested that 1.68 % of codon sites were positively selected (ω = 155.10, P = 0.009 using the Holm-Bonferroni method). No branches of CHS were subject to episodic positive selection (P < 0.05) in the BSR analysis, suggesting that there is no positive selection acting on CHS of Taiwanese lineages. The PAML and HyPhy analyses suggested that the evolution of CHS is probably under selective constraint and has slow amino acid replacement rates, while multiple episodes of positive Darwinian selection seem to drive UFGT diversity.
Results of Property Informed Models of Evolution (PRIME) showed site-specific amino acid properties in adaptive modification (positive values) and conservation (negative values)
Conant et al.’s method 
Atchley et al.’s method 
Secondary structure factor
Biologists have debated on adaptive evolutionary morphological change are more likely to occur in structural enzyme or regulator [17, 56]. Signatures of positive selection on R2R3-MYB genes that function on regulating ABP genes have been found in Scutellaria . In this study, CHS and UFGT are genes that are upstream and downstream of the ABP, respectively. They are structural enzymes and responsible for the biosynthesis of anthocyanin. Our results supported the previous suggestions of heterogeneous evolutionary rates in these downstream and upstream genes [2–4, 57]. However, most previous studies suggested that pervasive purifying selection pressures on these genes restricted anthocyanin diversity, except for certain episodic positive selections or relaxed constraints on CHS or UFGT genes . In our analyses, genetic signatures for positive selection were detected in UFGT by the site model and free-ratio model produced by PAML, the FEL and REL models produced by HyPhy, the pairwise dN/dS analysis, and the sliding window analysis. Five UFGT codons, 47 N, 131 L, 132S, 162 M, and 163R, were detected in two or all of site-specific models and in the sliding-window analysis of ω (Fig. 5b). In particular, 131 L, 132S, 163R were also detected in FEL constructed by HyPhy, which was suggested to be less false-positively than the other two models . These codons mostly corresponded to the absence or presence of an α-helix in the secondary structure (Fig. 5d). The regions display a structure with β-sheet interspersed among α-helixes that are known as Rossmann-like folds, and form a cleft for substrate binding . The absence or presence of an α-helix may alter the fold structure and may further affect substrate specificity . However, PRIME inferred that while these positively selected UFGT codons did not have any estimated property changes, the 31st codon had an altered charge/isoelectric point (Table 2). This suggested that it is selective pressure drives protein conformation divergence in UFGT, rather than amino acid properties.
The inference that positive selection drives divergence in UFGT conformation, rather than functional divergence, in Scutellaria species is plausible because the anthocyanins are a group of metabolic products involved in stress tolerance, UV resistance, and pollination [14–16, 59]. Such ecologically important products are usually functionally conserved (translational robustness), but show changes in functional (translational) efficiency [60, 61]. A conformation change in a gene sometimes contributes to the adaptive function of meiotic recombination [62, 63]. However, the small recombination rate of UFGT (C = 0.0024 per site and no evidence of recombination by the KH test) excludes the probability of adaptive allele replacement through gene conversion. These small or no recombination and gene conversion events and results of HKA test also reject the hypothesis of balancing selection, which usually retains ancestral polymorphism by recombination and is also present in long-lasting trans-species polymorphism [64–66]. Therefore, we suggest that the multiple amino acid replacements with higher ω values among UFGT lineages (Fig. 5d, Additional file 1: Figure S5) are the relicts of selection for advantageous mutations accumulated throughout evolutionary history.
The average pairwise ω values for at least UFGT in the Taiwanese species were larger than those for non-Taiwanese species, and their dN and dS values are significantly smaller than those of the non-Taiwanese species (Additional file 1: Figure S2 and Table S3), which suggests recent divergent selection. The recent divergence is consistent with the inference of recent positive selection of skullcaps according to the phenomenon of corolla colour transitions on terminal or Taiwanese branches (Additional file 1: Figure S1). Our previous study has shown recently paleoclimate-related species divergences of Taiwanese skullcaps within 0.2 Mya . Here, we provided the other evidence regarding the recent natural selection on genes of ecological traits on these island skullcaps. However, the dN between Taiwanese species is not as small as the dS, especially for UFGT, when compared with long-term divergent species (i.e., non-Taiwanese species) (Additional file 1: Figure S2). This indicates rapid accumulation of non-synonymous mutations in Taiwanese species, and also implies that non-synonymous mutations are not be continually accumulated, which would reduce the high mutational load [67, 68]. In contrast, synonymous mutations are easily accumulated when mutational load was low . If this suggestion is true, we could expect a decrease in the rate of non-synonymous mutation compared to the synonymous mutation rate between species with longer divergent times (e.g., between non-Taiwanese species or between Taiwanese and non-Taiwanese species), which would reduce the damage caused by mutational load. To test this hypothesis, we drew dN-vs.-dS plots. The results showed that the non-Taiwanese lineages (nT/T and nT/nT) had shallower slopes than the Taiwanese lineages (T/T) for both CHS and UFGT. Some even reached a plateau phase (i.e., rate saturation of amino acid change) (Fig. 4). In contrast, the recently diverged Taiwanese species had relatively steep slopes for both CHS and UFGT, which suggested that there was a rapid accumulation of non-synonymous mutations at the beginning of species divergence. The results for the recent selection of varied amino acids by recently diverged species supports the hypothesis of late-burst diversification in Taiwanese Scutellaria species  and local diversifying selection in Taiwanese species.
Geographic structuring variation has been illustrated in both genes, mild proportion of variation can be explained by geographic structure (between Taiwanese and non-Taiwanese, 20.73 % and 17.4 % in CHS and UFGT respectfully, Additional file 1: Table S6). The minor but significantly structure of variation in downstream ABP genes has also been reported in Ipomoea . Besides, the higher dN/dS ratio for ABP genes in Taiwanese species (Fig. 3, especially UFGT in Additional file 1: Table S3) implies that local selection is diversifying Multiple-time origins for Taiwanese species  is an alternative explanation for the high divergence in functional genes. The ω distribution (Fig. 3) was compared to ascertain which of these two hypotheses (i.e., local diversifying selection vs. multiple-time originations) was more likely. If the selection hypothesis were true, a skew towards higher ω values in Taiwanese species rather than non-Taiwanese species would be expected. Alternatively, if the multiple-origin hypothesis is true, we would expect the ω distribution to be similar between Taiwanese and non-Taiwanese species (i.e., functional gene divergence was driven by similar selective pressures). There is a slight skew towards higher ω values in CHS and an obvious low frequency of small ω values in the Taiwanese species (Fig. 2), which supported local diversifying selection.
Taken together, the Taiwanese species had higher amino-acid replacement rates in UFGT (Additional file 1: Table S3), which suggested that rapid divergence counteracted the selective sweeps of different selective pressure in ABP genes amongst these phylogenetically close species. Species within islands have higher encounter rates than continental species or species between islands, which increases the competitive pressure on species expansion (cf. ). Physiological divergence accelerates niche divergence and reduces the cost of species competition . This explains the rapid divergence of ABP genes in Taiwanese species compared to other species (Fig. 3). Both the CHS and UFGT genes sequenced in this study were expressed in the inflorescence buds, according to the tissue-specific RNA-seq of S. indica, S. tashiroi, S. playfairii, and S. taiwanensis (unpublished data), which indicated that these two genes are involved in flower-color diversity . The UFGT expression level is approximately 0.8 ~ 1.3 % that of CHS (unpublished data). It has been suggested that this highly expressed gene evolved slowly [72–74] and that there is a positive correlation between dN and expression specificity [72, 74]. Hence we speculate that the positive selection pressure on UFGT and the rapid divergence in Taiwanese Scutellaria species are probably related to petal color divergence, which also affects pollination [15, 16, 75]. In addition to UFGT, positive selection signals were also detected in the R2R3 type MYB11 and MYB16 genes in Taiwanese Scutellaria species , which help regulate the development of inflorescence axillary meristems [76, 77], and are the microRNA involved in filament development and pollen maturation , respectively. This study, together with previous studies, indicates that genes involving in floral or inflorescence types and colors may be key adaptive characters and may be associated with the rapid diversification of plant species.
The evolution of ABP genes is an important issue because of its ecological relevance to stress resistance and pollination [14–16, 59]. Several studies have indicated genetic conservation in genes that are upstream of ABP genes, whereas there has been a relaxation of selective constraints in downstream genes [2–4]. The results of this study further suggested that UFGT was positively selected in Scutellaria, especially in those phylogenetically closed island species in Taiwan. Such gene divergence was due to changes in protein conformation rather than variations in DNA-binding domains, which prevents mutational load when there is a rapid accumulation of non-synonymous mutations. A previous investigation suggested that there was positive selection of the R2R3-MYB transcription factors, which regulate the transcription of ABP genes. This study supplements the translational information on positive selection signatures for ABP genes. Although less direct evidence supports the hypothesis that there are positive selection pressures on the Taiwanese species, significantly higher ω values imply that rapid divergence of the ABP genes reduced the competitive pressures caused by positive selection within island species. Furthermore, our study extended previous reports with conclusion of relaxation of selective constraints in downstream ABP genes [2, 4], and illustrated that geographic heterogeneity may drive different evolutionary scenarios between different order of the same pathway genes. In this context, our study strongly suggested that floral diversity is an integrative mechanism combining both intrinsic (genetic) and extrinsic (ecological) factors.
The authors thank members of the Liao’s laboratory for sampling assistance.
This research was financially supported by the National Science Council in Taiwan (NSC 102-2621-B-003-005-MY3) and was also subsidized by the National Taiwan Normal University (NTNU), Taiwan.
Availability of data and materials
All relevant data are within the paper and its Supporting Information files.
Conceived and designed the experiments: PCL. Performed the experiments: YWC BHH. Contributed reagents/materials/analysis tools: YWC BHH CLH JG. Analyzed the data: YWC BHH PCL. Wrote the paper: PCL. All authors participated in the discussion, read and approved the final manuscript.
The authors declare that they have no competing interests.
Ethics approval and consent to participate
The authors declare that all experiments described herein comply with the law of government in which they were carried out. All plant materials used in this study are not protected species and are legally obtained and/or collected.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Waser NM. The adaptive nature of floral traits: ideas and evidence. In: Real L, editor. Pollination Biology. Orlando: Academic Press, Inc; 1983. p. 242–86.Google Scholar
- Rausher MD, Miller RE, Tiffin P. Patterns of evolutionary rate variation among genes of the anthocyanin biosynthetic pathway. Mol Biol Evol. 1999;16(2):266–74.View ArticlePubMedGoogle Scholar
- Lu YQ, Rausher MD. Evolutionary rate variation in anthocyanin pathway genes. Mol Biol Evol. 2003;20(11):1844–53.View ArticlePubMedGoogle Scholar
- Rausher MD, Lu YQ, Meyer K. Variation in constraint versus positive selection as an explanation for evolutionary rate variation among anthocyanin genes. J Mol Evol. 2008;67(2):137–44.View ArticlePubMedGoogle Scholar
- Streisfeld MA, Liu D, Rausher MD. Predictable patterns of constraint among anthocyanin-regulating transcription factors in Ipomoea. New Phytol. 2011;191(1):264–74.View ArticlePubMedGoogle Scholar
- Streisfeld MA, Rausher MD. Relaxed constraint and evolutionary rate variation between basic helix-loop-helix floral anthocyanin regulators in Ipomoea. Mol Biol Evol. 2007;24(12):2816–26.View ArticlePubMedGoogle Scholar
- Chang SM, Lu YQ, Rausher MD. Neutral evolution of the nonbinding region of the anthocyanin regulatory gene Ipmyb1 in Ipomoea. Genetics. 2005;170(4):1967–78.View ArticlePubMedPubMed CentralGoogle Scholar
- Kriangphan N, Vuttipongchaikij S, Kittiwongwattana C, Suttangkakul A, Pinmanee P, Sakulsathaporn A, Suwimon R, Suputtitada S, Chanvivattana Y, Apisitwanich S. Effects of sequence and expression of eight anthocyanin biosynthesis genes on floral coloration in four Dendrobium hybrids. Horticult J. 2015;84(1):83–92.View ArticleGoogle Scholar
- Bradley JM, Davies KM, Deroles SC, Bloor SJ, Lewis DH. The maize Lc regulatory gene up-regulates the flavonoid biosynthetic pathway of Petunia. Plant J. 1998;13(3):381–92.View ArticleGoogle Scholar
- Huang B-H, Pang E, Chen Y-W, Cao H, Ruan Y, Liao P-C. Positive selection and functional divergence of R2R3-MYB paralogous genes expressed in inflorescence buds of Scutellaria species (Labiatae). Int J Mol Sci. 2015;16(3):5900–21.View ArticlePubMedPubMed CentralGoogle Scholar
- Yuan YW, Sagawa JM, Frost L, Vela JP, Bradshaw HD. Transcriptional control of floral anthocyanin pigmentation in monkeyflowers (Mimulus). New Phytol. 2014;204(4):1013–27.View ArticlePubMedPubMed CentralGoogle Scholar
- Nakatsuka T, Yamada E, Saito M, Fujita K, Nishihara M. Heterologous expression of gentian MYB1R transcription factors suppresses anthocyanin pigmentation in tobacco flowers. Plant Cell Rep. 2013;32(12):1925–37.View ArticlePubMedGoogle Scholar
- Whittall JB, Voelckel C, Kliebenstein DJ, Hodges SA. Convergence, constraint and the role of gene expression during adaptive radiation: floral anthocyanins in Aquilegia. Mol Ecol. 2006;15(14):4645–57.View ArticlePubMedGoogle Scholar
- Clark ST, Verwoerd WS. A systems approach to identifying correlated gene targets for the loss of colour pigmentation in plants. Bmc Bioinformatics. 2011;12:343.View ArticlePubMedPubMed CentralGoogle Scholar
- Ylstra B, Busscher J, Franken J, Hollman PCH, Mol JNM, Vantunen AJ. Flavonols and fertilization in Petunia hybrida: localization and mode of action during pollen tube growth. Plant J. 1994;6(2):201–12.View ArticleGoogle Scholar
- Mol J, Grotewold E, Koes R. How genes paint flowers and seeds. Trends Plant Sci. 1998;3(6):212–7.View ArticleGoogle Scholar
- Hoekstra HE, Coyne JA. The locus of evolution: evo devo and the genetics of adaptation. Evolution. 2007;61(5):995–1016.View ArticlePubMedGoogle Scholar
- Yuan Y, Wu C, Liu YJ, Yang J, Huang LQ. The Scutellaria baicalensis R2R3-MYB transcription factors modulates flavonoid biosynthesis by regulating GA metabolism in transgenic tobacco plants. PLoS ONE. 2013;8(10):e77275.View ArticlePubMedPubMed CentralGoogle Scholar
- Chiang YC, Huang BH, Liao PC. Diversification, biogeographic pattern, and demographic history of Taiwanese Scutellaria species inferred from nuclear and chloroplast DNA. PLoS ONE. 2012;7(11), e50844.View ArticlePubMedPubMed CentralGoogle Scholar
- Chiang TY, Schaal BA. Phylogeography of plants in Taiwan and the Ryukyu archipelago. Taxon. 2006;55(1):31–41.View ArticleGoogle Scholar
- Hsieh T-H, Huang T-C. Notes on the Flora of Taiwan (20)-Scutellaria (Lamiaceae) in Taiwan. Taiwania. 1995;40(1):35–56.Google Scholar
- Yamagishi M, Toda S, Tasaki K. The novel allele of the LhMYB12 gene is involved in splatter-type spot formation on the flower tepals of Asiatic hybrid lilies (Lilium spp.). New Phytol. 2014;201(3):1009–20.View ArticlePubMedGoogle Scholar
- Ghalambor CK, McKay JK, Carroll SP, Reznick DN. Adaptive versus non-adaptive phenotypic plasticity and the potential for contemporary adaptation in new environments. Funct Ecol. 2007;21(3):394–407.View ArticleGoogle Scholar
- Bkowska A, Kucharska AZ, Oszmiański J. The effects of heating, UV irradiation, and storage on stability of the anthocyanin–polyphenol copigment complex. Food Chem. 2003;81(3):349–55.View ArticleGoogle Scholar
- Yuan Y, Qi L, Yang J, Wu C, Liu Y, Huang L. A Scutellaria baicalensis R2R3-MYB gene, SbMYB8, regulates flavonoid biosynthesis and improves drought stress tolerance in transgenic tobacco. Plant Cell Tissue Organ Cult. 2014;120(3):1–12.Google Scholar
- Steyn WJ, Wand SJE, Holcroft DM, Jacobs G. Anthocyanins in vegetative tissues: a proposed unified function in photoprotection. New Phytol. 2002;155(3):349–61.View ArticleGoogle Scholar
- Tallis MJ, Lin Y, Rogers A, Zhang J, Street NR, Miglietta F, Karnosky DF, De Angelis P, Calfapietra C, Taylor G. The transcriptome of Populus in elevated CO2 reveals increased anthocyanin biosynthesis during delayed autumnal senescence. New Phytol. 2010;186(2):415–28.View ArticlePubMedGoogle Scholar
- Aizza LCB, Dornelas MC. A genomic approach to study anthocyanin synthesis and flower pigmentation in passionflowers. Journal of Nucleic Acids. 2011;2011:371517.View ArticlePubMedPubMed CentralGoogle Scholar
- Brewer MS, Carter RA, Croucher PJP, Gillespie RG. Shifting habitats, morphology, and selective pressures: Developmental polyphenism in an adaptive radiation of Hawaiian spiders. Evolution. 2015;69(1):162–78.View ArticlePubMedGoogle Scholar
- Streisfeld MA, Kohn JR. Contrasting patterns of dloral and molecular variation across a cline in Mimulus aurantiacus. Evolution. 2005;59(12):2548–59.View ArticlePubMedGoogle Scholar
- Edgar RC. MUSCLE: a multiple sequence alignment method with reduced time and space complexity. Bmc Bioinformatics. 2004;5:1–19.View ArticleGoogle Scholar
- Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32(5):1792–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: Molecular Evolutionary Genetics Analysis Version 6.0. Mol Biol Evol. 2013;30(12):2725–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010;59(3):307–21.View ArticlePubMedGoogle Scholar
- Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61(3):539–42.View ArticlePubMedPubMed CentralGoogle Scholar
- Matzke NJ. Model selection in historical biogeography reveals that founder-event speciation is a crucial process in island clades. Syst Biol. 2014;63(6):951–70.View ArticlePubMedGoogle Scholar
- Hardy CR. Reconstructing ancestral ecologies: challenges and possible solutions. Divers Distrib. 2006;12(1):7–19.View ArticleGoogle Scholar
- BioGeoBEARS: biogeography with Bayesian (and likelihood) evolutionary analysis in R scripts [http://CRAN.R-project.org/package=BioGeoBEARS/]. Accessed 1 Mar 2016.
- Yang ZH. PAML 4: Phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91.View ArticlePubMedGoogle Scholar
- Pond SLK, Frost SDW, Muse SV. HyPhy: hypothesis testing using phylogenies. Bioinformatics. 2005;21(5):676–9.View ArticlePubMedGoogle Scholar
- Liao P-C, Chung J-D, Chen C-L, Hwang C-J, Sung Y-H, Chang Y-T, Hwang S-Y. Natural variation, functional divergence, and local adaptation of nucleotide binding site sequences in Rhododendron (Ericaceae). Tree Genet Genomes. 2012;8(4):879–93.View ArticleGoogle Scholar
- Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10(3):564–7.View ArticlePubMedGoogle Scholar
- Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25(11):1451–2.View ArticlePubMedGoogle Scholar
- Garnier J, Osguthorpe DJ, Robson B. Analysis of the accuracy and implications of simple methods for predicting the secondary structure of globular proteins. J Mol Biol. 1978;120(1):97–120.View ArticlePubMedGoogle Scholar
- Kloczkowski A, Ting KL, Jernigan RL, Garnier J. Combining the GOR V algorithm with evolutionary information for protein secondary structure prediction from amino acid sequence. Proteins. 2002;49(2):154–66.View ArticlePubMedGoogle Scholar
- Mcdonald JH, Kreitman M. Adaptive protein evolution at the Adh locus in Drosophila. Nature. 1991;351(6328):652–4.View ArticlePubMedGoogle Scholar
- Hudson RR, Kreitman M, Aguadé M. A test of neutral molecular evolution based on nucleotide data. Genetics. 1987;116(1):153–9.PubMedPubMed CentralGoogle Scholar
- Pond SLK, Posada D, Gravenor MB, Woelk CH, Frost SDW. Automated phylogenetic detection of recombination using a genetic algorithm. Mol Biol Evol. 2006;23(10):1891–901.View ArticleGoogle Scholar
- Hudson RR. Estimating the recombination parameter of a finite population model without selection. Genet Res. 1987;50(3):245–50.View ArticlePubMedGoogle Scholar
- Pond SLK, Frost SDW. Not so different after all: A comparison of methods for detecting amino acid sites under selection. Mol Biol Evol. 2005;22(5):1208–22.View ArticleGoogle Scholar
- Murrell B, Wertheim JO, Moola S, Weighill T, Scheffler K, Pond SLK. Detecting individual sites subject to episodic diversifying selection. Plos Genet. 2012;8(7), e1002764.View ArticlePubMedPubMed CentralGoogle Scholar
- Conant GC, Wagner GP, Stadler PF. Modeling amino acid substitution patterns in orthologous and paralogous genes. Mol Phylogenet Evol. 2007;42(2):298–307.View ArticlePubMedGoogle Scholar
- Atchley WR, Zhao JP, Fernandes AD, Druke T. Solving the protein sequence metric problem. Proc Natl Acad Sci U S A. 2005;102(18):6395–400.View ArticlePubMedPubMed CentralGoogle Scholar
- Wayne ML, Simonsen KL. Statistical tests of neutrality in the age of weak selection. Trends Ecol Evol. 1998;13(6):236–40.View ArticlePubMedGoogle Scholar
- Nielsen R, Yang Z. Likelihood models for detecting positively aelected amino acid sites and applications to the HIV-1 envelope gene. Genetics. 1998;148(3):929–36.PubMedPubMed CentralGoogle Scholar
- Carroll SB. Evo-devo and an expanding evolutionary synthesis: a genetic theory of morphological evolution. Cell. 2008;134(1):25–36.View ArticlePubMedGoogle Scholar
- Toleno DM, Durbin ML, Lundy KE, Clegg MT. Extensive evolutionary rate variation in floral color determining genes in the genus Ipomoea. Plant Spec Biol. 2010;25(1):30–42.View ArticleGoogle Scholar
- Wang XQ. Structure, mechanism and engineering of plant natural product glycosyltransferases. FEBS Lett. 2009;583(20):3303–9.View ArticlePubMedGoogle Scholar
- Winkel-Shirley B. Biosynthesis of flavonoids and effects of stress. Curr Opin Plant Biol. 2002;5(3):218–23.View ArticlePubMedGoogle Scholar
- Park C, Chen XS, Yang JR, Zhang JZ. Differential requirements for mRNA folding partially explain why highly expressed proteins evolve slowly. Proc Natl Acad Sci U S A. 2013;110(8):E678–86.View ArticlePubMedPubMed CentralGoogle Scholar
- Drummond DA, Bloom JD, Adami C, Wilke CO, Arnold FH. Why highly expressed proteins evolve slowly. Proc Natl Acad Sci U S A. 2005;102(40):14338–43.View ArticlePubMedPubMed CentralGoogle Scholar
- Proudfoot AEI, Rose K, Wallace CJA. Conformation-directed recombination of enzyme-activated peptide fragments: a simple and efficient means to protein engineering. Its use in the creation of cytochrome C analogues for structure-function studies. J Biol Chem. 1989;264(15):8764–70.PubMedGoogle Scholar
- Bashkirov VI, Surikov NN, Prozorov AA. Effect of substances altering the DNA molecular conformation on plasmid transformation, plasmid recombination and plasmid transduction in Bacillus subtilis. Genetika. 1982;18(11):1788–92.PubMedGoogle Scholar
- Charlesworth D. Balancing selection and its effects on sequences in nearby genome regions. Plos Genet. 2006;2(4):379–84.View ArticleGoogle Scholar
- Klein J, Sato A, Nagl S, O’hUigin C. Molecular trans-species polymorphism. Annu Rev Ecol Syst. 1998;29:1–21.View ArticleGoogle Scholar
- Kreitman M, Akashi H. Molecular evidence for natural selection. Annu Rev Ecol Syst. 1995;26:403–22.View ArticleGoogle Scholar
- Pybus OG, Rambaut A, Belshaw R, Freckleton RP, Drummond AJ, Holmes EC. Phylogenetic evidence for deleterious mutation load in RNA viruses and its contribution to viral evolution. Mol Biol Evol. 2007;24(3):845–52.View ArticlePubMedGoogle Scholar
- Wielgoss S, Barrick JE, Tenaillon O, Wiser MJ, Dittmar WJ, Cruveiller S, Chane-Woon-Ming B, Medigue C, Lenski RE, Schneider D. Mutation rate dynamics in a bacterial population reflect tension between adaptation and genetic load. Proc Natl Acad Sci U S A. 2013;110(1):222–7.View ArticlePubMedGoogle Scholar
- Li W-H, Gojobori T, Nei M. Pseudogenes as a paradigm of neutral evolution. Nature. 1981;292(5820):237–9.View ArticlePubMedGoogle Scholar
- Svenning JC, Gravel D, Holt RD, Schurr FM, Thuiller W, Munkemuller T, Schiffers KH, Dullinger S, Edwards TC, Hickler T, et al. The influence of interspecific interactions on species range expansion rates. Ecography. 2014;37(12):1198–209.View ArticlePubMedPubMed CentralGoogle Scholar
- Westoby M, Falster DS, Moles AT, Vesk PA, Wright IJ. Plant ecological strategies: Some leading dimensions of variation between species. Annu Rev Ecol Syst. 2002;33:125–59.View ArticleGoogle Scholar
- Yang L, Gaut BS. Factors that contribute to variation in evolutionary rate among Arabidopsis genes. Mol Biol Evol. 2011;28(8):2359–69.View ArticlePubMedGoogle Scholar
- Pal C, Papp B, Hurst LD. Highly expressed genes in yeast evolve slowly. Genetics. 2001;158(2):927–31.PubMedPubMed CentralGoogle Scholar
- Williamson RJ, Josephs EB, Platts AE, Hazzouri KM, Haudry A, Blanchette M, Wright SI. Evidence for widespread positive and negative selection in coding and conserved noncoding regions of Capsella grandiflora. Plos Genet. 2014;10(9), e1004622.View ArticlePubMedPubMed CentralGoogle Scholar
- Malcomber ST, Preston JC, Reinheimer R, Kossuth J, Kellogg EA. Developmental gene evolution and the origin of grass inflorescence diversity. Adv Bot Res. 2006;44:425–81.View ArticleGoogle Scholar
- Keller T, Abbott J, Moritz T, Doerner P. Arabidopsis REGULATOR OF AXILLARY MERISTEMS1 controls a leaf axil stem cell niche and modulates vegetative development. The Plant Cell Online. 2006;18(3):598–611.View ArticleGoogle Scholar
- Müller D, Schmitz G, Theres K. Blind homologous R2R3 Myb genes control the pattern of lateral meristem initiation in Arabidopsis. The Plant Cell Online. 2006;18(3):586–97.View ArticleGoogle Scholar
- Millar AA, Gubler F. The Arabidopsis GAMYB-like genes, MYB33 and MYB65, are MicroRNA-regulated genes that redundantly facilitate anther development. Plant Cell. 2005;17(3):705–21.View ArticlePubMedPubMed CentralGoogle Scholar