Progress in phylogenetics, multi-omics and flower coloration studies in Rhododendron

Accurate species delimitations and robust phylogenies are crucial for botanists, ecologists, and horticulturists alike to study and utilize plant resources. The genus Rhododendron (Ericaceae) is particularly appreciated in horticulture for its magnificent flowers^[1−5]. At the same time, with about 1,000 species, it is the largest genus of woody plants in the northern hemisphere. Rhododendron (L.) was first described by the Swedish botanist Carl Linnaeus (Carl von Linné, 1707−1778) in 1753 with only five species. In addition, Linnaeus described the closely related genera Azalea and Rhodora with six and one species, respectively. Subsequently, a larger number of Rhododendron, Azalea and Rhodora species were described, including those from the temperate zones of the northern hemisphere. The first tropical species R. malayanum was described in 1822, then the German botanist Carl Ludwig Blume (1796−1862) established the genus Vireya with four species from Java and Celebes. However, the related genera Azalea, Vireya, and Rhodora were combined into the genus Rhododendron^[5]. The re-defined genus Rhododendron was divided into eight sections by the Baltic-German botanist Karl Maximovich (1827−1891) in 1870, which had a significant impact on later classification systems, albeit with different names or categories^[3].

In the late 19^th and early 20^th centuries, the extensive collections of Rhododendron in southwest China led to the discovery of numerous species and the establishment of many series. These achievements were compiled in a handbook entitled 'The species of Rhododendron'^[6], which included 43 series and about 850 species of temperate and subtropical zones. However, no formal classification system was established. It was not until 1949 that Sleumer proposed a comprehensive taxonomic system for Rhododendron, comprising eight subgenera and 13 sections, and presented a monograph covering all Rhododendron species^[5], especially those from tropical Indochina, Malesia and Australia. Sleumer's system was later revised by Chamberlain et al., the main difference being the treatment of the subgenus Therorhodion^[1], which had been removed from Rhododendron by Sleumer^[5]. Chamberlain's treatment was largely confirmed by further phylogenetic analysis based on morphological characters and molecular data^[2], and the molecular evidence supported the position of subgenus Therorhodion within the genus Rhododendron. Finally, morphologically based classifications of Rhododendron have led to a consensus taxonomic framework that recognized the five major subgenera Azaleastrum, Hymenanthes, Pentanthera, Rhododendron, Tsutsusi, and the three minor ones Candidastrum, Mumeazalea, and Therorhodion^[1].

Recently, the traditional classification systems have been challenged by new molecular findings. Based on molecular and morphological data, Goetsch et al. further revised Chamberlain's system, whereby the subgenera Candidastrum, Mumeazalea, and Tsutsusi were merged into subgenus Azaleastrum and the polyphyletic subgenus Pentanthera was dissolved by transferring the sections Pentanthera and Rhodora to the subgenus Hymenanthes and the sections Sciadorhodion and Viscidula to the subgenus Azaleastrum^[1,2].

More recently, Xia et al. reconstructed a molecular phylogeny of Rhododendron comprising 200 species that represent all subgenera, sections, and nearly all multi-species subsections, using 3,437 orthologous nuclear genes generated from transcriptome data^[4]. They provided the first robust and dated backbone phylogeny based on genome-level data across the genus. Xia et al. recognized five subgenera and 11 sections and resolved the relationships of most taxonomically problematic groups of Rhododendron and its close relatives^[4] (Fig. 1). For example, subgenus Choniastrum was treated as a section of the subgenus Tsutsusi, and the genus Menziesia was merged into section Sciadorhodion under the subgenus Tsutsusi. They recognized the section Schistanche, which comprises 312 species, over 95% of which are endemic to the Malay Archipelago. A comparison with the recently published plastome phylogeny for Rhododendron revealed conflicts between nuclear and plastid phylogenetic topologies, suggesting that reticulation events may have occurred in the deep lineages of the genus^[7].

Figure 1. 

Phylogenetic reconstruction in Rhododendron L. Sect., Section; Subg., Subgenus. (Modified from a previous publication^[4]).

