The GRAS protein sequences in Arabidopsis were downloaded from the TAIR database¹ (Tian et al., 2004) and those in rice and Amborella trichopoda were collected and published by Alberto et al. (Cenci and Rouard, 2017). The GRAS proteins of the three plant species were used as reference queries to identify GRAS proteins in garlic using BLASTP (Chen et al., 2020). Subsequently, the preliminarily identified proteins were subjected to conserved domain region analysis using the CDD program in NCBI² (Lu et al., 2020) and MEME³ (Bailey et al., 2009). After manual verification, only proteins with lengths greater than 200 amino acids and specific GRAS domains were selected for subsequent analysis. The basic features of GRAS proteins, including molecular weight (MW), coding sequence length (CDS), and isoelectric point (pI), were predicted using the ExPASy software⁴ (Gasteiger et al., 2005), and subcellular localizations were predicted by Localizer 1.0.4 (Sperschneider et al., 2017).

The annotation of garlic genome (Sun et al., 2020), locastion of GRAS genes on the chromosomes/scaffolds was visualized using TBtools (v1.098696; Chen et al., 2020). Multiple sequence alignments of all GRAS proteins were conducted using the JalView (version 2.10.3) software to investigate the protein domains (Waterhouse et al., 2009). The conserved motifs were analyzed using two online tools: MEME (Bailey et al., 2009) and batch CD-Search (Lu et al., 2020). Protein structures and motifs were visualized using TBTools (Chen et al., 2020).

Tissue-specific expression pattern of AsGRAS genes was analyzed by extracting the expression data from the reported transcriptome analysis of garlic (Sun et al., 2020) and an expression heatmap was obtained using TBtools (Chen et al., 2020). Zhu et al. (2019) completed the mRNA sequencing of enlarged bulbs from 102 accessions and introduced an associated transcriptomic method to identify genes related to transcript abundance and bulb traits (Zhu et al., 2019). Based on the expression data of these 102 accessions and the associated transcriptomic method, we performed a correlation analysis of bulb traits and transcript abundance of garlic DELLA genes. Orthologous gene pairs between AsGRASs and AtGRASs were identified and the interaction network of orthologous GRAS genes was predicted using the STRING 11.5 database (von Mering et al., 2005).

Approximately 100 garlic cloves (cv. Ershuizao) were planted in a pot and grown in a greenhouse with a 12 h light period. After 1 month, half of the garlic seedlings were treated by spraying exogenous GA3 and the other half served as the control group. Briefly, GA3 (200 mg/l) and sterile distilled water were sprayed on the leaves of treatment and control groups, respectively. Leaf samples were collected in triplicates from both groups after 0, 4, and 16 h of treatment and were pre-ground to a powder in liquid nitrogen. Total RNA extraction was performed using SteadyPure plant RNA extraction kit AG21019 [Accurate Biotechnology (Hunan) Co., Ltd.] and resultant RNA population was qualified by gel electrophoresis followed by first-strand cDNA biosynthesis [Evo M-MLV RT premix for qPCR AG11706, Accurate Biotechnology (Hunan) Co., Ltd.]. qRT-PCR analysis was performed using CFX96 Real-Time PCR Detection System (Bio-Rad Laboratories, United States) with SYBR green Premix Pro Taq HS qPCR Kit AG11701 [Accurate Biotechnology (Hunan) Co., Ltd.]. Primer sequences are listed in Supplementary Table S1.

A total of 61 GRAS candidate genes were obtained from the garlic genome according to the BLAST results. Of these, 15 encoded a protein with <200 amino acids residues and/or without GRAS domains and were filtered. Finally, 46 GRAS members were identified which typically comprised of five conserved domains: LHRI, VHIID, LHRII, PFYRE, and SAW (Supplementary Figure S1; Supplementary Table S2). The 46 identified GRAS genes were distributed on seven of eight garlic chromosomes (except for chromosome 1) and their distribution varied among chromosomes, ranging from 2 to 15 genes per chromosome (Figure 1).

