Evolution analysis of FRIZZY PANICLE (FZP) orthologs explored the mutations in DNA coding sequences in the grass family (Poaceae)

FRIZZY PANICLE (FZP), an essential gene that controls spikelet differentiation and development in the grass family (Poaceae), prevents the formation of axillary bud meristems and is closely associated with crop yields. It is unclear whether the FZP gene or its orthologs were selected during the evolutionary process of grass species, which possess diverse spike morphologies. In the present study, we adopted bioinformatics methods for the evolutionary analysis of FZP orthologs in species of the grass family. Thirty-five orthologs with protein sequences identical to that of the FZP gene were identified from 29 grass species. Analysis of conserved domains revealed that the AP2/ERF domains were highly conserved with almost no amino acid mutations. However, species of the tribe Triticeae, genus Oryza, and C4 plants exhibited more significant amino acid mutations in the acidic C-terminus region. Results of the phylogenetic analysis showed that the 29 grass species could be classified into three groups, namely, Triticeae, Oryza, and C4 plants. Within the Triticeae group, the FZP genes originating from the same genome were classified into the same sub-group. When selection pressure analysis was performed, significant positive selection sites were detected in species of the Triticeae and Oryza groups. Our results show that the FZP gene was selected during the grass family’s evolutionary process, and functional divergence may have already occurred among the various species. Therefore, researchers investigating the FZP gene’s functions should take note of the possible presence of various roles in other grass species.


INTRODUCTION
The grass family (Poaceae) consists of many plant species that are essential to mankind, including cereal crops (the main source of grain foods), forage grasses (a valuable source of livestock feed), sugarcane (a key industrial raw material), bamboos, and turf grasses. Grasses produce tillers as well as florets on spikelets, and variation in the number and arrangement of both tillers and spikelets contributes to the great diversity of grass inflorescence structures. In cereal crops, the grains harvested from spikes serve as the main source of grain foods for humans. During the gradual domestication of wild cereal species to cultivars, dramatic changes have occurred in spike morphology (Lee et al., 2020). However, it is unclear whether genes that control the characteristics of spike morphology were selected during the evolutionary process.
FRIZZY PANICLE (FZP) is a key gene that controls spikelet differentiation in rice (Oryza sativa). The FZP gene is closely associated with rice yields as it prevents the formation of axillary bud meristems and permits the formation of floral meristems in rice spikelets. It encodes an AP2/ERF transcription factor and is the rice ortholog of the maize branched silkless1 (BD1) gene (Bai et al., 2017;Chuck, 2002;Komatsu et al., 2003). Komatsu et al. (2003) analyzed the sequences of different FZP mutants and found that amino acid mutations produced severe phenotypes in the conserved AP2/ERF domain and the acidic C-terminus region. This indicates that both the AP2/ERF domain and conserved C-terminus region are necessary for proper FZP function. Maize bd1 mutants are characterized by indeterminate branches in place of female spikelets and a series of lateral spikelets in tassels. Besides, bd1 mutants display similar tassel phenotypes but differ in severity in the ear. Specifically, weaker phenotypic variations initiate fewer branches in the ear and produce a few fertile florets, whereas the most severe phenotypic variations result in the conversion of most initial spikelets to indeterminate branches. In Brachypodium distachyon, the inflorescence is an unbranched spike with a limited number of lateral spikelets and terminal spikelets. The more spikelets1 (mos1) mutant of B. distachyon has more axillary meristems that can develop into branches with higher-order spikelets produced from inflorescence meristems. It is believed that MOS1 is an ortholog of bd1 and FZP (Derbyshire & Byrne, 2013). The inflorescences or spikes of barley (Hordeum vulgare) and common wheat (Triticum aestivum) are unbranched (Bonnett, 1935;Bonnett, 1936) and usually produce one spikelet per rachis node (Dobrovolskaya et al., 2015). The inflorescence meristems of barley and wheat also differ in the formation of spikelets (Derbyshire & Byrne, 2013). Dobrovolskaya et al. (2015) analyzed the structures and functions of three wheat FZP orthologs (WFZP) and found that coding mutations of WFZP-D caused the supernumerary spikelets (SS) phenotype. The most severe phenotypic effect was produced when WFZP-D mutations were combined with a frameshift mutation in WFZP-A.
The results described above demonstrate that FZP and its orthologs play a crucial role in the spike development process of grasses. Mutations in proteins encoded by the gene cause variations in spike shape, which indicates the highly conservative nature of the gene. Grasses possess extremely diverse spike morphologies, and the high grain yield of modern cereal cultivars is the result of long-term domestication of wild grass species. However, it is unclear whether FZP orthologs were selected during the evolutionary process of grasses. In the present study, we identified the FZP gene's orthologs in selected grass species and performed evolutionary analysis using bioinformatics methods to determine the mutations in coding DNA sequences (CDSs) during the evolution process of grasses. Results of this study may serve as a theoretical basis for understanding the functional divergence of the gene among different grass species.

