Comprehensive study of rice YABBY gene family: evolution, expression and interacting proteins analysis

Sequence alignment and phylogenetic analysis

YABBY sequences of rice, Arabidopsis, tomato (Solanum lycopersicum) and maize (Zea mays) were obtained from the Rice Annotation Project Database (https://rapdb.dna.affrc.go.jp/) (Sakai et al., 2013), the TAIR10 database (https://www.arabidopsis.org/) (Lamesch et al., 2012), the SOL Genomics Network (https://solgenomics.net/) (Fernandez-Pozo et al., 2015), and the MaizeGDB database (https://www.maizegdb.org/) (Portwood et al., 2019), respectively. Accession numbers of YABBYs were listed in Table S1. Multiple sequence alignment was performed using ClustalX (Larkin et al., 2007) with default values. Then neighbor-joining (NJ) phylogenetic tree and maximum-likelihood (ML) tree were constructed using MEGA7.0 (Kumar, Stecher & Tamura, 2016). For NJ tree construction, the parameters were set to the poisson model, 1,000 bootstrap replicates, and the bootstrap value more than 50% was listed at the branches, and those for ML tree are Jones-Taylor-Thornton (JTT) Model.

Gene structure, motifs and conserved domains prediction

Exon-intron structures of OsYABBYs were confirmed from the Rice Annotation Project Database (Sakai et al., 2013), and the conserved domains of OsYABBYs were identified by the uniprot database (https://www.uniprot.org/) (MacDougall et al., 2020). MEME website (https://meme-suite.org/meme/tools/meme) (Nystrom & McKay, 2021) was used to predict the motifs of OsYABBYs with the following parameter settings: the maximum number of motifs, 10; the minimum width and maximum width of motifs, 5 and 20; and default parameters. Eventually, gene structures, motifs, and conserved domains were visualized by IBS1.0.3 (Liu et al., 2015).

Plant growth conditions, RNA extraction and qRT-PCR

Rice variety Zhonghua 11 (ZH11, Oryza sativa ssp Japonica/geng) was cultivated in plant growth room. Seven-day-old seedlings were used for sampling. Roots, stems, and leaves were taken followed by a quick grind in liquid nitrogen. Similarly, reproductive tissues like spikelet (40–50 mm), embryo (7–10 days after pollination), stamen (1 day before flowering), pistil (1 day before flowering), palea/lemma (1 day before flowering) at reproductive stage were sampled and immediately stored at −80 °C. Total RNA from different tissues was isolated by TransZol (Transgen ET101-01; TransGen Biotech Co., Ltd., Beijing, China) reagent, and detected by NANODROP 1000 spectrophotometer to determine the extraction quality and concentration, then reversed transcribed into cDNA with M-MLV reverse transcriptase (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s instructions. Specific primers for all OsYABBYs were designed using online programs (https://sg.idtdna.com/scitools/Applications/RealTimePCR/) (Table S2). The RT-qPCR was performed on the QuantStudio 6 Flex real-time PCR instrument (Applied Biosystems, Foster City, CA, USA) using TB Green® Premix Ex Taq™ II reagent (Takara RR820A; Takara, Kusatsu, Japan). OsActin1 was used as an internal reference to standardize gene expression levels, and each cDNA was subjected to three biological replicates. Relative expression values were calculated using the 2^−ΔΔCt method (Livak & Schmittgen, 2001) and the heatmap of gene expression levels was drawn using TBtools (Toolbox for Biologists) v1.09876 software (Chen et al., 2020a).

PPIs network prediction

OsYABBY interaction proteins were predicted by STRING website (https://string-db.org) (Szklarczyk et al., 2021). Predicted interacting protein score was set to a minimum of 0.4, predicted number of direct interacting proteins was set to no more than 50, and the number of secondary interacting proteins was set to zero. Interaction information was visualized by Cytoscape 3.7.1 (Shannon et al., 2003). Nodes represented proteins, the connections between nodes were represented by edges, which hinted the interactions between these biological molecules (Bader & Hogue, 2003).

Molecular docking simulation

Based on PPI prediction of interacting proteins, RGAP (Rice Genome Annotation Project) database (http://rice.uga.edu/index.shtml) (Kawahara et al., 2013) was used to obtain protein sequences, and AlphaFold website (https://www.alphafold.ebi.ac.uk/) (Varadi et al., 2022) was used to get PDB files containing 3D structure information. ZDOCK SERVER (http://zdock.umassmed.edu/) (Chen, Li & Weng, 2003) was applied to dock prediction between two proteins by inputting the PDB files. The spatial structures were visualized by PyMOL (Yuan et al., 2016), and the closest distance between possible interacting proteins was measured by the measurement plug-in.

Yeast two-hybrid

The full length of OsYABBY cDNA was amplified and cloned into the bait vector (pGBKT7), and the CDS of OsWOX3A and OsWUS were obtained and cloned into prey vector (pGADT7). The fused pGBKT7 and pGADT7 vectors were transformed into yeast cells (AH109). Transformed cell growth status on the SD/-Leu/-Trp/-His/-Ade medium was used to determine whether the proteins interact with each other. The primers used to construct pGBKT7 and pGADT7 vectors were listed in Table S2.

Luciferase complementation imaging assay (LCI)

OsYABBYs were fused with the amino-terminal of luciferase to construct the N-luc vector, and OsYABBY-interacted protein coding genes were fused with the carboxyl-terminal of luciferase to construct the C-luc vector. The two fusion proteins were simultaneously expressed in tobacco (Nicotiana benthamiana) through agrobacterium-mediated transformation. If the two proteins interacted with each other, the two parts of luciferase would be recombined into whole luciferase, which could oxidate luciferin and generate bioluminescence detected by Tanon5200 (Sun, Zheng & Zhu, 2017).

Phylogenetic analysis of YABBY family proteins

To reveal the phylogenetic relationship of YABBY proteins, amino acid sequences of Zinc finger and YABBY domains of 37 YABBYs, including eight from rice, six from Arabidopsis, nine from tomato and 14 from maize were aligned and constructed two phylogenetic trees. In NJ (neighbor-joining) tree (Fig. 1), all these YABBYs from four species were clearly divided into five groups: FIL/YAB3, CRC, YAB5, INO and YAB2, based on the similarity of amino acid sequences. OsYABBYs existed in four groups excluding YAB5 group, and ZmYABBYs was also not in YAB5 group, which was consistent with the phylogenetic analysis of other monocot YABBYs (Romanova et al., 2021; Jie et al., 2022) and suggested that the members in YAB5 group might had unique function of dicotyledons. In both CRC and INO groups, dicotyledonous and monocotyledonous YABBYs clustered into a subgroup, respectively, indicating that the function of these YABBYs might have obvious differentiation. However, the differentiation of dicotyledonous and monocotyledonous YABBYs in the FIL/YAB3 and YAB2 groups was not as significant as that in CRC and INO groups. In the ML tree (Fig. S1), similar results to the NJ tree was shown, indicating that the classification of YAB proteins of the NJ tree was reliable.

Gene structure, motif and conserved domain identification

Gene structure can provide effective information for understanding gene evolution and function. Structural analysis revealed that FIL/YAB3 and CRC group genes (OsYABBY3/4/5, DL) contained seven exons, INO and YAB2 group genes (OsYABBY7, OsYABBY1/2/6) had six exons (Fig. 2A). Members of the same group were more similar to each other than to members of other groups, indicating that each group of OsYABBYs might have specific functions.

The function of a protein is determined by its domains, so the motifs and domains of OsYABBYs were analyzed. Six motifs were predicted among all YABBY proteins (Fig. 2B, Fig. S2). OsYABBYs in the same group had similar motifs. For example, YABBY3/4/5 all had motif5, but other YABBYs did not. It might be a motif unique to FIL/YAB3 group (Fig. 2B). Motif1 and motif4 constituted the YABBY domain of OsYABBYs, and motif2 and motif3 constituted the zinc finger domain (Fig. 2C). Both the zinc finger domain and the YABBY domain were considered to be responsible for binding DNA, in which YABBY domain was the main binding site and was also required for YABBY’s nuclear localization (Gross, Broholm & Becker, 2018). However, not every OsYABBY had motif4. Absence of motif4 in some OsYABBYs suggested that the sequences of YABBY domain might be different. Changes in sequence often lead to changes in function. In the site directed mutagenesis of AtCRC, it was found that mutations in C2C2 zinc finger domainat affected the DNA-binding ability, and mutations in YABBY domain destroyed 3D structure and DNA-binding ability (Gross, Broholm & Becker, 2018). Within the amino acid alignment of the domains between OsYABBYs and AtCRC (Fig. S3), it could be seen that a few sites mutated in YABBY domain, suggesting that OsYABBYs might have different DNA-binding ability.

Expression level of OsYABBYs in various tissues

Plant-specific TRs, such as OsWOX11 (Zhao et al., 2009), OsLBD6 (OsIG1) (Zhang et al., 2015), OsIDD10 (Xuan et al., 2013), play vital roles during organogenesis and exhibit specific expression profiles. The expression data of OsYABBYs in various tissues (Figs. 3A and 3B) provided by the RiceXpro website (http://ricexpro.dna.affrc.go.jp/) were used to construct a heat map (Fig. 3C, Table S3). Except for the constitutive expression of DL, most OsYABBYs were predominantly expressed in reproductive tissues (Fig. 3C).

To verify the reliability of the microarray data of the RiceXpro website, transcriptional level of OsYABBYs in vegetative and reproductive organs was detected by RT-qPCR (Fig. 4). All OsYABBYs were expressed at a higher level in reproductive organs (spikelet, embryo, stamen, palea/lemma and pistil) than those in vegetative organs (root, stem and leaf). OsYABBY1/2/6 (YAB2 group) and OsYABBY7 (INO group) displayed high expression in palea/lemma, suggesting these genes might redundantly regulate palea/lemma development. OsYABBY1 expressed widely in various tissues, which was in accordance with previous report (Dai et al., 2007b) and indicated OsYABBY1 might participate in different developmental processes. FIL/YAB3 group genes (OsYABBY3/4/5) had higher transcript in spikelet and weaker in root, demonstrating that they might co-regulate spikelet growth. Different from the constitutive expression of OsDL in microarray data, RT-qPCR analysis showed tissue-specific expression pattern of OsDL, which was mainly expressed in pistil, palea/lemma and spikelet. These results suggested that different OsYABBYs might be responsible for the development of different reproductive tissues.

Prediction and verification of OsYABBY interacting proteins

In order to further elucidate OsYABBY functions, the interacting proteins of OsYABBYs were investigated through STRING website (https://string-db.org). The results revealed that some OsYABBYs could interact with each other (Fig. S4), and also form complex with other proteins (Fig. S5). Since the complexes that formed between YABBYs had been reported (Sieber et al., 2004; Stahle et al., 2009), the more other proteins that could interact with YABBYs needed to be elucidated. As TRs, YABBYs function in the nucleus, their interaction candidate proteins should also be localized in the nucleus. Based on this standard, 17 non-YABBY proteins were screened, including transcription factors (TFs), chromosome recombination factors, elongator-associated factor, zinc transporter, gibberellin regulatory protein, enzymes and chaperone (Fig. S5A, Table S4). Among these, the number of TFs accounted for 70% of all proteins, and the WOX (WUSCHEL-related homeobox) family proteins accounted for 30% of TFs, indicating that OsYABBYs and OsWOXs might interact with each other (Fig. S5B). It was reported that some OsWOX genes (OsWUS, OsWOX3A, OsWOX9A) were highly expressed in reproductive organs and function in reproductive organ development (Cho et al., 2013; Cheng et al., 2014), which was similar to expression pattern and role of OsYABBYs. Based on interaction prediction and expression profiling analysis of OsWOXs and OsYABBYs, interaction verification of OsYABBYs and OsWOXs was further conducted by molecular 3D docking, Y2H and LCI experiments.

The 3D docking of OsYABBYs and OsWOXs was executed by ZDOCK program, which was a docking program to predict several protein complexes via PSC (Pairwise Shape Complementarity) of input protein structure. The docking results were visualized by PyMOL software (Fig. 5), and the closest distances measured by PyMOL and the closest distances between hydrogen bonds were organized into Table S5. Two groups of proteins YABBY3-WUS and YABBY5-WUS were the closest in space, with distances of 1.0 and 0.9 Å, respectively. The closest distances of YABBY3-WOX3A and YABBY5-WOX3A are 3.2, and 3.1 Å, respectively. The reliability of the STRing website results was confirmed by the distance between proteins, and the information between hydrogen bonds can also predict these protein interactions.

The Y2H suggested that OsYABBYs except OsYABBY7 could interact with OsWOX3A, and OsYABBY3, OsYABBY5 could interact with OsWUS (Fig. 6A). LCI experiment further verified the results of Y2H (Fig. 6B). It could be seen from the Y2H and LCI that OsYABBY5-OsWUS, OsYABBY4-OsWOX3A, OsYABBY6-OsWOX3A had strong interactions, while DL-OsWOX3A had weaker interaction (Fig. 6). The results showed that there was indeed possibility of interaction between OsYABBYs and OsWOXs, which would lay a foundation of interpreting OsYABBYs regulatory mechanisms.