1. Introduction

Roses are well-known ornamental plants that are globally beloved for their highly diverse flowers and flowering times. Flowering at a proper time is a key developmental switch that is essential for reproductive success and the life cycle (Amasino and Michaels, 2010; Baurle and Dean, 2006). Roses feature high diversity of flowering behaviors, including continuous flowering, occasionally re-blooming, and single seasonal blooming; hence, roses serve as a model for studying flowering time diversity in plants (Bendahmane et al., 2013; Dong et al., 2017). Flowering time is under complex but elegant regulation by exogenous and endogenous signaling pathways, which involve a set of gene regulatory networks (Andres and Coupland, 2012; Fornara et al., 2010). In both rose and strawberry, KSN, a homolog of Arabidopsis thaliana TFL1, has been proposed to play an essential role in regulating continuous flowering behavior (Iwata et al., 2012; Randoux et al., 2012, 2014). Duplication and functional diversification of COP1-like genes might also contribute to flowering time regulation in roses as well as other Rosaceae plants (Sun et al., 2020).

Rosa chinensis 'Old Blush' (OB) is a founder genotype in modern rose domestication and hence has been explored for many rose genetic studies (Byrne et al., 2007; Crespel et al., 2002; Hibrand Saint-Oyant et al., 2018; Li et al., 2015, 2019; Shupert et al., 2007; Spiller et al., 2011). However, Rosa wichuraiana 'Basye's Thornless' (BT) harbors important traits that differ from OB. For example, OB plants always grow erect; in contrast, before its annual flowering between April and June, BT shoots switch from prostrate growth in the vegetative stage to erect growth in the reproductive stage. This switch from prostrate to erect growth may help BT adapt to the light requirement and endogenous developmental signals upon flowering transition (Figs. 1A and S1). However, no information is available on the genetic regulatory mechanisms that underlie this switch.

Light is one of the key environmental factors controlling flowering time (Liu et al., 2017, 2020; McCormac and Terry, 2002; Munnik and Nielsen, 2011; Tang et al., 2013). When plants sense a reduction in the ratio of red to far-red light, due to competition for light from the neighbor plants or other types of shade, they initiate a couple of adaptive responses, which are termed shade avoidance syndrome, including rapid shoot elongation, erect shoot, early flowering (Franklin and Whitelam, 2005). Previous studies have shown that FHY3 and FAR1, two homologous transcription factors essential for PhyA-mediated far-red light signaling, can negatively regulate Arabidopsis flowering time under both long-day and short-day conditions by promoting the expression of Early Flowering 3 and 4 (ELF3 and ELF4) (Li et al., 2011; Lin et al., 2007; Liu et al., 2019). Under fluctuating light conditions, FAR1 and FHY3 can alter their physical interaction with SPL3/4/5 transcription factors, which are key regulators in the aging pathway of flowering time control, to modulate the expression of downstream floral integrators such as AP1, FUL, LFY and miR172c, thus modifying the flowering time behaviors of A. thaliana (Xie et al., 2020).

FAR1 and FHY3 regulate a wide range of biological processes, including light signal transduction, photomorphogenesis (Liu et al., 2019; Siddiqui et al., 2016; Tang et al., 2013; Wang et al., 2016; Zhang et al., 2019), circadian clock and flowering time (Johansson and Staiger, 2015; Li et al., 2011; Liu et al., 2020; Ritter et al., 2017; Xie et al., 2020), shoot and floral development (Li et al., 2016), chloroplast and chlorophyll biosynthesis (Ouyang et al., 2011; Tang et al., 2012; Wang et al., 2016), starch synthesis (Ma et al., 2017), ABA (Tang et al., 2013) and oxidative stress responses (Ma et al., 2016) as well as immunity (Wang et al., 2016). Arabidopsis also has several FAR-RELATED SEQUENCE (FRS) and FRS-RELATED FACTOR (FRF) proteins that can modulate flowering time (Bulik-Sullivan et al., 2015; Gao et al., 2013; Ma and Li, 2018; Ritter et al., 2017; Tang et al., 2013). However, the conservation and diversification of FAR1/FRS-like genes in other species, especially in woody plants, such as roses, remains to be explored in detail.

In this study, we identified FAR1/FRS-like genes in two rose genotypes, BT and OB, for which high-quality whole-genome sequences are available (Hibrand Saint-Oyant et al., 2018; Raymond et al., 2018; Zhong et al., 2020). We found that several gene duplication mechanisms, including whole-genome duplication (WGD), segmental duplication, tandem duplication and transposon element (TE) mediated duplication, especially the dispersed duplication were simultaneously involved in the expansion of FAR1/FRS-like genes in BT. Finally, we propose that the dispersed FAR1/FRS-like genes expansion and their expression pattern might associate with the regulatory of shoot-growth behavior in BT.