Xia et al. have shown that phylo-transcriptomics is an extremely efficient approach for reconstructing phylogenies, even for plant groups with rapid radiations such as Rhododendron^[4] (Fig. 1). In the updated phylogeny, subgenus Therorhodion, previously considered a distinct genus, is the earliest diverged branch and sister to all other Rhododendron. This group includes two species distributed in the boreal tundra of East Asia and North America. Then there is the subgenus Tsutsusi, a monophyletic clade consisting of the sections Sciadorhodion, Choniastrum, Azaleastrum, Brachycalyx and Tsutsusi and comprising around 130 species, which are mainly distributed in East Asia. The sections Choniastrum and Azaleastrum were previously recognized as separate subgenera based on morphological and/or molecular evidence^[2]. In section Sciadorhodion, R. albiflorum from western North America is a sister species to species from northeast Asia^[4]. Species of section Tsutsusi are important genetic resources for evergreen azaleas, which are popular as ornamental shrubs or pot azaleas. Some traditional Chinese and Japanese varieties have been cultivated for more than 400 years. Among the 50 Tsutsusi species native to Japan, some evergreen species such as R. kaempferi, R. macrosepalum, R. indicum, and R. ripense are highly decorative, and there are hundreds of azalea cultivars selected from natural populations and their hybrids^[8,9]. The various molecular data is enabling the potential to unravel the complex phylogeny of the genus Rhododendron based on high-throughput sequencing. In a recent study, around 800,000 high-quality SNPs were obtained with restriction-site associated DNA sequencing (RAD-seq)^[10]. Based on these large sets of SNPs, Shen et al. decoded the taxonomic relationships within the genus Rhododendron, reconstructed ancestral states of Rhododendron, and showcased the feasibility of applying RAD-seq to unravel intricate evolutionary relationships despite a significant amount of missing data^[10]. With the further development of high-throughput sequencing, it is believed that deep whole-genome DNA sequencing will bring about many new discoveries in the future.

About 200 years ago, Belgian pot azaleas with a relatively narrow genetic base were developed from collector material from China and/or Japan, with several forms of R. simsii serving as potential parents, and R. indicum, R. scabrum and R. mucronatum may also have contributed to the modern R. simsii hybrids^[8,11]. Azalea cultivars in China are categorized into four main types: (1) Eastern azaleas, which have small leaves and mostly originate from Japan; (2) Western azaleas (Belgian azaleas), including plants with large flowers and double and semi-double petals returning to China from Europe; (3) Hairy azaleas, which are characterized by large simple flowers and finely pubescent stems and leaves and are cultivated in China; and (4) Summer azaleas, which originate from Japan with cultivars blooming in early summer. However, the delimitation of these types is not supported by molecular evidence, and different cultivars may have the same name or a single plant may have multiple names^[12].

There are over 500 species in the largest subgenus Rhododendron. These species are mainly distributed in the Himalayan-Hengduan Mountains (HHM) and the Malay Archipelago and can be identified by the scales on the back of their leaves. R. hypoleucum and R. tomentosum from circumboreal regions form a clade that is the first to diverge^[4], and the sister clade comprises two reciprocal lineages: section Schistanthe, which occurs mainly in the Malay Archipelago, and section Rhododendron, which includes the Asian species (mainly in the HHM) and two European species (R. ferrugineum and R. hirsutum). Subgenus Pentanthera was previously considered as a section of subgenus Hymenanthes^[2], but Xia et al. showed that it is related to subgenus Hymenanthes^[4] and includes more than 20 deciduous species disjunctly distributed between North America, East Asia, and Western Europe. The subgenus Hymenanthes, on the other hand, comprises about 300 evergreen species, which are mainly distributed in the Hengduan Mountains.

Analyses of whole-genome resequencing data revealed that the relationships among species within the subgenus Hymenanthes were largely inconsistent with phylogenies inferred from different methods or datasets, especially for species from southwestern China, suggesting extensive hybridization and genomic admixture during the historical radiation of the subgenus, and repeated isolation and hybridization are thought to have promoted rapid diversification of species in mountainous regions^[13]. In addition, most taxonomic treatments based on morphological characters have been rejected in this subgenus, and some subsections, such as the subsections Fortunea and Argyrophylla, have not been supported as monophyly^[13].

