Genome-wide identification and expression profile of YABBY genes in Averrhoa carambola

Data sources

A genome sequence (accession: SAMC124704) and a GFF3 file (accession number: GWHABKE00000000) of Ave. carambola were downloaded from the China National Center for Bioinformation (https://bigd.big.ac.cn/). The transcriptome data of leaves, inflorescences, bracts, and fruits at different developmental stages were sequenced by Biomark Biotechnology Cooperation, Transcriptome alignment and assembly are performed using Hisat (Kim, Langmead & Salzberg, 2015) and Stringtie2 (Kovaka et al., 2019), respectively. The raw data can be found under the following accession numbers: SRR17036676–SRR17036678, SRR15860439–SRR15860441, SAMN21399627–SAMN21399632.

Identification and physicochemical properties of YABBY genes in star fruit

To identify all candidate YABBY genes in star fruit, a local BLASTP search with a threshold e-value of 1e⁻¹⁰ was performed using A. thaliana YABBY protein sequences as query sequences. Identity and coverage areas (>50%) were used as filtering criteria to eliminate inappropriate YABBY genes. In addition, a pfam seed model (PF04690) was obtained from the online database (http://pfam.xfam.org/) and was used for building a hidden Markov model (HMM) file using HMMER3 software with default parameters. The HMM search program was performed to search YABBY genes from .hmm file generated in the previous step. The results of HMM and BLASTP were compared, and repeats were removed manually. Truncated peptides and proteins at the same chromosomal position were eliminated. To further verify the reliability of the selected sequences, the CDD (https://www.ncbi.nlm.nih.gov/) was used for domain analysis to ensure the presence of the YABBY domain in each candidate AcYABBY gene. In addition, the protein isoelectric point and molecular weight of AcYABBY genes (AcYABBYs) were predicted using the ProtParam tool (Gasteiger et al., 1999; https://web.expasy.org/protparam/).

Analysis of conserved motifs, conserved domains, and subcellular location prediction

To identify shared motifs and structural divergences among the proteins encoded by the YABBY genes, YABBY protein sequences were subjected to analysis using online software MEME 5.1.1’s (Bailey et al., 2006; http://meme-suite.org/tools/meme), and the number of motifs was set to 15. NCBI Batch CD-Search was used to quantitatively predict the structural domains and conserved sites of genes, and DNAMAN v9.0 was used to visualize the domain. Subcellular localizations were predicted using the three protein subcellular location prediction tools Softberry (http://www.softberry.com/), PSORT (https://www.genscript.com/wolf-psort.html?src=leftbar), and LocTree3 (https://www.rostlab.org/services/loctree3/).

Multiple sequence alignment and phylogenetic analysis

Multiple sequence alignments of YABBY proteins were performed using the MUSCLE alignment function in MEGA 7.0 software with default settings (Kumar et al., 2018). The aligned sequences were saved with a .fasta extension by choosing Export Alignment from the Data menu for exporting the file. After this, the fasta file was converted to a .phy file using EasyCodeML (Gao et al., 2019) and was submitted to CIPRES (https://www.phylo.org/) for phylogenetic analyses (Miller, Pfeiffer & Schwartz, 2010). A Maximum Likelihood (ML) method phylogenetic tree was constructed through RAxML (Stamatakis, 2014) under a GTRGAMMA substitution model with 1,000 bootstraps. The resulting phylogenetic data was exported as a Newick file to produce a phylogenetic tree in the FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/), and polar tree layout was chosen to visualize a ring phylogenetic tree.

Gene structure analysis and chromosomal localization

The GFF3 file of the identified AcYABBYs was submitted to the online GSDS v2.0 (http://gsds.gao-lab.org/) (Hu et al., 2014) for gene structure analysis, and the output was saved in SVG format, the convertio, an online image manipulation tool was used to convert SVG format to JPEG format. To determine the chromosomal position of AcYABBYs, the sequences of AcYABBYs were mapped to the genome of star fruit, according to the coordinates of each YABBY gene on the genome.

Prediction of secondary structures of AcYABBYs

To investigate the secondary structure of AcYABBYs, the whole AcYABBY protein was submitted and predicted using http://bioinf.cs.ucl.ac.uk/.

Prediction of promoter elements

To identify putative cis-acting elements in the promoter, we used TBtools to obtain 2,000 bp gene sequence upstream of the promoter codon from the star fruit genome sequence (Chen et al., 2020a). Cis-acting elements in the promoter region were analyzed using PlantCARE (Higo et al., 1999). Microsoft Excel 2010 was applied to the screen promoter elements to produce a histogram.

Colinear and selective pressure

To identify the pattern of gene duplication, MCscanX (Wang et al., 2012) was used to analyze YABBY genes in star fruit vs A. thaliana, and star fruit vs grape. Two piars of blast result files from star fruit vs A. thaliana, A. thaliana vs star fruit and star fruit vs grape, grape vs star fruit were merged. MCscanX was then used to examine the merged blast file and the merged gff3 file that met the specifications. The dual_synteny_plotter script program was used to visualize the merged gff3 file, the collinearity file obtained in the previous step, and the ctl file. To assess the selection pressure of genes encoding YABBY proteins, the ratio of Ka/Ks (an indicator of selective pressure) was used to evaluate its evolutionary pressure. MAFFT 7.0 (Nakamura et al., 2018) was applied to align the protein sequence of YABBY genes in star fruit, and then the pal2nal.pl and csplit scripts were used to compare the aligned protein sequence with the CDS sequences of AcYABBYs and to split the aligned file. The split files were then used to form new gene pair files, and the parseFastaIntoAXT.pl script was used to convert the fasta file into an axt file. The KaKs_Calculator software (Wang et al., 2009) was used to calculate Ka/Ks values of the gene pairs.

Expression analysis of AcYABBYs

To obtain more information regarding the roles of YABBY genes in star fruit, RNA sequencing data were used to quantify the expression levels of YABBY genes in different organs and developmental stages from star fruit, and the raw count number of mapped reads for each AcYABBY (reads per kilo base per million mapped reads) was calculated using StringTie (Kovaka et al., 2019). A gene expression profile heatmap was then produced using Tbtools (Chen et al., 2020a).

RNA extraction and Reverse Transcription qPCR (RT-qPCR)

Since the ripening cycle of star fruit is about 2 months, we use the fruits 20, 40, 60 days after pollination (F-DAP20, F-DAP40, F-DAP60) to represent the early, middle and late stages of fruit development. Leaves, inflorescences, buds, fruits from F-DAP20, F-DAP40, F-DAP60 were collected and frozen in liquid nitrogen for 30 min, and were stored at −80 °C until RNA extraction which was performed as quickly as possible (Fig. S1). Total RNA was extracted using the RNAsimple Plant Kit (Tiangen Biochemical Technology Company, Beijing, China) following the manufacturer’s instructions.

RT-qPCR was performed to further confirm the reliability of the expression profile results using all genes (from AcYABBY1 to AcYABBY8). Total RNA of all collected samples was extracted using the TIANGEN DP441 Reagent (TIANGEN, Beijing, China) following the manufacturer’s instructions. RT-qPCR analysis was performed using a Roche detection system (Roche, Switzerland) with SYBR Green assays. Gene-specific primers for the RT-qPCR analysis of eight selected genes and a reference gene are listed in Table S1. The reaction conditions were 4 min at 95 °C and 40 cycles of 10 s at 95 °C and 40 s at 60 °C. The specificity of the amplicon for each primer pair was verified by melting curve analysis. The melting curve conditions were 95 °C for 15 s, 60 °C for 1 min, 95 °C for 30 s, 60 °C for 15 s, and each sample showed only one melting temperature peak. All experiments were performed in three biological replicates, and each replicate was measured thrice. The log2 fold change was calculated using the 2−ΔΔCT method with Yangtao2017018, a homolog of beta-actin (ACTB), as a reference gene to normalize the target gene expression and to correct for variation between samples.

Identification and physicochemical properties of YABBY genes in star fruit

Local BLAST and HMM analyses were performed, and genes such as pseudogenes, premature stop codons genes, or those without a complete YABBY domain were removed. Finally, eight remaining putative functional YABBY genes with high confidence were identified in Ave. carambola. All genes contained the YABBY domain (Fig. S2). In addition, the YABBY genes varied substantially in the length of encoded protein sequences, ranging from 169 to 235 bp (Table S2). The molecular weight of these AcYABBYs ranged from 18,053.34 to 25,197.26 Da. Most of the AcYABBYs exhibited alkaline isoelectric points higher than 8.00, with the highest being 9.56 for AcYABBY7, while two proteins had acidic isoelectric points below 8.00, of which AcYABBY2 had the lowest at 5.22.

Analysis of conserved motifs and subcellular location predictions

To identify motifs in the AcYABBYs, motifs of YABBY proteins in A. thaliana and star fruit were analyzed using the online analysis tool MEME, and 15 motifs were set. Motif 1 and motif 2 encoded the YABBY and zinc finger domains, respectively, which are the most conserved domains with 50 and 48 amino acids, respectively (Fig. S3). All YABBY genes in A. thaliana and star fruit had these two motifs. Motif 3, motif 5, and motif 8 were specifically distributed in the FIL/YAB3 subclade of A. thaliana and star fruit, but motif 6, motif 11 were specifically distributed in the FIL/YAB3 subclade of A. thaliana or star fruit, respectively (Fig. S3). Motif 4 and motif 7 were present in the YAB2 and YAB5 subclades. Motif 9 and motif 13 were specifically distributed in the INO clade. Motif 10 and motif 14 were specifically distributed in the YAB2 clade. The motif 12, motif 15 were specifically distributed in the YAB5, CRC subclade, respectively. Notably, though AcYABBY3, AcYABBY4, and YAB2 belonged to YAB2 clade, only ACYABBY3 had the same motif structure as YAB2, and ACYABBY4 has no other motifs except the common motif 1 and motif 2 of eight YABBY genes (Fig. S3). In addition, the results from subcellular location predictions showed that all YABBY proteins were most likely to be located in the nucleus, suggesting that they function on the nucleus or directly on the nucleus.

Classification and phylogenetic analyses of AcYABBYs

The phylogenetic analysis results showed that star fruit AcYABBYs can be divided into five subclades, namely, FIL/YAB3 (two members, AcYABBY5 and AcYABBY6), YAB2 (two members, AcYABBY3 and AcYABBY4), YAB5 (two members, AcYABBY1 and AcYABBY8), INO (one member, AcYABBY2), and CRC (one member, AcYABBY7), and the total number of AcYABBYs is far below that in Z. mays (Table S4, Fig. 1). Specifically, the two monocot plants rice and maize have more YAB3-like genes (three and five members, respectively) than the other examined dicotyledons, in which A. thaliana, tomato, grape, and star fruit have six, nine, seven and eight members, respectively (Table S4). Similar to the other two tested species with fleshy fruit, star fruit had two YAB2-like members, which is less than that in maize (five members) and rice (three members) but more than that in A. thaliana (only one member). Interestingly, the YAB5-like genes were absent in the monocot plants rice and maize but present in other plants, and their number seemed to be conservative (one or two). Similarly, the FIL/YAB3 subclade consisted of two parts: monocots and dicots. Among all the tested plants, both tomato and maize have two CRC-like genes, there was only one CRC-like gene present in the other plants, and only one INO-like gene was found in all examined species.

Structure and chromosome distribution of YABBY genes in star fruit

A total of eight AcYABBYs in star fruit were distributed on five chromosomes (Fig. S4). Chr.02, Chr.06, and Chr.08 contained equal numbers of AcYABBY (only one member), whereas Chr.04 contained two AcYABBYs and Chr.09 contained three AcYABBYs members. In addition, a pair of homologous genes, AcYABBY1 and AcYABBY8, was located on the same chromosome, but two other homologous gene pairs (AcYABBY3 and AcYABBY4, AcYABBY5 and AcYABBY6) were located on different chromosomes. Two floral-specific genes (INO-like gene AcYABBY2 and CRC-like gene AcYABBY7) were located on the same chromosome.

The results suggested that all of the AcYABBYs have introns and exons. YAB2-like (AcYABBY3 and AcYABBY4) and CRC-like (AcYABBY7) AcYABBYs had fewer introns than the other AcYABBYs, while the lengths of introns in the YAB2-like (AcYABBY3 and AcYABBY4) AcYABBYs were greater than those in other AcYABBYs, at over 2 kb (Fig. S5). Among the eight AcYABBYs, there were six exons in YAB2-like (AcYABBY3 and AcYABBY4), CRC-like (AcYABBY7), and one YAB1-like member AcYABBY8, compared to the seven exons in other AcYABBYs (Fig. S5). In addition, regardless of whether it was YAB1-like, YAB2-like, or YAB5-like, the gene structure showed a high degree of similarity in spite of slight differences existing in the introns and exons, which may be because of their phylogenetic relationship. AcYABBYs contained highly conserved gene structures, especially regarding orthologous genes (Fig. S5).

Secondary structure analysis of YABBY proteins

The secondary structure analysis showed that all the star fruit AcYABBY proteins contained α-helix, strand, and coil structures (Fig. S6). Homologous genes usually have generally similar secondary structures and stable secondary structure types at nearly the same position as the base pair (bp). For example, AcYABBY3 and AcYABBY4 had a helix and a strand between 117 bp and 123 bp, and between 12 bp and 15 bp, respectively.

Cis-acting elements located in YABBY gene promoters

To explore the possible regulatory functions of AcYABBYs, the 2,000 bp promoter regions of eight AcYABBYs were submitted to the PlantCARE database for the identification of putative cis-elements. A total of 14 types of elements were identified: the zein metabolism regulation element, light responsiveness element, MeJA-responsiveness element, anaerobic induction element, abscisic acid responsiveness element, salicylic acid responsiveness element, meristem expression element, seed-specific regulation element, circadian control element, low-temperature responsiveness element, auxin-responsive element, endosperm expression element, defense and stress responsiveness, and wound-responsive element (Fig. S7). All AcYABBYs contained multiple promoters and all contained light responsiveness element, while most subfamilies obtained salicylic acid responsiveness element, except for the YAB3 subfamily (Fig. S7). Most of the AcYABBYs contained more than five elements, the smallest was AcYABBY5 with three and the largest was AcYABBY6 with nine elements (Fig. S7). Interestingly, the abscisic acid responsiveness element, and the anaerobic induction element were absent in AcYABBY5 and AcYABBY3, respectively, but their homologous genes and other AcYABBYs all contained these elements (Fig. S7). In addition, some promoter elements were specific to a single AcYABBY, including wound-responsive element specific to AcYABBY7, endosperm expression element specific to AcYABBY6, and low-temperature responsiveness element specific to AcYABBY6 (Fig. S7). Moreover, some elements were specific to two AcYABBYs, for example, the seed-specific regulation element specific to AcYABBY4 and AcYABBY6, the auxin-responsive element specific to AcYABBY1 and AcYABBY6, the circadian control element specific to AcYABBY3 and AcYABBY4, and the zein metabolism regulation element specific to AcYABBY3 and AcYABBY8 (Fig. S7).

Synteny and evolutionary analysis of AcYABBYs

We investigated collinear relationships among the orthologous YABBY genes from star fruit, grape, and A. thaliana to find putative evolutionary events. There are six, seven, and eight YABBY genes in A. thaliana, grape, and star fruit, respectively. The results showed that all of YABBY genes showed one-to-one corresponding relationship in A. thaliana - star fruit, and star fruit-grape (Fig. 2).

In the process of evolution, genes usually face various selection pressures, such as positive selection (Ka/Ks > 1), neutral selection (Ka/Ks = 1) and purifying selection (Ka/Ks < 1), understanding the selection pressure of genes plays an important role in studying the evolution of genes (Khan et al., 2019). To explore evolutionary constraints among the homologous AcYABBYs, we calculated Ka, Ks, and the ratio for Ka/Ks. In our study, 28 gene pairs were identified using KaKs_Calculator. Ka, Ks, and Ka/Ks values between each pair are listed in Table S5. The Ka of AcYABBYs ranged from 0.14116 to 0.591704, the Ks from 0.877205 to 3.51706, and the Ka/Ks from 0.104914 to 0.313231, suggesting that all of them had undergone strong purifying selection (Table S5).

Expression patterns of YABBY genes in star fruit

To further elucidate the function of the star fruit YABBY gene family, eight AcYABBYs were used to produce an expression profile heat map using Tbtools, and three biological replications were performed for each tissue to obtain accurate expression levels. AcYABBY2 expression was extremely low and almost undetectable in all tissues (Fig. 3). AcYABBY7 was at an extremely low expression level in the leaf and fruit at different developmental stages but exhibited a relatively high expression level in the inflorescences, bud (Fig. 3). Overall, expression levels of YAB2-like, YAB5-like, and FIL-like genes in reproductive organs were slightly higher than those in the leaf, and F-DAP40 expression levels of nearly all AcYABBYs expect for AcYABBY6 and AcYABBY7 were either lowest or highest during fruit development (Fig. 3).

To further validate the accuracy of transcriptome sequencing, all AcYABBYs in the star fruit during the developmental stages of fruit were subjected to RT-qPCR. The expression levels of AcYABBYs can be divided into four categories based on the trend of relative expression (Fig. 4). AcYABBY1, AcYABBY2, AcYABBY5, AcYABBY6, and AcYABBY8 showed a trend of an initial decrease followed by an increase, and the expression level of F-DAP20 was highest. AcYABBY7 showed a diametrically opposite trend. Expression levels of AcYABBY4 gradually increased during fruit development, AcYABBY3 also showed a trend of an initial decrease followed by an increase, and the expression levels of F-DAP60 was highest (Fig. 4).