2. Materials and methods 2.1. Plant materials and transcriptomic profiling

The high-quality chromosome-level genome assembly for the BT genotype was constructed by ourselves (Zhong et al., 2020), and the genome sequences for haplo-OB was used for comparisons (Raymond et al., 2018). BT and OB plants were grown in the glasshouses at the Flower Research Institute of Yunnan Academy of Agricultural Sciences (Kunming, Yunnan, China). Leaf materials were collected in November (Nov) and March (Mar). Shoot apical meristem materials (SAM) with minimum leaf materials were sampled in March. RNA extraction followed by strand-specific Illumina sequencing has been described previously (Li et al., 2018).

2.2. Identification and phylogenetic analysis of FAR/FRS-like genes in roses

Protein amino acid sequences were downloaded for FAR1/FRS-like genes from TAIR10 (https://www.arabidopsis.org/), and were then used for BLAST identification of homologous sequences in BT and OB genomes with iTAK pipeline (-f 6) (Zheng et al., 2016). MAFFT was applied for sequence alignment with default parameters (Katoh et al., 2002). Phylogeny reconstructions were carried out using the Neighbor-Joining (NJ) method and 100 replicates of bootstrap simulation in RAxML 8.2.11 (stamatakis2014">Stamatakis, 2014). Chromosome distributions and annotations were plotted by TBtools (Chen et al., 2020) with house R scripts. Gene structures were plotted using Gene Structure Display Server 2.0 to show their exon/intron composition information (Hu et al., 2015) (http://gsds.cbi.pku.edu.cn//index.php).

2.3. Duplications and syntenic analysis of FAR1/RS-like genes

MCSanX and related functional blocks were used to compare chromosome structural variation (Wang et al., 2012). The FAR1/FRS duplication type and collinearity were predicted with the duplicate_gene_classifier function according to protocols described in the pipeline manuals. The phylogenetically clustered gene sets were defined as syntenic only when a minimum of five genes were collinear. The remaining genes were considered non-collinear. The non-collinear genes without phylogenetic clustering signal were designated as BT expansion genes.

2.4. Exon-intron structure, conserved protein domains and motif analyses

Conserved protein domains were analyzed with the NCBI CD-Search Tool (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) using the PSSM model (maximum number of hits < 500). Conserved motifs in FAR1/FRS-like genes were identified with MEME (http://meme-suite.org/tools/meme; optimum motifs, 6-15; maximum number of motifs, 8). TBtools was used for exon-intron visualization.

2.5. Identification of insertion events between BT and OB

To find large insertion (100–50, 000 bp) events in BT compared to OB genome, each BT scaffold sequence was aligned to the OB genome with blastn. Only alignments with ≥60% identity and longer than 100 bp were kept for further analysis. Insertion density was plotted as the insertion events every 100 kb in 1 Mb window with ggplot2 in R with home scripts.

2.6. Transposon element (TE) analysis

Transposon elements were analyzed with RepeatMasker (V4.1.0) (Saha et al., 2008) using plant RepBase database (https://www.girinst.org/) and Dfam 3.1 database (www.dfam.org). TEs located within 2000 bp up- and down-stream of the FAR1/FRS-like genes were counted and compared to whole genome level counts per protein-coding gene. Statistical differences between the expansion and syntenic FAR1/FRS-like genes, and the rest of the whole genome level was tested with the functionWilcoxon rank sum test in R.

2.7. Expression analysis for FAR1/FRS-like genes and their potential interactor encoding genes

RNA-seq reads for shoot apical meristem (SAM), just opened young leaves in March (leaf_Mar) and November (leaf_Nov), which had been described previously (Li et al., 2018), were mapped to BT and OB genomes with HISAT2 (Kim et al., 2015) and assembled with StringTie (Pertea et al., 2015). Uniquely mapped reads (default parameters expect for -stranded = no) were used to calculate the relative expression as fragments per kilobase per million (FPKM) with ballgown package in R, and further compared for the expression pattern in different tissues and developmental stages. Normalization was carried out using the ZeroToOne method and scaled with log2. We collected proteins that interact with FAR1/FRS in A. thaliana with STRING online (https://www.expasy.org/). Protein sequences of these genes were extracted and used as seed sequences to blastp the BT and OB genomes (e-value 10^-7, minimum identity 0.35, minimum length coverage 0.6). Ortholog gene pairs were then identified using OrthoFinder2 (emms2019">Emms and Kelly, 2019) in combination with the blastp, phylogeny, and synteny-based manual correction.

3. Results 3.1. BT expands significantly its FAR1/FRS-like genes than OB

Via sequence similarity analysis to the 17 Arabidopsis FAR1/FRS family genes, we identified 91 and 52 FAR1/FRS-like genes in BT (Zhong et al., 2020) and OB (Raymond et al., 2018), respectively (Figs. 1B and S2; Tables S1 and S2). In BT, these genes were almost randomly dispersed on the seven chromosomes, with Chr2, 5, 6 harboring 17, 15, 16 genes, respectively. BT expanded significantly the FAR1/FRS-like genes, especially on Chromosomes 4, 5, 6, in comparison to OB genotype as well as other Rosaceae plants (Yates' chi-square test, P < 0.001; Fig. 1B; Table S2). Gene-collinearity and phylogenetic reconstruction analysis using the OB genotype as a reference indicated that 31 of the FAR1/FRS-like genes in BT are the result of gene expansion events (Fig. 1C).

3.2. Random duplication dominates BT expansion of RwFAR1FRS-like genes

Next, we examined the potential mechanisms underlying the FAR1/FRS-likes expansion in BT. We classified these genes into four types of duplication events: whole-genome-duplication (WGD; two events), tandem duplication (nine events), proximal duplication (six events) and dispersed duplication (69 events). For the apparently randomly-duplicated gene sets, the gene structure and organization, including untranslated regions (UTRs), exons and introns, differed significantly from each other (Fig. S3). This pattern indicates that these were real tandemly duplicated gene sets not by products of erroneous annotation or assembly. The detection of identical numbers of WGD and tandem duplication events for OB suggests that random dispersion played a dominant role (29 among 31) in the expansion of BT RwFAR1/FRS-like genes (Fig. 1D; Table S3). However, random insertion of 100–50, 000 bp DNA fragments did not correlate with the expansion.

3.3. Dispersed transposon insertion accompanies the expansion of RwFAR1/FRS-likes

Because FAR1/FRS-like genes are derived from transposases, we next asked whether expansion of RwFAR1/FRS-like genes in BT is correlated with transposon element distribution (Hudson et al., 1999; Lin et al., 2007; Lin and Wang, 2004). We counted the TE insertion events within two kilobases up- and down-stream of BT expansion (non-syntenic) and syntenic RwFAR1/FRS-like genes. TE distance in expansion RwFAR1/FRS-like genes was not significantly different than in other genes throughout the genome (Fig. S4). Interestingly, the number of TEs was significantly higher in the expansion RwFAR1/FRS-like genes than in both the syntenic RwFAR1/FRS-like genes and in other genes in the genome (Fig. 1E; Table S4). This increased number of TEs indicates that the expansion of RwFAR1/FRS-like genes may have been accompanied by transposon activity.

3.4. Expanded RwFAR1/FRSs feature diversified protein structures and motifs

We examined the characteristics of proteins encoded by the genes of the RwFAR1/FRS-like expansion. The gene structure of syntenic and expansion RwFAR1/FRS-like genes (e.g., exon-intron numbers and organization) were very similar (Figs. S5 and S6). Motif analysis revealed that proteins encoded by expansion RwFAR1/FRS-like genes have between two to eight motif types, similar to proteins of syntenic RwFAR1/FRS-like genes (Fig. S7). Proteins encoded by expansion RwFAR1/FRS-like genes contain FAR1, FHY3 and DDE_Tnp_ISL3 domains individually or in combination. These patterns of protein motifs do not differ from those of proteins encoded by syntenic RwFAR1/FRS-like genes (Fig. S6), which corroborates the hypothesis that expansion RwFAR1/FRS-like genes were dispersed via random transposon insertion.

3.5. Lineage-specific expansion of RwFAR1/FRS-like genes mainly occurs in clade I and II

To understand whether the RwFAR1/FRS-like gene expansion had any further features, we performed a phylogenetic analysis with protein sequence alignment. This analysis grouped the FAR1/FRSs into five clades (I–V) with all Arabidopsis proteins clustered in clades III-V and both FAR1 and FHY3 sitting within clade IV (Fig. 2A). Both BT and OB had similar numbers of FAR1/FRSs to Arabidopsis within clades III-V. However, BT (54) doubled its FAR1/FRSs in clades I and II in comparison to OB (25; Yates' chi-square test, P < 0.001) with most of them randomly dispersed on the seven BT chromosomes (Fig. 2A). Twenty-three among the 31 (74.2%) expanded RwFAR1/FRS-like genes were grouped in clades I and II (Fig. 2A), indicating that BT featured a lineage-specific expansion of RwFAR1/FRS-like genes.

3.6. RwFAR1/FRS-like genes are expressed at relatively higher and diverged levels than OB ones

In BT, several genes showed a relatively high expression (FPKM values more than 30) in SAM (13), March leaves (8) and November leaves (9). No OB genes were expressed at relatively high levels. Several OB genes were expressed at medium levels in SAM (23), March leaves (1), and November leaves (1); in contrast, there were a greater number of BT genes expressed at medium levels in SAM (24), March leaves (18), and November (18) (Fig. 2B; Table S5). Three BT genes (Rw6G036920 in clade I, Rw1G014580 in clade III, and Rw3G000170 in clade V) were only expressed in leaves (Fig. 2B; Table S5). Notably, BT paralogs were always expressed at relatively higher levels than OB genes in those expressed FAR1/FRS-like genes. Additionally, the BT expansion genes showed a variety of expression patterns, including expression restricted to either SAM or leaves, expression in both tissues, or no expression (Fig. 2A). Only two genes (Rw2G034530 in clade II and Rw3G000160 in clade V) showed differential expression between March and November leaves, indicating that these genes might have specialized function. Taken together, these expression patterns indicate that FAR1/FRS-like gene expression differs significantly in BT and OB.

3.7. Expression of FAR1/FRS-like protein interaction partners differs between BT and OB

Arabidopsis FAR1 and FHY3 have several known interaction partners and their interactions play essential roles in light-signaling processes (Fig. S8). However, among these potential interacting partners, only homologous genes encoding SPL4/5-like proteins showed differential expression between BT and OB in different tissues (Fig. 2C; Fig. S9; Table S6). Because the molecular interaction between FAR1/FHY3 and SPL4/5 plays an important role in flowering time regulation in Arabidopsis (Xie et al., 2020), we hypothesize that the coordinated variation in expression between RwFAR1/FRS-like genes and RwSPL4/5-like genes might play a role in mediating the switch from prostrate-to-erect growth in shoots, and therefore flowering time in BT plants (Fig. 2D).

4. Discussion

In this study, transposon element (TE)-mediated duplication was detected, a pattern has been reported in many other plants (Cannon et al., 2004; Dodsworth et al., 2016; Flagel and Wendel, 2009; Li et al., 2017; Panchy et al., 2016; Qiao et al., 2019). We found that FAR1/FRS-like genes diversified in their gene structure, motifs, and protein domain compositions after TE-associated dispersed expansion. More interestingly, we observed strong variation in gene expression between rose genotypes for this important family of transcription factors. Previous studies have shown that various FAR1/FRS-like genes have different DNA-binding activities and functional features involved in regulating plant development and environmental adaptation (Li et al., 2016; Ma and Li, 2018; O'Malley et al., 2016; Ritter et al., 2017); thus, we speculate that variation in gene expression may account for shoot growth differences and flowering time changes between rose genotypes. Confirmation of this hypothesis requires further functional characterization. FAR1/FRS-like genes play important roles in adaptation to varying light and other environmental conditions in A. thaliana, but the evolutionary pattern of this important family of transcriptional regulators has been poorly investigated especially in woody plants. Our study thus represents one of the first examples linking FAR1/FRS-like gene evolution with mechanisms in the generation and maintenance of plant diversity (Grierson et al., 2011; Pennisi 2005).

In conclusion, the use of high-quality chromosome-level genome sequences (Hibrand Saint-Oyant et al., 2018; Raymond et al., 2018; Zhong et al., 2020) has increased our understanding of the evolutionary patterns of key transcription factors, the FAR1/FRS-like family. Our findings provide an important starting point for dissecting the molecular genetic regulation of many morphological novelties in non-model plants, especially woody plants (Dong et al., 2017). Combining data described in this research with further genetic analyses (e.g., QTL and GWAS technologies) may also offer important information for marker-assisted breeding of non-model crops with interesting traits.

Author contributions

J-Y H and M-C Z conceptualized the project; X-D J and W–H C collected the samples, extracted the genomic DNA and total RNA. M-C Z analyzed and visualized the data. J-Y H and M-C Z drafted the manuscript with contributions from all authors. All authors have read and approved the final manuscript.

Declaration of competing interest

The authors declare no competing interest.

Acknowledgements

This work was funded by the Strategic Priority Research Program of the Chinese Academy of Sciences to J-Y H (XDB31000000); the CAS Pioneer Hundred Talents Program to J-Y H (292015312D11035) and Yunnan Recruitment Program of Experts in Science to J-Y H. We acknowledge the support from Dan Wang and Yi-Bo Sun. We appreciate Shubin Li, Shulan Chen, Hongyuan Yu and Yuanlin Lv for the help with rose cultivation. This work is partially facilitated by the Germplasm Bank of Wild Species of China and Kunming Botanical Garden, Kunming Institute of Botany, Chinese Academy of Sciences.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2020.11.002.