Flower color evolution

Flower colors have been well studied by early evolutionary biologists in ecological and traditional genetic contexts, and it has been suggested that the attractive color display is mainly determined by three main groups of pigments, namely anthocyanins, carotenoids and betalains^[14−17]. Rhododendrons are valued for their great ornamental beauty and a wide range of flower colors: red, purplish pink, purple, yellow, white, and various intermediate colors. Until now, the mechanism underlying flower coloration in Rhododendron has rarely been studied, with the main focus on cultivated azaleas. Co-pigmentation between anthocyanins and flavonols has been found to be responsible for petal coloration in species and cultivars of the subgenus Tsutsutsi, with color intensity being influenced by the amount of anthocyanins, and increasing amounts of flavonols leading to a color change from red to purple^[18]. In many sections of Rhododendron, the carotenoids are not involved in the coloration of the petals but appear in the form of spots and blotches^[18].

Flower color is often correlated with other floral traits, leading to the common use of the term 'pollination syndromes' to refer to adaption to pollinator regimes^[19]. Recently, Du et al. studied the flower colors of 30 representative Rhododendron species from seven subgenera of the genus and identified seven anthocyanins (cyanidin 3-O-galactoside, cyanidin 3-O-arabinoside, cyanidin 3,5-di-O-glucoside, cyanidin 3-O-arabinoside-5-O-glucoside, delphinidin 3-O-glucoside, delphinidin 3-O-arabinoside-5-O-glucoside, malvidin 3-O-arabinoside-5-O-glucoside) were responsible for the different colors of Rhododendron flowers, leading to the conclusion that cyanidin-type anthocyanins were responsible for the red coloration and delphinidin- and malvidin-type anthocyanins for the purple coloration^[20]. However, they also identified 21 flavonol derivatives that may not play a major role in flower coloration^[20]. Their results improve our understanding of flower coloration, genetic variation and environmental adaptation in Rhododendron. It seems that Rhododendron has comparatively uniform pigment types but different amounts and components in different classification levels^[20], which is supported by morphological and molecular evidence. Xia et al. found that red, purplish pink, purple and yellow flowers occurred repeatedly in the major clades and subclades of the Rhododendron phylogeny, indicating widespread parallel evolution of flower colors in this genus^[4]. Variations in flower colors among closely related species are often associated with rapid shifts in specific pollinators^[21−23], highlighting the rapid adaptation of plants to pollinators and species diversity in the species-rich genus Rhododendron driven by biotic factors.

Over the past 20 years, the sequencing of plant genomes has increased rapidly, leading to an increase in the quantity and quality of publicly available genomic resources^[24]. Genomic information on Rhododendron, one of the largest woody plant genera, has been acquired and disseminated by several research institutions. Since the publication of the draft genome assembly of R. delavayi in 2017, a total of 14 genomes from 12 species have been sequenced^[13,25−34]. Supplemental Table S1 shows that most genome assemblies (all but one) have reached the chromosome level, with three assemblies approaching a near-complete T2T level. Recently, the R. vialii genome assemblies have reached the haplotype level^[35]. The improvement of assembly quality, especially in terms of completeness, is in line with the development of sequencing technologies. PacBio HiFi sequencing technology generates long-read sequencing datasets with average read lengths of 10−25 kb and an accuracy of over 99.5%^[36,37]. Combined with Hi-C sequencing technology and newly developed assembly techniques, it will become easier to obtain more high-quality Rhododendron genome assemblies. This undoubtedly represents a great opportunity for Rhododendron species with high genome heterozygosity.