The length of the putative GRAS proteins ranged from 202 to 705 amino acids and the predicted molecular weights ranged from 22 to 80 kDa (Table 1). The theoretical pI varied from 4.87 to 9.76. Except for Asa4G03041.1, all the proteins had a hydropathicity value of <0, indicating the hydrophilic nature of these GRAS proteins. Additionally, only seven proteins showed an instability index value of <40.0 (i.e., threshold value; Table 1), suggesting that most of the garlic GRAS proteins probably had an unstable structure. Subcellular localization prediction identified between 16 and 2 GRAS proteins located in the nucleus and chloroplast, respectively (Table 1). Interestingly, the Asa3G01451.1 protein was predicted with position signals in both the nucleus and chloroplast.

Phylogenetic analysis was performed for 46 AsGRAS members, together with 123 known GRAS proteins from Arabidopsis (33 members), O. sativa (56 members), and A. trichopoda (34 members). The results indicated that 46 AsGRAS members were assigned to 16 of the 17 known GRAS subfamilies (Figure 2). The protein structure investigation indicated that 46 AsGRAS proteins had conserved domains of the GRAS or GRAS superfamily, comprising of a conserved motif structure with LHRI-VHIID-LHRII-PRYRE-SAW domains (Table 1; Figure 3A). However, protein motifs showed larger differences between members from different subfamilies than from the same subfamily, suggesting a diversification of motif sequences among the garlic GRAS subfamily (Figure 3B; Supplementary Figure S2). The details of structural organization revealed that 32 of 46 AsGRAS genes had no introns; in particular, all genes of the LISCL, DELLA, SCL3, SHR, and HAM subfamilies lacked any introns, whereas all genes of the PAT1, RAM1, and SCR subfamilies had more than two exons (Figure 3C).

The expression patterns of 46 AsGRAS genes in seven tissues and eight bulb-developmental stages were investigated and revealed significant differences (Figure 4). Notably, DELLA-like Asa4G02090.1 showed the highest expression level in all the tissues except roots, whereas nine AsGRAS genes were not expressed at all in the 15 samples including, Asa6G00928.1, Asa4G03041.1, Asa6G04430.1, Asa6G06640.1, Asa2G00446.1, Asa8G00797.1, Asa6G02143.1, Asa6G00303.1, and Asa7G03967.1.

To understand the biological function of AsGRAS TFs further, their protein–protein interactions were predicted using the ortholog-based method. The results showed that 46 AsGRAS members were orthologs of 20 Arabidopsis GRAS proteins that constituted an interacting network (Figure 5). According to the prediction, AsGRAS proteins interacted with PIF, BZR1, BIN2, and GID1 proteins, which are involved in GA and brassinolide signal transduction mechanism, and lighting-responsive PHYA and PHYB proteins. Furthermore, some members of SCL subfamilies showed a putative interaction with cell cycle regulatory DPA and E2F proteins, indicating a potential role of these garlic GRAS proteins in the regulation of cell cycle.

As a subfamily of the AsGRAS family, the DELLA subfamily has been frequently studied as these proteins are the central repressors of gibberellin (GA) response (Park et al., 2013). Bioinformatics prediction identified five garlic DELLA members as Asa2G00237.1, Asa2G00240.1, Asa4G02090.1, Asa6G05069.1, and Asa6G07138.1. Of these DELLA-like genes, Asa2G00237.1 and Asa2G00240.1 showed a tandem repeat distribution in the genomic regions, with a complete identical coding sequence; whereas Asa6G05069.1 and Asa6G07138.1 displayed almost identical sequences in their coding regions.

The bulb is the main consumed organ of garlic and consists of several cloves. GA3 treatment can increase the number of cloves in garlic (Liu et al., 2019). In this study, we investigated the expression response of DELLA-like genes to GA3 treatment. Because a complete identical coding sequence makes a challenge to distinguish their transcripts, we performed the expression analysis for these two genes as a whole. Consequently, we found that Asa4G02090.1 and Asa2G00237.1/Asa2G00240.1 showed distinct expression response to GA3 treatment (Figure 6; Supplementary Table S3). Interestingly, the expression of these three genes in response to GA3 treatment was observed during the transition from dark to light, but not in constant light or dark conditions. Furthermore, the associated transcriptomic analysis revealed a significant correlation between the transcript abundance of DELLA-like genes, Asa4G02090.1, and bulb weight and diameter (p = 0.003784 and 0.001473, respectively). Also, the expression of Asa2G00237.1/Asa2G00240.1 displayed a significant association with the clove number trait in 102 garlic accessions (p = 0.03147; Figure 7), indicating a role of these three DELLA-like genes in bulb growth. Collectively, our results indicate that GA3 could result in varied expression of DELLA-like Asa4G02090.1 and Asa2G00237.1/Asa2G00240.1, thereby influencing bulb growth.