Analysis of conserved motifs and multiple sequence alignment of the FZP gene in Poaceae species
Using the online MAFFT (Multiple Alignment using Fast Fourier Transform) multiple sequence alignment program (https://www.ebi.ac.uk/Tools/msa/mafft/), multiple sequence alignment was performed for the protein sequences and CDSs of the Poaceae species. Alignment results obtained from MAFFT were corrected using MACSE (for CDS) (Ranwez et al., 2011), a built-in alignment software equipped with the PhyloSuite v1.2.1 software (Zhang et al., 2020), and used in subsequent analysis. Use Jalview v2.11.1.4 software (Waterhouse et al., 2009) to beautify the multiple sequence alignment files.

Phylogenetic analysis of the FZP gene in Poaceae species
The Models function in the MEGA X (Molecular Evolutionary Genetics Analysis X) software (Kumar et al., 2018) was used to search for the optimum model for the generation of a maximum likelihood (ML) tree, used full-length CDS alignment sequence of FZP orthologs gene (El-Sappah et al., 2021a). Subsequently, the Phylogeny function was used for the construction of the ML tree and improvement of its visual appearance.
Use Gblocks 0.91b (Talavera & Castresana, 2007) to select the sequence conserved region of the above full-length CDS alignment sequence, MrModeltest 2.3 (Nylander, 2004) selects the optimized nucleotide substitution model and calculates the relevant parameters, Mrbayes 3.04b (Ronquist et al., 2012) reconstructs Bayesian tree of FZP genes. In the process of rebuilding the Bayesian tree, 4 Markov chains are set up, starting with a random tree, running a total of 2 000 000 generations, sampling once every 100 generations and disposing of 25% of the aging samples, and constructing them based on the remaining samples consensus tree, and compute the posterior probability.

Selection pressure analysis of the FZP gene in Poaceae species
We used codeml (Yang, 2007), as implemented in a preset model of EasyCodeML v1.2 software (Gao et al., 2019), to test for natural selection in each of the three sub-groups with the Branch-site model using the CDS alignment results and constructed ML tree from the processes described in the previous two sections. This model allows the ratio of nonsynonymous to synonymous substitution rates (ω) to vary among sites (between sites and branches) to test for site-specific evolution in the context of phylogeny. For each gene region, we evaluated which of codon substitution models (Model A vs. Model A null, Model A: fix_omega = 0, omega = 2; Model A null: fix_omega = 1, omega = 1) was favored using likelihood ratio tests with posterior probability (p > 0.95). In addition, we use the custom model of EasyCodeML v1.2 software to set the model parameters and use different ω values to test the robustness of the analyses result.
Predict the impact of amino acid substitutions or indels on the biological functions of FZP gene protein PROVEAN (Protein Variation Effect Analyzer) is a software tool that predicts whether an amino acid substitution or indel impacts the biological function of the protein. PROVEAN is useful for filtering sequence variants to identify nonsynonymous or indel variants that are predicted to be functionally important. We utilize the PROVEAN online server (http://provean.jcvi.org/seq_submit.php), which predicts whether an amino acid indel or substitution of FZP gene (Os07g0669500) affects the biological function. Software operation reference NCBI nr (September 2012) database.

Protein 3D structure prediction
We use the I-TASSER online server (https://zhanggroup.org/I-TASSER/) to predict the 3D structure of the target gene. Then we did use PyMOL v2.5.0 software (Schrodinger, 2015) to visualize according to the structure file of the prediction result (El-Sappah et al., 2021b).

Analysis of conserved domains of the FZP gene in Poaceae species
A highly conserved AP2/ERF domain is present in all 35 FZP orthologs ( Fig. 1A and Fig. S1). The domain comprised 59 amino acids, with mutations only occurring at positions 13 (QBI22216.1) and 39 (QBI22219.1) (Fig. 1A), which belong to Triticum monococcum and Triticum turgidum, respectively. This indicates that the domain was a core domain that had remained highly conserved throughout the evolutionary process. The amino acids of   8, 10, 16, 18, 25, 27 (asterisks) indicate that confer specific GCC box binding. In a specific Poaceae species category, the acidic domain of C-terminus was also relatively conserved during evolution. Visible differences exist in the acidic domain of C-terminus of these three categories of Poaceae, which serve as the main domain for distinguishing between them ( Fig. 1B and Fig. S1). For instance, position 3 was occupied primarily by Ala (

Phylogenetic analysis of the FZP gene in Poaceae species
Using MEGA X, the optimum nucleotide replacement model T92+G was identified to generate the ML tree for the full-length CDS sequence of FZP gene in Poaceae species (Fig. 2) (Data S3). Results indicate that the grass species could be classified into groups I, II, and III. The bootstrap values of the three groups were 98%, 100%, and 66%, respectively, which were >95% and indicative of high reliability. Group I consisted of species of the  (Milano et al., 2018), and P. miliaceum (Son & Sugiyama, 1992)) and two species of the genus Setaria (S. italica (Heimann et al., 2013) and S. viridis (Brutnell et al., 2010)). The group also consisted of two genes from Z. mays (BD1 and BD1B), one gene from S. bicolor (Paterson et al., 2009;Tazoe et al., 2016), and two genes from Saccharum spontaneum (Zhang et al., 2018).  (Studer et al., 2016). This shows that the FZP orthologs can effectively distinguish between C3 and C4 grasses that diverged earlier, but not between C3 grasses whose divergence occurred more recently. In other words, mutations had already occurred in the amino acid sequence of the gene during the early divergence of C3 and C4 grasses. Although amino acid mutations in the AP2/ERF domain have been rare, distinct mutations have occurred in the C-terminus acidic domain, as shown in Fig. 1. After multiple alignments, redundant sequences were cut to obtain a CDS sequence conserved region with a length of 438 bp. According to the AIC standard of MrModeltest software, the optimum nucleotide replacement model for BI tree is F81+G, and the rate of variation between sites rates = gamma (Lset nst = 1 rates = gamma; Prset statefreqpr = Dirichlet (1, 1, 1, 1)), The FZP gene Bayesian phylogenetic tree reconstructed by Bayesian inference (BI) tree is shown in Fig. S2 and supplementary data 4. The topology is consistent with the ML tree.

Detection of selection pressure of the FZP gene in Poaceae species
The three branches; Triticeae, Oryza, and C4 plants, were set separately as the foreground branch (Fig. 2) and all other branches were set as the background branches for the analysis of selection pressure of the FZP gene in grasses using the branch-site model in EasyCodeML v1.2. When Triticeae was marked as the foreground branch, the likelihood ratio test (LRT) was highly significant (P = 0.000017894). Two positive selection sites were detected, with a significant positive selection site (61 I, posterior probability >95%) ( Table 2 and Fig. S1) located near the upstream region of the AP2/ERF domain. When the Oryza branch was marked as the foreground branch, the LRT was significant (P = 0.030827376). Three positive selection sites were detected, including a significant positive selection site (484 G, posterior probability >95%) (Table 2) located at the C-terminus. When the C4 plants branch was marked as the foreground branch, the LRT was not significant (P = 0.127622034). Six positive selection sites were detected, including 401 E sites within the acidic domain of C-terminus and one significant positive selection site (328 S, posterior probability >99%) ( Table 2 and Fig. S1) near the upstream region of the acidic domain of C-terminus.

Predict the impact of amino acid substitutions or indels on the biological functions of FZP gene protein
In order to understand whether an amino acid substitution or indel has an impact on the biological function of protein. We used PROVEAN online server (http: //provean.jcvi.org/seq_submit.php) to predict amino acid substitutions or indels of FZP gene (Os07g0669500) on the biological functions. The results show that the significant positive selection sites (sites with asterisks) detected in 3.4, PROVEAN scores are above −2.5, variants are considered ''neutral'' (Table 3). However, the Y69C (score = −8.925) and L95P (score = −6.242) amino acid substitutions in the AP2/ERF domain of QBI22216.1 and QBI22219.1 gene are considered ''deleterious'' (Table 3). The nine amino acid deletions in the C-terminus of the AAX83534.1 gene in Oryza longistaminata (Y310_H318del, score = −4.189) are considered ''deleterious'' (Table 3).     (Fig. S1). c Default threshold is −2.5, variants with a score equal to or below −2.5 are considered ''deleterious'', if above −2.5 are considered ''neutral''.

Protein 3D structure prediction
We did use I-TASSER online server (https://zhanggroup.org/I-TASSER/), using default parameters, predicted the protein 3D structure of rice ZFP gene (Os07g0669500), wheat WFZP-A gene (TraesCS2A02G116900.1), maize BD1 gene (GRMZM2G307119), T. monococcum QBI22216.1 gene and T. turgidum QBI22219.1 gene to understand whether amino acid substitutions or indels in the conserved amino acid domain of the FZP gene affect the 3D structure of the protein. The results (Fig. 3 and Table S1) show that the 3D structure of these five genes is similar to the Arabidopsis AtERF96 gene 3D structure (Protein Data Bank [PDB] ID: 5WX9) (Chen et al., 2020) (Fig. 3A). There are seven amino acid residues (amino acids with labels in the 3D structure at Fig. 3) that specifically bind to the DNA GCC box at AP2/ERF domain (red α-helix, β-sheet and random coil in 3D structure at Figs. 3A-3F). The 3D structure of Arabidopsis AtERF96 gene does not have an acidic domain at the C-terminus (Fig. 3A), and gramineous plants all have an acidic domain (the orange α-helix, β-sheet and random coil in the 3D structure at Figs. 3B-3F) at C-terminus (green coil/ α-helix in 3D structure at Figs. 3B-3F). The C-terminus acidic domain of gramineous and the AP2/ERF domain are close to each other, especially close to Figure 3 The 3D structure of Arabidopsis AtERF96 gene and FZP orthologous gene protein. The red and orange area in 3D structure indicate N-terminus AP2/ERF domain and C-terminus acidic domain of FZP genes, respectively. The C-terminus, α-helix, β-sheet and DNA were point out with black, purple, blue and orange arrows, respectively. The light blue amino acids with labels in the 3D structure indicate that they specifically bind to the DNA GCC box. The red amino acids with labels (P96 and C70) indicate that an amino acid substitution in the AP2/ERF domain are considered ''deleterious'' (PROVEAN predicted results). (A) Crystal structure of AtERF96 with GCC-box. (B-F) The 3D structure of Triticum monococcum QBI22216.1 gene, Zea mays BD1 gene (GRMZM2G307119), Triticum turgidum QBI22219.1 gene, Triticum aestivum WFZP-A gene (TraesCS2A02G116900.1) and Oryza Sativa L. spp. japonica FZP gene (Os07g0669500). The 3D structure of B-F graphics is visualized using model 1 by PyMOL software.
Full-size DOI: 10.7717/peerj.12880/ fig-3 the region of amino acid residues that specifically bind to the DNA GCC box (except for the QBI22216.1 gene). The substitution of P96 (L183P in Fig. S1) amino acid residue in the AP2/ERF domain of T. turgidum QBI22219.1 gene did not cause major changes in the 3D structure of the protein (Fig. 3D). It is worth noting that the substitution of C70 (Y157C in Fig. S1) amino acid residues in the AP2/ERF domain of T. monococcum QBI22216.1 gene causes a huge change in the 3D structure of the protein. The AP2/ERF domain cannot form an α-helix, and the C-terminus acidic domain was far away from AP2/ERF domain (Fig. 3B).

DISCUSSION
In this study, a total of 34 orthologs of the FZP gene were identified in 26 grass species. No complete or highly identity FZP orthologous genes have been searched in Triticum urartu, Eragrostis tef, and Oryza brachyantha, which are also gramineous plants. There are two possible reasons for these results: (1) The genome assemblies of T. urartu, E. tef, and O. brachyantha may be incomplete; and (2) an additional level of genetic regulation is required for inflorescence formation in grasses (Komatsu et al., 2003). Specifically, maize (Z. mays) possessed two FZP orthologs (BD1 and BD1B), located on chromosome 2. Mutations in the BD1 gene cause the formation of indeterminate lateral tillers in the ear and tassel, whereas the functions of the BD1B gene have not yet been reported. This indicates that the BD1B gene maybe a duplicate of BD1 generated during the maize genome duplication. The common wheat (T. aestivum) possessed three FZP orthologs (WFZP) that separately belonged to the A, B, and D genomes. The results of the evolutionary analysis revealed that the WFZP-A gene of common wheat was located on the same branch as the FZP ortholog of T. monococcum, the WFZP-B gene was located on the same branch as the FZP ortholog of T. dicoccoides, and the WFZP-D gene was located on the same branch as the FZP ortholog of A. tauschii. These results demonstrate that the evolutionary process of the FZP gene in common wheat was consistent with the evolutionary theories for common wheat. Modern theories have generally asserted that common wheat was domesticated in the Fertile Crescent approximately 10,000 years ago (Lev-Yadun, Gopher & Abbo, 2000) with its hexaploid genome formed through the hybridization of tetraploid wheat with the wild diploid grass A. tauschii (Abbo et al., 2014). Except for T. turgidum, T. dicoccoides, and Saccharum spontaneum, which possessed two FZP orthologs, all other species had a single FZP ortholog that only existed in plants, indicating that the evolutionary process of the genes occurred more recently. Our analysis of the conserved motifs of FZP orthologs in the grass species revealed that amino acid mutations were almost absent in the AP2/ERF domain, which indicates the highly conserved nature of the domain during the evolutionary process. It is known that the occurrence of nonsynonymous substitutions in the AP2/ERF domain leads to changes in gene function. However, the rich evolutionary information has been retained through significant mutations in the acidic C-terminus region of the gene in species of the tribe Triticeae, genus Oryza, and C4 plants, which can be used for phylogenetic analyses. Interestingly, except for species belonging to the genera Oryza, B. distachyon, Eragrostis curvula, and Panicum miliaceum, the FZP genes of all other grass species were located on chromosome 2. Previous research has indicated that chromosome 2 of foxtail millet (Setariaitalica) was formed by combining chromosomes 7 and 9 of rice. Besides, the researchers observed the occurrence of fusion events in the chromosomes of Sorghum (Sorghum bicolor), which led to the deduction that the chromosome fusion events most likely occurred before the divergence of sorghum and foxtail millet. These chromosomal recombination events form a key basis for species divergence and genetic variation in species (Bennetzen et al., 2012;Zhang et al., 2012). Besides maize and polyploid grasses, all other grasses only possessed a single FZP ortholog. Researchers have also found that chromosome 2 of foxtail millet was formed through the fusion of chromosomes 7 and 9 of rice and that similar fusion events had occurred in the chromosomes of Sorghum, which led to the deduction that the chromosome fusion events most likely occurred before the divergence of sorghum and foxtail millet. Researchers had also discovered a specific chromosome fusion event after the divergence of foxtail millet. These chromosomal recombination events form a key basis for species dissimilarity and genetic variation in species (Bennetzen et al., 2012;Zhang et al., 2012). During selection pressure analysis of FZP orthologs in Poaceae using the branch-site model, significant positive selection sites were distinguished in the Triticeae and Oryza branches (Table 2 and Fig. S1), which demonstrates that the FZP gene in the tribe Triticeae and genus Oryza was decidedly chosen to various extents during the divergence process and may have gone through functional divergence. This may clarify the substantial differences in spike morphology between the two groups of grass species, and specialists should observe such divergence when investigating the functions of the FZP gene.
To understand whether positive selection sites and other amino acid substitution or indel has an impact on the biological function of the protein. We use the PROVEAN online server (http://provean.jcvi.org/seq_submit.php), which predicts whether amino acid substitution or indel of FZP gene (Os07g0669500) has an impact on the biological function. The results show that the significant positive selection sites (PROVEAN scores are above −2.5) are considered ''neutral'' (Table 3 and Fig. S1). However, the Y69C (score = −8.925) and L95P (score = −6.242) amino acid substitutions in the AP2/ERF domain of QBI22216.1 and QBI22219.1 gene are considered ''deleterious'' (Table 3 and Fig. 1A). These results indicate that although positive selection sites have been detected in the Triticeae dumort and Oryza branches, the substitution of these sites may not affect the gene's function. AP2/ERF domain is highly conserved in gene evolution, and the substitution of amino acids may affect gene function. Previous studies have also proved this point. Komatsu et al. (2003) studied rice mutants (fzp-2 and fzp-7 ), they found that one of the 7 amino acids (R25L in Fig. 1A and R25W in Fig. 1A) that specifically bind to the GCC box (Allen et al., 1998;Chen et al., 2020;Hao, Ohme-Takagi & Sarai, 1998) in the AP2/ERF structure is replaced, which can lead to a severe frizzy panicle (fzp) phenotype. Even if the conserved amino acids that non-specifically bind to the GCC box in the AP2/ERF domain are replaced (T30M in Fig. 1A), it can result in a weaker fzp phenotype. It indicates that FZP gene may regulate its target gene through GCC box-mediated gene expression. Certain nucleotide variations at the end of the AP2/ERF domain may produce premature stop codons (Q52Stop in Fig. 1A), leading to a severe fzp phenotype. Dobrovolskaya et al. (2015) detected a mutation in the GCC-box binding site (E16K in Fig. 1A) and a premature stop codon (Q35Stop in Fig. 1A) in the AP2/ERF functional domain of FZP-D resulted in the supernumerary spikelets (SS) phenotype at bread wheat mutant lines (wfzp-d.1 and wfzp-d.2). Dobrovolskaya et al. (2015) also found similar mutants (Bd8202 and Bd8972) in Brachypodium distachyon (Bd8202 and Bd8972), showed an SS phenotype but produced viable seeds. Two transitions were detected: G232A in Bd8972 and G163A in Bd8202. These mutations lead to amino acid substitutions (G11S and A34T in Fig. 1A) in the conserved region of the AP2/ERF domain. Chuck (2002) analysis of mutant alleles (bd1-N and bd1-3) revealed a transition mutation (W10Stop in Fig. 1A) within the AP2/ERF domain that introduces a premature stop codon in bd1-N, and deletion of the 3 end (C-terminus) of the gene that is predicted to remove the last ten amino acids in bd1-3. The nine amino acid deletions (Y310_H318del in Fig. S1) in the C-terminus of the AAX83534.1 gene in Oryza longistaminata are considered ''deleterious'' (Table 3 and Fig. S1) in this research. These results indicate that the C-terminus is also very important for the biological function of the FZP gene. Komatsu et al. (2003) also believe that the C-terminus strongly indicate that this conserved region must also be necessary for proper FZP function in addition to the AP2/ERF domain through the study of FZP gene mutants.
Through previous studies on FZP gene mutants and our analysis of FZP orthologous genes, we believe that although the evolution of FZP genes in Gramineae species is relatively conservative, especially AP2/ERF domain, it is not ruled out that its function will change in some gramineous species. We also observed the C70 mutation in the AP2/ERF domain of the QBI22216.1 gene (T. monococcum) ( Fig. 3 and Table 3) in the prediction results of the protein 3D structure of the 5 FZP orthologous genes. This mutation occurred in the seven amino acids region that specifically binds to the GCC box lead to the 3D structure of the whole protein has undergone major changes, the AP2/ERF domain cannot form an α helix, and the C-terminus acidic domain is far away from the AP2/ERF domain, which may directly affect the biological function of the gene.

CONCLUSION
FZP, an essential gene that controls the development and the spikelet differentiation in the grass family (Poaceae), prevents the formation of axillary bud meristems and is closely associated with crop yields. It is unclear whether the FZP gene or its orthologs were selected during the evolutionary process of grass species, which possess diverse spike morphologies. Hence, in the present study, we identified the orthologs of the FZP gene in selected grass species and performed evolutionary analysis using bioinformatics methods to determine the mutations in coding DNA sequences during the evolution process of Poaceae. The phylogenetic analysis showed that the 29 Poaceae species could be classified into three groups, namely, Triticeae, Oryza, and C4 plants. Within the Triticeae group, the FZP genes originating from the same genome were classified into the same sub-group. Our analysis revealed that amino acid mutations were almost absent in the AP2/ERF domain, which indicates the highly conserved nature of the domain during the evolutionary process. We believe that our study makes a significant contribution to the literature because the results of this study may serve as a theoretical basis for understanding the functional divergence of the gene among different grass species.