Among the five subgenera of Rhododendron, Hymenanthes has the most sequenced species (eight published genomes from seven species), followed by Tsutsusi (four genomes), and Pentanthera with one species (two assemblies of R. molle) (Supplemental Table S1). No genome assemblies were reported for the other two subgenera, Therorhodion and the larger subgenus Rhododendron (Supplemental Table S1). Moreover, only one of the sequenced species is deciduous (R. molle), while the remaining 11 species are evergreen Rhododendrons (Supplemental Table 1). Clearly, more importance should be given to larger and more evenly sampled genome projects. Rhododendron plants exhibit a wide range of flower colors, broadly classified into red, white, purple, and yellow color series^{[13,16,18,20,38,39]}. The sequenced genomes covered four color series and provided the opportunity to investigate the genetic basis of flower color transitions at the genomic level (Supplemental Table S1). To fully utilize these resources, a website was developed to integrate and make these data accessible. Here we have compiled download links for different data types to facilitate access for researchers (Table 1). In addition, two Rhododendron-specific databases^[40] are accessible: RPGD (http://bioinfor.kib.ac.cn/RPGD/) and http://rhododendron.plantgenie.org/.

Genome sizes and divergence

Genome size is commonly referred to as the amount of DNA present in the cell nucleus^[41]. Genome size in land plants is influenced by several factors, such as whole genome duplication, polyploidization, and amplification and elimination of repetitive sequences^[42−44]. All sequenced Rhododendrons are diploid with 13 chromosomes in the reported genome assemblies. The variation in genome size is not significant and ranges from ~500 Mb in R. ripense to ~700 Mb in R. irroratum (Supplemental Table S1). Comparative genomic analyses revealed that the increase in genome size of R. molle (653.46 Mb) compared to R. simsii (528.64 Mb) was due to the recent proliferation of long terminal repeat retrotransposons (LTR-RTs), particularly Gypsy, resulting in an increase of 125 Mb (19%)^[26,31]. A positive correlation was observed between the length of repetitive sequences and genome size, suggesting that repetitive sequences may influence genome size variation in Rhododendron (Supplemental Table S1).

Rhododendron species show differentiation in inter-genomic sequences but are still conserved to a certain extent. With R. simsii genome assembly as a reference, a comparative analysis of the whole genome sequences revealed that 60% of the R. molle genome sequence exhibits inter-genomic collinearity, while 36% shows sequence differentiation^[26,31]. LTR-RTs are considered to be the main driver of interspecific sequence variation. In particular, Gypsy elements could play an important role in the centromeric regions of chromosomes and contribute to the rapid evolution of interspecific centromeric sequences^[26,31]. In the reported Rhododendron genome assemblies, LTR-RTs were identified as the predominant type of repetitive sequences, with Gypsy being the most abundant LTR-RTs element. It is worthwhile to investigate the evolutionary patterns of Gypsy-like LTR-RTs in Rhododendron and their effects on genome size and sequence variation.

Whole genome duplication (WGD) has been another focus in Rhododendron, thought to be another force contributing to changes in genome size and sequence. It is generally agreed that the genus Rhododendron has experienced two common WGD events: the first is the ancient whole genome triplication event (gamma) common to all core eudicots, while the second is the WGD event associated with the origin of the entire Ericales order. The occurrence of the most recent WGD is the subject of debate, with three main perspectives being advocated: (1) the Ericales^[26,45−48]; (2) the core Ericales^[49−51]; and (3) the Ericaceae and Actinidiaceae^[31,51]. The complexity of gene loss and retention patterns after WGD makes it difficult to determine the origin and phylogenetic position of Ericales. Numerous techniques have been used to study and identify polyploidy events. These include non-phylogenetic methods such as karyotyping and synteny^[52,53] as well as phylogenetic methods such as Ks-based analysis, least common ancestor (LCA) reconciliation, count-based, and gene networks approaches^[54−57]. Despite the application of these detection methods in various studies^{[46−50,58,59]}, the exact timing of the ancient WGD event in the order Ericales remains unclear.

A fundamental problem is the difficulty of separating WGD from Ericales differentiation. A recent work using a gene-tree reconciliation algorithm has revealed the allopolyploid origin of Ericales^[60]. However, no potential progenitors were identified in this study and the sample size was limited. While the proposed view is novel, ongoing genome sequencing projects in Rhododendron and its close relatives will eventually provide clarity on the evolutionary history of this duplication event.

Genes underpinning flower color

Genome sequencing and gene annotation techniques have opened up the possibility of identifying and recognizing genes for flower colors throughout the genome. And multi-omics have opened new perspectives for understanding structural and regulatory genes and their interaction, i.e., the molecular mechanisms underlying flower pigmentation^[61,62]. Anthocyanins have been identified as important pigments influencing color variations in Rhododendron based on transcriptomics and metabolome assessment of floral tissues at different developmental stages^[39,63−66]. Other molecules such as flavonoids and carotenoids are also expected to play a role in pigment production^[67,68]. Nie et al. showed that the yellow flower of R. molle can be attributed to (i) the synthesis of carotenoids and flavonols, (ii) the degradation of chlorophyll, and (iii) the complete absence of anthocyanin synthesis^[26]. The metabolic pathways for carotenoids/chlorophyll/flavonol/anthocyanin were reconstructed in R. simsii and R. molle using whole-genome analysis (Figs 2 & 3)^[26,31]. In these reconstructed metabolic pathways, a total of 15 enzymes related to anthocyanin synthesis were identified, including 11 crucial enzymes (PAL, C4H, 4CL, CHS, CHI, F3H, F3'H, F3'5'H, DFR, ANS, and F3oGT) involved in anthocyanin glycoside synthesis.

Figure 2. 

Metabolic pathways of anthocyanins and flavonols in R. simsii and R. molle. Fractions next to gene abbreviations indicate the ratio of the number of expressed genes for R. molle vs R. simsii; the numerator is the gene number for R. molle and the denominator is the gene number for R. simsii. Equal gene numbers are colored blue, excess ratio in the gene number for R. molle is colored yellow, while excess ratio in R. simsii colored red.

Figure 3. 

Metabolic pathway of (a) chlorophylls and (b) carotenoids in R. simsii and R. molle. Fractions next to gene abbreviations indicate the ratio of the number of expressed genes for R. molle vs R. simsii; the numerator is the gene number for R. molle and the denominator is the gene number for R. simsii. Equal gene numbers are colored blue, excess ratio in the gene number for R. molle is colored yellow, while excess ratio in R. simsii colored red.

Transcription factors (TFs), which influence flower coloration by regulating the expression of pigment-encoding genes, are considered another category of flower color genes. The involvement of transcription factor families such as MYB, bHLH (basic helix-loop-helix), WD40 (WD40-repeat-containing protein), ERF and WRKY in Rhododendron pigment synthesis have been widely reported^[69,70]. Wang et al. identified 68 RsWRKY gene members in R. simsii, which can be divided into 3 main groups and 7 subgroups^[70]. Wu et al. revealed the conserved features during the evolution of R2R3-MYB in Rhododendron^[69]. And Nie et al. identified the members and investigated the evolutionary patterns of five TF families in R. simsii and R. molle, respectively^[26]. These studies have reconstructed the metabolic pathways and identified key regulatory genes, providing insights into the intricate processes that control the development of flower color in Rhododendron. These processes can potentially be used in various fields, including breeding programs to improve floral traits, developing novel cultivars with desired color patterns, and even exploring innovative biotechnological approaches to pigment manipulation in Rhododendron and other woody plants.

Duplication and loss of genes related to flower color

At the whole-genome level, numerous gene duplication events were detected in both R. simsii and R. molle (Fig. 4a), and DSD (dispersed gene duplication) being the most frequent and possibly linked to the expansion of LTR-RTs, leading to differences in gene number between these two Rhododendron species^[26]. Among flower color-related genes, TD (tandem gene duplication) and PD (proximal gene duplication) were identified as the predominant duplication types (Table 2). The gene families associated with anthocyanin synthesis were larger than those involved in carotenoid and chlorophyll synthesis (Fig. 4b). Further analysis suggests that gene duplication, particularly TD and PD, may have facilitated the expansion of these pigment-related gene families. In R. simsii and R. molle, TD and PD accounted for 43.16% and 44.14% of anthocyanin-related genes, respectively. In contrast, these two duplication types accounted for only a small percentage of chlorophyll-related genes (22.22% and 26.56%), and an even smaller percentage of carotenoid-related genes (17.24% and 23.07%) (Table 2). It can be concluded that the expansion of anthocyanin-related gene families, especially TD and PD, in Rhododendrons is responsible for the production of the major pigments. Gene duplication contributes in different ways to the expansion of gene families for flower color. For the R2R3-MYB gene family, WGD contributed to more than 50% of the total number of genes in the three Rhododendron plants^[69], and WGD also contributed to 47.37% (27 of the 57) for WRKY.

Figure 4. 

Gene duplications and losses in R. simsii and R. molle. (a) Barplots show the counts between R. simsii and R. molle of coding genes and pseudogenes for the different duplication modes (WGD, TD, PD, TRD, DSD and SL: whole-genome, tandem, proximal, transposed and dispersed duplications, and singletons). Different suffixes on the left denote '-P' for Pseudogene, '-C' for Coding gene, and '-A' for all coding genes and pseudogenes. (b) Size distribution of enzyme coding gene families for the different pigment metabolic pathways. (c) Pseudogene numbers from the different metabolic pathways for anthocyanin, carotenoid and chlorophyll synthesis. We have linked a pseudogene to a specific enzyme, depending on which enzyme its parental coding gene encodes. (d) Barplots showing the counts of enzyme families from the different metabolic pathways for anthocyanin, carotenoid and chlorophyll synthesis. A fraction > 1 for enzyme encoding gene numbers (R. molle/R. simsii) defines as 'Expansion'; a fraction of 1 defines as 'Stable'; a fraction < 1 defines as 'Contraction'. (Modified from a previous publication^[26]).

The contraction of gene families involved in anthocyanin synthesis is accompanied by pseudogenization in R. simsii and R. molle. Genes related to carotenoid and chlorophyll synthesis underwent relatively little pseudogenization, with only three to four pseudogenes identified in each Rhododendron species (Fig. 4c). In contrast, more pseudogenes were found in anthocyanin-related gene families, notably in R. molle, where 25 anthocyanin-related pseudogenes were discovered (Fig. 4c). Pseudogenization has led to a significant contraction of anthocyanin gene families in R. molle, with 13 contracted families accounting for 81.25% of all 16 families (Fig. 4d). In contrast, there were fewer contracted families and pseudogenes in the chlorophyll and carotenoid metabolic pathways (Fig. 4c & d).

Research to date on flower coloration genes is essentially limited to a handful of species, which is in stark contrast to the intricate interspecific color diversity observed in Rhododendron^[66,71,72]. Interestingly, repetitive sequences could serve as key factors that facilitate the duplication and subsequent loss of color-determining genes. This unfolding perspective enriches our understanding of the evolution of both repetitive sequences and color diversity in Rhododendron. In future, deciphering the intricate evolutionary interplay between LTR-RTs, gene duplications and losses, and color diversity in different Rhododendron species, aligned with phylogenetic frameworks, may become a focus of research.

Gene regulation of flower coloration

The metabolic pathways and regulatory networks that influence phenotypic traits such as color are central to our understanding of plant biology. The genus Rhododendron, one of the most prolific groups of woody plants, aroused great interest in the scientific community. This interest stems mainly from the great diversity of flower colors in wild species and cultivars^[18,73]. So far, research on color differentiation in Rhododendron has primarily investigated the composition and concentration of pigments^[18,20,63]. However, the underlying genetic and regulatory mechanisms that determine these vivid variations remain largely obscure.

In recent years, a number of structural genes responsible for the biosynthesis of flavonoids have been isolated from a variety of wild and cultivated Rhododendron. The spatiotemporal expression patterns of key regulatory genes were analyzed^[13,39,68]. Pigment synthesis has been shown to be under the influence of specific transcription factors including MYB, bHLH, WD40, and their associated complexes, and members of these gene families have been identified and categorized at the whole-genome level in Rhododendron species such as R. simsii, R. molle, R. delavayi, R. irroratum, and R. ovatum^{[29,30,49,60,74]}. In addition, the WRKY and ERF transcription factor families may also play an important role in establishing the temporal gene regulatory network that determines flower color.

The analysis of gene expression networks is an important tool that has been used extensively to study the nuances of flower coloration and differentiation. It is also important to keep in mind that flower color formation is a process that unfolds over time and therefore requires the study of gene expression and regulatory mechanisms at different stages of petal development. We have developed an analytical approach 'InterSpecificTOGCN' (https://github.com/JeffreyNIEgithub/InterSpecificTOGCN) to comprehensively decipher the patterns of the time-ordered gene co-expression networks (TO-GCNs). Here, the approach was applied to analyze the time-series transcriptomes of five timepoints (here T1 to T5) in the process of pigment formation under different Rhododendrons (red flowered R. simsii and yellow flowered R. molle) (Fig. 5a). We described the specific and consensus hierarchical TO-GCNs for two Rhododendrons and reconstructed seven levels (here L1−L7) of TO-GCNs separately (Fig. 5b). Different levels were summarized into three major stages: initial (T1−T2; corresponding to L1−L3, flower appears white and green), transitional (T3; corresponding to L4, during green to yellow/red flower transition), and terminal (T4−T5; corresponding to L5−L7, for the determined yellow or red flower color) (Fig. 5c). The genes of carotenoid/flavonol biosynthesis related to yellow flower specificity formation were thought to be located in the terminal networks (L5 to L7) in R. molle-specific TO-GCNs, based on metabolome, gene function enrichment analysis and network complexity statistics for each stage^[26,31]. We further identified hub genes of SGR (Rhmol07G0115200), PSY (Rhmol10G0271900), PAL (Rhmol03G0128000), and F3'H|F3'5'H (Rhmol09G0029300) families (Fig. 6a), which all represent top graph degrees of R. molle-specific TO-GCNs. By examining the network topology and TF binding site (TFBS) predictions, these four hub genes might be regulated with two TFs, MYB_related (Rhmol11G0261800) and B3 (Rhmol04G0179100), as common regulators in a hierarchical manner in R. molle (Fig. 6b & c)^[26]. In R. simsii, two crucial structural genes, F3oGT (Rhsim01G0008100) and F3'H|F3'5'H (Rhsim09G0023900), and their potential regulators, MYB_related (Rhsim03G0160200) and ERF (Rhsim03G0176600), were predicted within the specifically reconstructed hierarchical subnetworks related to anthocyanins with the same analysis principles (Fig. 6d)^[26]. The genus Rhododendron displays a notable abundance of flower color intermediates, particularly in the case of red hues. It is of great importance to investigate the gene expression changes in the temporal development of these crucial flower color intermediates. Bioinformatics approaches such as TO-GCN have great potential to address these issues from the perspective of time-series gene co-expression network.

Figure 5. 

Time-ordered gene co-expression networks during flower coloration of R. molle and R. simsii. (a) Five timepoints (T1−T5) during flower pigmentation. (b) Predicted regulatory network and the resolved gene regulation among transcription factors (TFs; pink points) and enzymatic genes involved in pathways of carotenoid (orange points), chlorophyll (green points) and anthocyanin/flavonol (red points) biosynthesis. L1 to L7 indicate seven levels identified in three hierarchical time-ordered gene co-expression networks (R. molle-specific, R. simsii-specific and the consensus TO-GCNs). Edges between enzymatic genes were not shown. (c) Heatmaps of average normalized TPMs (z-score) at each timepoint of flowering at each level identified in three TO-GCNs. Three stages of flower coloring were identified as the initial (T1−T2), transitional (T3) and terminal (T4−T5) stages, based on the expression profile. Low to high expression is indicated by a change in color from blue to red.

Figure 6. 

Resolved hierarchical gene regulation during flower coloration of R. molle and R. simsii. (a) Degree of 34 enzymatic genes of the chlorophyll degradation and carotenoid/flavonol biosynthesis in two subnetworks. The degree cutoff for hub gene is 40. (b) Resolved hierarchical regulations for four hub genes PSY, SGR, PAL and F3'H|F3'5'H in R. molle. (c) Gene expressions (TPM) and TF binding site (TFBS) detected in the 2 kb upstream sequences of four hub genes and four potential regulators in R. molle. (d) Resolved hierarchical regulation for two hub genes F3oGT and F3'H|F3'5'H in R. simsii. (Modified from previous publications^[26,31]).

Furthermore, it is hypothesized that LTR-RTs may play an important role in the formation of interspecific gene regulatory patterns. Active LTR-RTs in R. molle could potentially induce flower color shifts by integrating in or near specific pigment biosynthesis genes and transcriptional regulatory factors, triggering changes in gene transcription and regulation^[7]. The applicability and extent of this influence may be difficult to validate in the laboratory; however, their potential impact should not be ignored. As widespread and highly active components in the genome, LTR-RTs have great potential to promote the emergence of the enormous flower color diversity in Rhododendron.

In addition, an Agrobacterium-mediated genetic transformation system for R. simsii has been developed. However, the long growth cycle of woody plants has hindered the experimental validation of genes for flower color in Rhododendron. Therefore, the need to explore innovative methods to validate gene function is particularly evident in the study of Rhododendron color traits. This is not only essential for ongoing Rhododendron-specific research, but is also of great value for the study of reproductive traits in other woody species. The possible application of LTR-RTs as natural mutagens deserves special attention and research.

As an ornamental horticultural plant, Rhododendron has a long history of artificial cultivation and has high visibility and a large industry chain worldwide^[75]. There is a common saying in the Western gardening community: 'Without China's Rhododendrons, there would be no diversity in Western gardens'. Rhododendron breeding faces almost all the challenges of genetic breeding of woody plants, including long cycles, weak breeding foundations, and limited breeding methods^[76]. The selection of flower color is of utmost importance in Rhododendron breeding. Wild Rhododendrons in nature exhibit a wide range of flower colors, while cultivated Rhododendrons are mostly red or white, with fewer yellow and blue varieties that are more valuable. Therefore, breeding yellow or blue Rhododendrons has become an important breeding goal^[77].

From a global perspective, cultivation and breeding studies on Rhododendrons started earlier and more systematically in Japan, North America, and Europe^[6,9]. China is considered the country with the richest resources of wild Rhododendron germplasm, with about 600 of the world's more than 1,200 species^[78]. Collecting and studying germplasm resources is the first step towards the conservation and utilization Rhododendron resources, and previous research has made significant progress in this area^[79−82]. A large number of new Rhododendron cultivars have been bred through traditional breeding methods such as selection and hybrid breeding, which will continue to be the main means of flower breeding in the future. However, advances in biotechnology, especially in molecular genetics and genetic engineering, offer new ideas and new opportunities for Rhododendron breeding and have great market potential.

The genus Rhododendron has undergone extensive interspecific hybridization within broad subgenera and a complex history of radiation evolution^[4,13,25,83]. The improvement in sequence length and accuracy has made genome assembly cheaper, faster, and more accurate^[84]. Functional genomics such as transcriptomics, metabolomics, proteomics, and phenomics have also become more accessible and comprehensive. Advances in these techniques will help researchers to better understand the regulatory mechanisms of gene expression, signal transduction pathways, metabolite synthesis and pathways, protein composition and interactions, and phenotypic variations in Rhododendron^[61,85]. Therefore, genomics-centered next-generation -omics will improve our understanding of Rhododendrons genetic diversity and evolution of in a more precise and efficient way. Due to the weak research base and the nature of Rhododendrons as woody plants, a complete breeding cycle for flower color can take up to 10 years. Flower color is a typical qualitative trait, and population genetics and linkage analysis are effective strategies to identify the genes responsible for flower color. Multi-omics-based breeding will therefore enable faster and more efficient methods for cultivating Rhododendron and other woody plants.

The authors confirm contribution to the paper as follows: study conception and design: Nie S, Yang FS, Mao JF; data collection: Nie S, Ma HY, Tian XC, Shi TL; analysis and interpretation of results: Nie S, Ma HY, Tian XC, Shi TL; draft manuscript preparation: Nie S, Yang FS, Ma HY; edited and improved the manuscript: Porth I, El-Kassaby YA and Mao JF. All authors reviewed the results and approved the final version of the manuscript.

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.