To identify the PEBP genes in various plant species, we performed BLASTP (E-value<1.0e-15) search against genomes of Arabidopsis thaliana, Spinacia oleracea, Oryza sativa and Beta vulgaris in Phytozome V13 (https://phytozome-next.jgi.doe.gov), and C. pallidicaule and C. suecicum in Chenopodium database (https://www.cbrc.kaust.edu.sa/chenopodiumdb/), using FT (AT1G65480.1), TFL (AT5G03840.1) and MFT (AT1G18100.1) protein sequences of Arabidopsis as the queries. Then those homologs were aligned with the Hidden Markov Model (HMM) profile of PEBP domain (PF01161) using Pfam search (E-value<1.0e-10) (Finn et al., 2016) to ensure those homologs habor a PEBP domain. Multiple sequence alignment of PEBP sequences from various species was performed using CLASTALW (Thompson et al., 2002). Phylogenetic tree was constructed using MEGA 11.0 (Tamura et al., 2021) based on the Neighbor-Joining method (Dohm et al., 2013) with a bootstrap value of 1000.

The General Feature Format (GFF) file and chromosome-scale genome sequence of quinoa were downloaded from Phytozome V13 database (https://phytozome-next.jgi.doe.gov/info/Cquinoa_v1_0). Based on these two files, the physical location on chromosomes, and intron and exon structures of PEBP genes were determined and visualized by using the two programs of Gene Location Visualize and Gene Structure View in TBtools (Chen et al., 2020). Multiple Em for Motif Elicitation (MEME) program (https://meme-suite.org/meme/tools/meme) was used to identify the conserved motifs in PEBP proteins setting the maximum motif count of 8. The motif analysis results were illustrated using the Gene Structure View program in TBtools (Chen et al., 2020).

The GFF file and chromosome-scale genome sequence of sugar beet (Beta vulgaris) were downloaded from Phytozome V13 database (https://phytozome-next.jgi.doe.gov/info/Bvulgaris_EL10_1_0). The GFF files and chromosome-scale genomes of quinoa and sugar beet were input into the Multiple collinear scanning tool (MCScanX) in TBtools (Chen et al., 2020) and the collinearity between quinoa and sugar beet genomes was analyzed. Then the collinearities between PEBP genes from quinoa and sugar beet were determined. Meanwhile, the duplication events among quinoa PEBP genes were analyzed. The results were visualized by using the Advanced Circos program in TBtools (Chen et al., 2020).

The 2000bp sequences upstream of the translation start site of PEBP genes were recognized as promoter regions. The promoter sequences were uploaded to the Plant Cis-Acting Regulatory Elements (PlantCARE) database to search cis-acting elements. The physical distribution of various cis-acting elements was displayed using Simple BioSequence Viewer in TBtools (Chen et al., 2020).

For tissue-specific expression pattern analysis, we used the RNA-seq data from Zou et al. (Zou et al., 2017). Raw data of various quinoa tissues, including 1-week-old seedling, stem, leaf, inflorescence from 6-week-old plants and dry seed, was downloaded from Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) of the BioProject PRJNA394651. The fragments per kilobase of transcript per million fragments mapped (FPKM) value of each gene was calculated with the previously described methods (Wu et al., 2019; Wu et al., 2020; Wu et al., 2021). For gene expression analysis in quinoa inflorescences at six developmental stages, we investigated the FPKM values of PEBP genes in our published transcriptome data (Wu et al., 2019). YP1, YP2, YP3 and YP4 represent the young non-branching panicles, whereas P1 and P2 stand for the panicles of elder branching stages. Raw data generated from six-stage inflorescence samples was downloaded from SRA (BioProject PRJNA511332), and was further processed with the bioinformatic pipeline methods described before (Wu et al., 2020). To know gene expression changes of PEBP genes during seed germination, we used the published transcriptome data of our laboratory (Wu et al., 2020). We compared the expression levels (FPKM value) of PEBP genes in control and Abscisic acid (ABA)-treated seeds 5h and 15h after imbibition. To investigate PEBP gene expression changes in response to night-break (NB) treatment, we used our published RNA-seq data generated from leaf samples treated by short-day and NB conditions (Wu et al., 2021). To investigate the PEBP gene diurnal expression pattern changes in response to photoperiod, two-week-old quinoa seedlings were cultivated under short-day (8h/16h, light/dark) and long-day (16h/8h, light/dark) conditions for sampling and RNA-seq analysis. The top two fully-expanded leaves of 3~5 plants were harvested at the time point of 17:00, 20:00, 23:00, 02:00, 05:00, 08:00, 11:00 and 14:00, with two biological replicates. All the leaf samples were subjected to RNA extraction, high-throughput sequencing and data analysis as previously described (Wu et al., 2019; Wu et al., 2021).

The expression profiles of 32 samples covering 16 time points of SD and LD were subjected to weighted gene co-expression network analysis (WGCNA) by applying R package, with the parameters set as following: gene expression threshold: FPKM≥1.0, power=1, minimum module size=30, minimum height for merging modules=0.25. As a result, 3972 genes were sorted into 6 co-expressed modules. CqFTL5 and CqFTL8 with were clustered into the blue module containing 934 co-expressed genes. Then all the 934 protein sequences in blue module were uploaded to PlantTFDB v5.0 website (Jin et al., 2017) (http://planttfdb.cbi.pku.edu.cn) for TF prediction.

By using BLASTP and Pfam search methods, 23 PEBP homologous genes were identified in quinoa genome. The shortest quinoa PEBP gene (AUR62033889, named CqFTL3), encoding for 88 amino acid residues ( Table 1 ), was predicted to harbor an incomplete PEBP domain. The longest PEBP gene (AUR62013052, named CqFTL2), encoding for 339 amino acid residues ( Table 1 ), was predicted to harbor two PEBP domains. These results are in line with the points in previous study (Jarvis et al., 2017; Štorchová, 2020), indicating that CqFTL3 may be a pseudogene and CqFTL2 may experience tandem duplication.

To evaluate the phylogenetic relationships between PEBP genes of quinoa, and closet relatives of its diploid ancestors, C. pallidicaule (AA) and C. suecicum (BB), and its relatives in Amaranthaceae family, spinach (Spinacia oleracea) and sugar beet (Beta vulgaris), and model dicot and monocot plants, Arabidopsis (Arabidopsis thaliana) and rice (Oryza sativa), a Neighbor-Joining phylogenetic tree containing 83 sequences ( Supplementary File 1 ) was inferred using MEGA 11.0 (Tamura et al., 2021). As displayed in Figure 1 , 23 quinoa PEBP genes were sorted into three major subfamilies. Quinoa contains 5 MFT clade members, named CqMFT1 to CqMFT5, 11 FT-like (FTL) clade members, named CqFTL1 to CqFTL11, 7 TFL1-like (TFL) clade members, named CqTFL1 to CqTFL7 ( Figure 1 ). C. pallidicaule contains 3 MFTs, 4 FTLs and 3 TFLs, and C. suecicum contains 2 MFTs, 5 FTLs and 3 TFLs ( Figure 1 ).

Then, we drew a chromosomal location map of PEBPs. As illustrated in Figure 2 , 23 PEBP genes are distributed on 12 chromosomes with the exception that CqTFL6 is anchored on Chr00. Sub-genome A and sub-genome B have nearly equal number of PEBPs. Chr06 of sub-genome B and Chr07 of sub-genome A contain the largest number of CqFTLs and CqTFLs, respectively.

Diagram of PEBP gene structures shows that, out of 23 PEBPs, most (17 of 23) harbor 4 exons and 3 introns ( Figure 3 ). CqFTL2 contains the most (8 exons and 7 introns), whereas CqFTL3 contains the least (2 exons and 1 intron) number of exons and introns ( Figure 3 ). MEME program was used to identify the conserved motifs in quinoa PEBP proteins. A total of 8 conserved motifs were identified in PEBP proteins ( Figure 3 ). Motif composition schematic analysis indicates that motif 1, 2, 3, 4 and 5 were constitutively occurred in most of the PEBP genes ( Figure 3 ). Yet, several PEBP proteins lack some specific motifs. For example, CqFTL2 lacks motif 1, 2 and 3; CqFTL8, CqFTL9 and CqFTL10 lack motifs 2 and 3; CqFTL11 lacks motif 4; CqMFT1 lacks motif 5; CqMFT4 lacks motifs 2, 4, and 5 ( Figure 3 ). In addition, we found that motif 6 was occurred in most of the FTL and TFL subfamilies, whereas was absent in MFT subfamily ( Figure 3 ). Motif 7 instead of motif 8 was possessed by most of the FTL subfamily members, while motif 8 was exclusively occurred in majority of TFL and MFT subfamilies ( Figure 3 ). As PEBP proteins usually interact with other proteins (Taoka et al., 2013; Tsuji et al., 2013), such as FD and 14-3-3 proteins, to exert their roles, we speculated that these distinct motif compositions of quinoa PEBPs may result in varied protein-protein interactions and function diversities.

Multiple sequence alignment suggests that the conserved amino acid residue tyrosine (Y) or histidine (H) in motif 1 is a key site distinguishing FTL and TFL clades ( Figure 4 ), in consistent with the findings in PEBP families of other plants (Wang et al., 2015). Besides, as the occurrences of motif 7 and motif 8 in FTL and TFL/MFT subfamilies are mutually exclusive ( Figure 4 ), we speculated that in these two motifs there may be specific amino acid residues distinguishing FTL and TFL/MFT clades. As expected, we found that, the third amino acid residue in motif 7 and 8 was Glycine (G) in FTL genes, whereas it was Alanine (A) in TFL and MFT genes ( Figure 4 ). Thus, this key site (G/A) may provide as a novel site for investigation of the converse functions of FTL and TFL/MFT clades.

Diploid sugar beet contains 8 PEBP genes in the genome. As a close relative of sugar beet, theoretically tetraploid quinoa (AABB) should has doubled number of PEBP homologs. However, we identified 23 PEBPs, nearly 3 times of that in B. vulgaris ( Table 1 ). We wondered what caused the disproportional PEBP family expansion. Gene number ratios between quinoa and sugar beet in MFT and TFL clades are 5:2 and 7:3, respectively ( Table 1 ). Surprisingly, gene number ratio in FTL clade is 11:3 ( Table 1 ), remarkably higher than the plausible ratio 2:1. We analyzed the collinearities of PEBP genes between sub-genomes, and found that a total of 7 orthologous gene pairs from sub-genome A and B, including CqMFT4/CqMFT5, CqFTL1/CqFTL2, CqFTL4/CqFTL5, CqFTL6/CqFTL9, CqTFL1/CqTFL3, CqTFL4/CqTFL5 and CqTFL6/CqTFL7, displayed collinear relationships ( Table 1 ; Figure 5A ). In addition, CqTFL1 and CqTFL5 also displayed collinearity ( Figure 5A ), indicating inner sub-genome duplication. Collinearity analysis between quinoa and sugar beet genomes showed that, 13 out of 23 PEBP genes of quinoa were orthologous to 7 out of 8 PEBPs of sugar beet ( Table 1 ; Figure 5B ). Due to the allotetraploidy of quinoa genome, most of the sugar beet PEBPs have two sister genes, from quinoa sub-genome A and sub-genome B, respectively. There are 7 orthologous clusters, including CqMFT4/CqMFT5-BvMFT1, CqFTL1/CqFTL2-BvFT1, CqFTL4/CqFTL5-BvFT2, CqFTL6-BvFTL3, CqTFL1/CqTFL3-BvTFL1, CqTFL1/CqTFL4/CqTFL5-BvTFL2 and CqTFL6/CqTFL7-BvBFT1 ( Table 1 ; Figure 5B ). Most of these orthologous groups are in line with the gene pairs from sub-genome A and sub-genome B. There are 3 copies of quinoa TFL gene in the CqTFL1/CqTFL4/CqTFL5-BvTFL2 cluster ( Figure 5B ), probably rising from the segmental duplication events between Chr17 and Chr18 in sub-genome B ( Figure 5A ).

The rest of PEBPs, including CqMFT1, CqMFT2, CqMFT3, CqFTL3, CqFTL7, CqFTL8, CqFTL10, CqFTL11 and CqTFL2, lack syntenic regions. We analyzed the gene physical locations, and found some of them were distributed in close vicinity to the other PEBP genes ( Table 1 ; Figure 2 ). Among those genes, 4 tandem repeats were found, including CqMFT3/CqMFT4 on Chr01, CqMFT2/CqMFT5 on Chr02, CqFTL9/CqFTL10/CqFTL11/CqFTL7 on Chr06 and CqTFL2/CqTFL3 on Chr07 ( Table 1 ; Figure 2 ). Notably, all these tandem repeats are located on distal telomeric ends of chromosomes ( Figure 5A ), indicating possible higher frequency of tandem duplications in distal telomeric segments. In general, these evidences indicated that tandem duplication is the major mechanism underlying PEBP family expansion.

Then, we analyzed the rise of gene duplication of PEBPs by combining phylogenetic analysis with collinearity analysis. According to the genome map of sugar beet, BvMFT1 and BvMFT2 are in tandem repeat location. Meanwhile, we found the orthologs, in the relatives of sub-genome A (AAA07574/AAA07573 in C. pallidicaule) and sub-genome B (BBB10201/BBB10201 in C. suecicum), are also in tandem repeat locations. Thus, we deduced that the tandem duplications of CqMFT3/CqMFT4 and CqMFT2/CqMFT5 happened before the divergence of the ancestor of Chenopodium from sugar beet, far before the tetraploidization of quinoa. Though CqMFT1 was in non-synteny with sugar beet and C. suecicum, however, it is sister to the MFT genes of Arabidopsis, and C. pallidicaule and spinach. Phylogenetic analysis indicated that the tandem duplication of CqTFL2/CqTFL3 was happened after the ancestor of sub-genome A diverged from C. pallidicaule. CqFTL3, a sister of CqFTL4 in the phylogenetic tree, but is in non-syntenic region and lacks 4 motifs of CqFTL4, indicating that gene amplification may occur through transposable elements. CqFTL8, despite located nearby CqFTL6, but was not detected as a tandem repeat with CqFTL6 ( Figure 5A ). Meanwhile, CqFTL8 is in non-syntenic region and lacks motif 3, 6, 2, 7 compared with the other CqFTL genes ( Figure 3 ), indicating it might arise from small-scale transposition. Synteny was found in CqFTL9 of sub-genome A and CqFTL6 of sub-genome B ( Figure 5A ), whereas only CqFTL6, rather than CqFTL9, has synteny with BvFTL3 of sugar beet ( Figure 5B ), suggesting the tandem duplication in CqFTL9/CqFTL10/CqFTL11/CqFTL7 possibly was generated after the divergence of the ancestor of sub-genome B from sugar beet.

Cis-acting elements in gene promotor have important regulatory roles mediating transcriptional activation and repression, and various cis-acting elements controlling specific progresses have been identified. To predict and understand the versatile functions of PEBP genes, we explored the cis-acting elements in promotors. Two kilo base pairs upstream of the translational start site of PEBP genes were submitted to PlantCARE database for cis-acting elements prediction. As displayed in Figure 6 , various cis-acting elements, related to light, phytohormones, cold, drought and circadian clock responses, were identified. Light responsive elements (LREs) take a great proportion among various cis-acting elements in PEBP promoters, suggesting that PEBPs may be tightly associated with light biological events. A total of 8 kinds of LREs were identified. LRE_Box4 and LRE_G-box were the two major elements distributed in all PEBP promoters. Some cis-acting elements showed gene-specific distribution patterns. CqFTL8 promoter harbors the most abundant LRE_Box4. CqFTL11 promoter harbors more LRE_Sp1. LRE_G-box rather than LRE_Box4 was more abundant in the promoters of MFT clade. More low-temperature responsive elements (LTRs) were enriched in CqMFT1 promoter, indicating the possible roles of CqMFT1 in sensing cold stress. Circadian regulatory elements (Circadian) were only distributed in the promoters of CqFTL7, CqFTL11, CqMFT2, CqTFL2 and CqTFL3. Abscisic acid responsive elements (ABREs) were abundantly located in promoters of CqMFT2 and CqMFT3, indicating these two gene may be involved in ABA signaling. In addition, we noticed distributions of the cis-acting elements in promoters of some gene pairs were also in synteny. The distributions of ABRE, AuxRR-core, LRE_Box4, LRE_G-box, LRE_GT1-motif, LRE_Sp1 and LRE_TCT-motif in CqFTL1 promoter were in synteny with that in CqFTL2 promoter. This is in line with the collinearity between these two genes. The distributions of the cis-acting elements including Circadian, GRE_P-box, LRE_Box4, LRE_chs-CMA1a and LRE_GT1-motif in CqTFL2 promoter were in highly collinear relationship with that in CqTFL3 promoter. This may be caused by the tandem duplication of CqTFL2/CqTFL3. The specific distributions of cis-acting elements in PEBP promoters suggest their possible roles in different biological events.

To characterize the expression patterns of PEBPs in various quinoa organs, we analyzed the RNA-seq data and compared the FPKM values of 23 PEBP genes in seedling, leaf, stem, inflorescence and seed. Out of 23 PEBP genes, 18 were expressed in at least one organ tissue ( Figure 7 ; Supplementary File 2 ). Generally, expressions of PEBPs in various tissues were clade-specific. FTL clade genes were more enriched in leaf, stem and inflorescence, and TFL clade genes were abundant in seedling, stem and inflorescence, whereas expressions of the MFT clade genes were relatively higher in seed ( Figure 7 ). It seems like that expressions of the MFT clade genes are more conserved, and all the MFT gene pairs were expressed in similar patterns ( Figure 7 ). By contrast, expressions of the orthologous gene pairs in FTL and TFL clades were likely differing ( Figure 7 ). For CqFTL1/CqFTL2, CqFTL1 was highly enriched in inflorescence, whereas CqFTL2 was ubiquitously expressed in all tissues ( Figure 7 ). For CqFTL4/CqFTL5, CqFTL4 was expressed mainly in leaf, and CqFTL5 was abundant in leaf and stem ( Figure 7 ). Another branch containing CqFTL7 and CqFTL8 was expressed specifically in leaf ( Figure 7 ). For CqTFL1/CqTFL3, higher expressions of CqTFL1 were detected in both stem and inflorescence, and CqTFL3 was specifically abundant in stem ( Figure 7 ). The duplicated gene of CqTFL2, CqTFL3, was also expressed mainly in stem ( Figure 7 ). For CqTFL6/CqTFL7, CqTFL6 was specifically expressed in seedling, whereas CqTFL7 transcripts were abundant in seedling and inflorescence ( Figure 7 ). CqTFL4/CqTFL5 harbored similar expression patterns, both were specifically expressed in seedling. The differing expression patterns in gene pairs indicate the possible diversified roles of orthologs from sub-genome A and sub-genome B.

Branching of inflorescence is a key factor influencing plant architecture and yield. Numerous evidences have suggested the important roles of PEBP genes in inflorescence branching and yield control (Teo et al., 2014; Benlloch et al., 2015). To investigate which quinoa PEBP members are potential regulator for panicle branching, we analyzed the expression changes of PEBPs across 6 developmental stages, before and after panicle branching. Expressions of 19 PEBPs were detected in at least one stage ( Figure 8 ; Supplementary File 2 ). Differing from the diversified expressions in tissues above, most of the orthologous gene pairs in FTL and TFL clades showed similar expression patterns. For FTL clade, with the exception that CqFTL3 and CqFTL4 were expressed ubiquitously at non-branching stages (YP1 to YP4) ( Figure 8 ), expressions of CqFTL1/CqFTL2, CqFTL5, CqFTL6/CqFTL9 and CqFTL8, ascended with the development of inflorescence and were abundant at branching stages (P1 and P2) ( Figure 8 ). For TFL clade, expression levels of gene pairs CqTFL1/CqTFL3, CqTFL4/CqTFL5, and the duplicated gene of CqTFL3, CqTFL2 were relatively higher at non-branching stages (YP1 and YP2) while descended with the development of inflorescence ( Figure 8 ), whereas the gene pairs CqTFL6/CqTFL7 were relatively abundant in branching stages (P2) ( Figure 8 ). Most of the MFT clade genes were ubiquitously expressed across all stages ( Figure 8 ). The specific expressions of CqTFL1/CqTFL2/CqTFL3 and CqTFL4/CqTFL5 at non-branching stages, and the enrichment of CqFTL1/CqFTL2, CqFTL5, CqFTL6/CqFTL9, CqFTL8 and CqTFL6/CqTFL7 at branching stages, suggest that those PEBPs may participate in panicle architecture regulation.

Pre-harvest sprouting (PHS) is a knotty problem that influences the yield and nutritional qualities of quinoa (Wu et al., 2020). Understanding the mechanisms underlying seed germination will facilitate breeding PHS-resistant elites of quinoa. We investigated the expressions of PEBPs during seed germination before and after ABA treatment. A higher proportion of MFT clade genes (4 out of 5, 80%), whereas relative lower percentages of FTL clade (6 out of 11, 55%) and TFL clade genes (1 out of 7, 15%) were detected during seed germination ( Figure 9 ; Supplementary File 2 ). We found CqMFT2, CqMFT3 and CqMFT4 were largely repressed as germination went on ( Figure 9 ). Further, we found that the transcriptional changes of CqMFT2, CqMFT3 and CqMFT4 during germination were attenuated when treated by ABA ( Figure 9 ). This is in agreement with the enrichment of the cis-acting element ABRE in MFT gene promoters, indicating those three genes may respond to ABA to regulate seed germination.

Night break (NB) has a repressive effect on the flowering of short-day plant (SDP) mostly by repressing florigen-encoding genes (Ishikawa et al., 2005). Like that in other SDPs, in our previous study we also noticed the repressive effects of NB on quinoa flowering. To know which PEBP genes may be involved in NB responses, we analyzed the transcriptome data of quinoa leaf samples collected from SD and NB conditions. As illustrated in Figure 10 , FT clade genes were more active than TFL and MFT clade genes. A total of 9 PEBP genes were detected, of which 8 were FTL clade genes ( Figure 10 ; Supplementary File 2 ). Obviously, CqFTL3, CqFTL5, CqFTL7, CqFTL8 and CqFTL9 were expressed at higher levels under SD, whereas were rapidly down-regulated by NB treatment ( Figure 10 ). The opposed expression patterns of the 5 FTL genes under SD and NB indicate their possible positive effects in quinoa flowering. In addition, we noticed the differing responses to NB between orthologous gene pairs. For example, in CqFTL4/CqFTL5, CqFTL5 was sensitive to NB, whereas CqFTL4 was only slightly repressed by NB (less than 2 folds). In CqFTL6/CqFTL9, CqFTL9 was largely repressed by NB, while CqFTL6 were down-regulated regardless of light conditions.

As quinoa belongs to SDPs, SD induces whereas long day (LD) represses quinoa flowering. To evaluate the relationships between PEBPs and day lengths, we performed a comprehensive analysis on PEBP expressions over 24 h under SD and LD. The quinoa plants were grown in growth chamber supplied with SD or LD. Before bolting stage, the top-fully-expanded leaves of two-week-old seedlings were collected for RNA-seq. We detected the expressions of 9 FTL clade genes, but only 2 TFL and 1 MFT clade genes were detected in leaves at this stage ( Figures 11A-L ; Supplementary File 2 ). The 24 h expression profile showed that, the expressions of several quinoa PEBP genes were diurnally rhythmic and differed under different day lengths ( Figures 11A, E, F, H, L ). Clearly, expression patterns between the orthologous gene pairs CqFTL1/CqFTL2 ( Figures 11A, B ), CqFTL4/CqFTL5 ( Figures 11D, E ), CqFTL6/CqFTL9 ( Figures 11F, I ) were different. Expressions of CqFTL1, CqFTL4, CqFTL5 and CqFTL8 were considerably higher under SD than that under LD ( Figures 11A, D, E, H ). By contrast, CqFTL6 was induced by LD rather than SD ( Figure 11F ). Expression levels of CqFTL7 and CqFTL9 were slightly higher under SD than that under LD ( Figures 11G, I ). Expression levels of CqFTL3 were comparable between under SD and under LD ( Figure 11C ). CqFTL1 and CqFTL8 expressions were peaked at dawn time, and CqFTL5 and CqFTL9 were abundantly expressed at 2 h after dawn break ( Figures 11A, E, H, I ). CqFTL2 was highly expressed throughout the day under both SD and LD ( Figure 11B ). The expressions of CqTFL1 and CqTFL2 were comparable between under SD and LD ( Figures 11J, K ). CqMFT2 was induced by SD and peaked at the end of day ( Figures 11L ). Together, these results indicated that CqFTL1, CqFTL4, CqFTL5 and CqFTL8 may act as SD-type flowering inducer, whereas CqFTL6 may act as LD-type flowering inducer. Recently, Patiranage et al. also investigated the diurnal expression patterns of FTL genes, but only limited to 6 FTL genes of quinoa at bolting stage (Patiranage et al., 2021). As florigen encoding FT-like genes are highly induced in leaves before floral transition, we chose to investigate the diurnal expressions at vegetative stage before bolting. The classical florigen genes are likely to diurnally expressed and are sensitive to night break (Ishikawa et al., 2005; Meng et al., 2011). By combining the responses to night break and photoperiodic expressions of quinoa FTL genes, we speculated that CqFTL5 and CqFTL8 are the major florigen-encoding genes in quinoa.

As floral integrators, FT genes usually are regulated by a complex network consisting of a lot of transcription factors (TFs) (Tsuji et al., 2013). To know which TFs may lay in the signal cascades upstream of FT genes in quinoa, WGCNA was performed to identify the co-expressed genes. The expression profiles of 32 samples covering 16 time points of SD and LD were used for WGCNA. A total of 6 co-expression modules were obtained. We found only CqFTL5 and CqFTL8 were predicted to occur in those modules. CqFTL5 and CqFTL8 were clustered into the same module-blue module, which contains 934 co-expressed genes. Then those genes were submitted to PlantTFDB v5.0 for TF prediction. A total of 64 TFs, 6.85% of blue module genes, were predicted to be co-expressed with CqFTL5 and CqFTL8. A heatmap was drawn according to their expression profiles ( Figure 12 ). Clearly, those TFs were clustered into two major clusters. One cluster displayed LD-repressive and SD-inducible expression patterns, consistent with that of CqFTL5 and CqFTL8, whereas, the other cluster displayed opposed expression patterns ( Figure 12 ). Thus, those TFs of the two clusters were predicted to be putative promoters and repressors of CqFTL5 and CqFTL8. Notably, compared with the overall TF percentage of blue module genes, remarkable higher rates of family members were found in CO-like (CONSTANS), DBB (Double B-box), GeBP, HSF (Heat Stress Transcription Factors) and bZIP families ( Table 2 ). DBB encodes for a double B-box zinc finger protein, and is regulated by circadian rhythm and participates in light signal transduction during photomorphogenesis. Previous studies indicated that CO-like encodes for a B-box zinc finger protein, and functions as an important mediator of circadian clock to control FT during floral transition (Putterill et al., 1995; Izawa et al., 2002). Thus, like that in other plants, the CO-FT signal pathway is presumably conserved in quinoa flowering. The central roles of MADS-box genes during floral transition have been extensively investigated (Kim et al., 2008; Ryu et al., 2009). We identified that two MADS-box genes were co-expressed with CqFTL5 and CqFTL8. Together, we speculated that these TFs may modulate the transcription levels of CqFTL5 and CqFTL8 in leaves to control quinoa flowering time.