The garlic genome is one of the most complex genomes with a large size (16.9 Gb), high heterozygosity (1.69%), and a high ratio of repetitive sequences (91.3%). Recently, the assembly of garlic genome was completed using a repertoire of five advanced sequencing methodologies (Sun et al., 2020), which made it feasible to identify the genes important for garlic growth and development in Allium crops. Consequently, several gene families, such as GH19, PR, and KNOX, have been systemically characterized (Anisimova et al., 2021; Filyushin et al., 2021; Zhang et al., 2021). The GRAS family is one of the most important transcription factor families that participate widely in the plant growth regulation and development. However, GRAS genes were not identified in garlic so far. In the present study, we identified 46 GRAS members in the garlic genome. The GRAS members were slight fewer than in rice (57; Tian et al., 2004) and maize (86; Guo et al., 2017), indicating a slight contraction of this gene family in garlic genome. According to the difference of grouping criterion, the number of subfamily varied across previous studies, such as eight groups identified in Arabidopsis, rice, and maize (Tian et al., 2004; Guo et al., 2017), and 17 in the report of Cenci and Rouard (2017). To obtain a fine classification, this study analyzed the distribution of garlic GRAS members in 17 reported subfamilies and revealed that 16 of these 17 groups had garlic members. Furthermore, we systematically characterized these GRAS regulators by studying their sequence alignment, gene structure, expression, chromosomal distribution, protein domains encoded by them and their respective interactions, and subcellular localizations. The identification and characterization of AsGRAS members provide an important basis for further investigation of their function in the future.

The garlic bulb is composed of several cloves that are essentially buds in morphology. Previous studies have shown that GAs play a key role in bud formation (Zhang et al., 2020) and the DELLA protein is an essential repressor of gibberellin signal transduction (Ikeda et al., 2001; Sun and Gubler, 2004; Park et al., 2013). For example, SLR1 is the DELLA protein of rice which promotes bud outgrowth by inhibiting the degradation of MOC1 protein, leading to an increase in tiller number (Liao et al., 2019). Similarly, StGA20 induces tuber formation by reducing gibberellin activity in potatoes (Carrera et al., 2000). In addition, exogenous application of GA3 could increase tillers in Welsh onion (Yamazaki et al., 2015) and clove numbers in garlic (Liu et al., 2019). Therefore, gibberellin has a potential role in the development of garlic bulbs. However, the mechanism that gibberellin regulates bulb formation remains unclear. In the present study, we identified three DELLA genes, Asa4G02090.1 and Asa2G00237.1/Asa2G00240.1, whose expression was observed in response to exogenous GA3 treatment and was associated with bulb traits in 102 garlic accessions. These results indicate the possible role of these three DELLA-like genes in the bulb growth. This observation provides an important evidence to explain a previous finding that spraying exogenous GAs could increase the clove number in garlic (Liu et al., 2019).

Light and gibberellins (GAs) are two essential signals that trigger plant developmental processes and their signaling pathways show great overlap (Feng et al., 2008). Light signals can promote the accumulation of DELLA proteins by reducing the GA levels (Achard et al., 2007). In Arabidopsis and potato, the expression of GA 20-oxidase responds to photoperiod to regulate the GA biosynthesis in vivo (Carrera et al., 2000; Porri et al., 2014). PIL5 is a light-labile bHLH TF that interacts with phytochrome, directly binds to the promoter of DELLA, and increases the expression of the DELLA gene in the dark (Oh et al., 2007). In garlic, long daylight can promote bulb enlargement (Wu et al., 2016). GA3 treatment dramatically stimulates lateral bud formation but inhibits the growth of garlic plants and bulbs (Liu et al., 2019). Probably, enlarged size and growth of garlic bulbs are regulated by gibberellin signaling and photoperiod. However, the combined underlying mechanism of gibberellin and daylight in garlic bulb development remains poorly understood. In this study, although none of garlic DELLA-like genes showed an expression response to light treatment, the GA3 response of three DELLA-like genes were promoted in the dark environment, indicating that a potential cross-talk between the signals of GA3 and light. These findings provide insights into the response and adaptation of garlic crops to environmental changes.