Genome-Wide Identiﬁcation and Expression Pattern of the GRAS Gene Family in Pitaya ( Selenicereus undatus L.)

Simple Summary: The GRAS gene family plays a critical role in regulation of growth, defense, light and hormone responses. We identiﬁed 45 GRAS genes in the pitaya genome and categorized them into nine respective subfamilies: PAT1, SHR, LISCL, HAM, SCR, RGL, LAS, DELLA and SCL-3. Among these 45 HuGRAS family members, we reported nine candidate genes that played key roles in the growth and development of the pitaya plant


Domain Analysis of HuGRAS Proteins
All 45 HuGRAS protein sequences were subjected to finding of conserved domains using the NCBI conserved domain tool, https://www.ncbi.nlm.nih.gov/Structure/cdd/ wrpsb.cgi (accessed on 13 October 2022), against Pfamv34.0-19178pSSMs. The retrieved data was used to draw the structure of the GRAS domain using TBTools [59]. A motif finder, https://www.genome.jp/tools/motif/ (accessed on 17 October 2022), was also used to compute the concerned motifs.

Pattern and Distribution of Conserved Motifs
The MEME suite program was used to analyze the 45 HuGRAS genes and identify conserved motifs. The default parameters were used to identify the maximum 10 motifs at https://meme-suite.org/meme/tools/meme (accessed on 13 October 2022).

Network Analysis of HuGRAS Proteins
All 45 HuGRAS genes were subjected to collection of their interactions with other genes using the pitaya genome website, http://www.pitayagenomic.com/coexpression (accessed on 30 October 2022) [56]. The Cytoscape tool was used to build the network using information from all of the identified and interacting genes [62].

Cis-Acting Element Analysis in HuGRAS Promoter Sequences
Promoter sequences were retrieved from the pitaya genome file using TBTools (2000 bp upstream of the start codon) [59]. The PlantCARE database was used to retrieve cisacting regulatory elements [63].

Plant Materials
The "Shuangse Dahong" pitaya variety was used in this experiment, and flower buds were collected from the germplasm resource of Hainan-Shengda Modern Agriculture Development Company, Qionghai, Hainan, China. All plants were grown in field conditions.

RNA Isolation and Real-Time Quantitative PCR Expression Analysis
Utilizing the RNAprep Pure Plant Kit, total RNA was extracted (TIANGEN, Beijing, China). The plant material used for RNA extraction included stems (one-month-old stems, one-year-old stems and two-year-old stems, designated S1Ms, S1Ys and S2Ys, respectively), flower buds (FBs), pericarp (PeriC) and pulp. A NanoDrop 2000C spectrophotometer was used to measure the concentration of the samples (Thermo Fisher Scientific, Waltham, MA, USA). DNase I was used to remove genomic DNA from a total of 1 g of RNA from each sample before being utilized as a template for reverse transcription to create the desired amount of cDNA (QuantiTect Reverse Transcription Kit; Qiagen, Shanghai, China). The RNA sample for each qRT-PCR was standardized using the actin-gene-expression level in S. undatus L. Three biological and three technical replications of each sample were used in the qPCR, with ACTIN serving as the internal control. The SYBER Green Master Mix (Novogene, Shanghai, China) was used, along with the LightCycler 480 real-time PCR system (Applied Biosystem, St. Louis, MO, USA). qRT-PCR results were analyzed using the double-delta CT method [64,65].

Genome-Wide Identification of the GRAS Family in Pitaya
Through genome-wide analysis, 45 candidate genes were retrieved from the pitaya genome, and these genes were designated HuGRAS-1 to HuGRAS-45. Basic physical and chemical properties of the HuGRAS genes, including each gene's chromosome number, position on the chromosome, CDS length, protein length, protein molecular weight, isoelectric point (pI) and GRAVY, are summarized in Table 1. Protein length and molecular weight varied greatly, ranging from 97 (HU08G00229.1) to 809 AA (HU01G00472.1), and molecular weight ranged from 15-95 kDa. All 45 pitaya HuGRAS proteins were predicted to be hydrophilic because their representative GRAVY values were less than 0, ranging from −0.006 (HU06G00376.1) to −0.436 (HU02G01570.1). All HuGRAS proteins comprised varying degrees of pI values, ranging from 5.4 (HU10G00709.1) to 9.5 (HU08G00229.1), with an average value of 7.2.

Phylogenetic Analysis of the GRAS Gene Family
To construct a phylogenetic tree, protein sequences from the characterized species were retrieved from previous studies, including those on Arabidopsis, Medicago truncatula, tomatoes, rice, soybeans and maize. Characterized GRAS protein sequences from

Phylogenetic Analysis of the GRAS Gene Family
To construct a phylogenetic tree, protein sequences from the characterized species were retrieved from previous studies, including those on Arabidopsis, Medicago truncatula, tomatoes, rice, soybeans and maize. Characterized GRAS protein sequences from six species and the Phytozome website were used to retrieve all GRAS proteins from their respective genomes. Collectively, 380 protein sequences were used to draw a phylogenetic tree, including the 45 HuGRAS protein sequences from the pitaya genome and 335 protein sequences from the other six species (Supplementary File S1). In the resulting phylogenetic tree, the GRAS genes were divided into nine subfamilies: PAT1, SHR, LISCL, HAM, SCR, RGL, LAS, DELLA and SCL3 ( Figure 2). All 45 HuGRAS proteins were grouped as follows: twelve in PAT1; ten in LISCL; five in HAM; four each in SHR, SCL3 and SCR; three in DELLA; two in LAS; and one in RGL. The PAT1 subfamily contained 12 types of HuGRAS protein and was the largest subfamily of GRAS protein, while RGL had only one type of HuGRAS protein and was one of the smallest subfamilies of GRAS protein.

HuGRAS-Protein Sequence Alignments and Conserved Motifs
A conserved-motif analysis of each GRAS protein was carried out using the MEME tool. Ten conserved motifs were identified. The C terminal regions of the HuGRAS proteins contained highly conserved domains. Motif 2, motif 6 and motif 7 were found in almost all of the HUGRAS proteins. However, motif 5 was not found in HuGRAS-35, HuGRAS-39, HuGRAS-42 or HuGRAS-45 ( Figure 3). Most of the GRAS proteins carried similar motifs within the group, with very few motif differences. These findings also helped us to understand the close evolutionary relationships of the same protein group. The known-motif of the amino-acid-sequence is exhibited in Figure S1.

HuGRAS-Protein Sequence Alignments and Conserved Motifs
A conserved-motif analysis of each GRAS protein was carried out using the MEME tool. Ten conserved motifs were identified. The C terminal regions of the HuGRAS proteins contained highly conserved domains. Motif 2, motif 6 and motif 7 were found in almost all of the HUGRAS proteins. However, motif 5 was not found in HuGRAS-35, HuGRAS-39, HuGRAS-42 or HuGRAS-45 ( Figure 3). Most of the GRAS proteins carried similar motifs within the group, with very few motif differences. These findings also helped us to understand the close evolutionary relationships of the same protein group. The knownmotif of the amino-acid-sequence is exhibited in Figure S1.

Gene Structure and Distribution of HuGRAS Genes on Chromosomes
To find the gene structures, the intron and exon structures of all of the HuGRAS genes were aligned ( Figure S2). Among all 45 HuGRAS genes, most of the HuGRAS gene sequences showed two sequences of introns and one sequence of exons. The majority of the HuGRAS genes had a similar pattern of exons, indicating the phylogeny and evolution of their gene family, except for HuGRAS-8 and HuGRAS-42.
The HuGRAS genes were physically located on 11 chrs in the pitaya genome. All GRAS genes were mapped on pitaya chrs based on the information available at the pitaya genome website, http://www.pitayagenomic.com/ (accessed on 1 September 2022). The chr length and position of each gene of the pitaya genome is presented in Supplementary File S2. Forty-five HuGRAS genes were unevenly distributed on 11 chrs. Most of the HuGRAS genes were found on chr 02 and chr 08. Chr 10 had one GRAS gene but chr 09 had no HuGRAS genes ( Figure 4).

Gene Structure and Distribution of HuGRAS Genes on Chromosomes
To find the gene structures, the intron and exon structures of all of the HuGRAS genes were aligned ( Figure S2). Among all 45 HuGRAS genes, most of the HuGRAS gene sequences showed two sequences of introns and one sequence of exons. The majority of the HuGRAS genes had a similar pattern of exons, indicating the phylogeny and evolution of their gene family, except for HuGRAS-8 and HuGRAS-42.
The HuGRAS genes were physically located on 11 chrs in the pitaya genome. All GRAS genes were mapped on pitaya chrs based on the information available at the pitaya genome website, http://www.pitayagenomic.com/ (accessed on 1 September 2022). The chr length and position of each gene of the pitaya genome is presented in Supplementary File S2. Forty-five HuGRAS genes were unevenly distributed on 11 chrs. Most of the HuGRAS genes were found on chr 02 and chr 08. Chr 10 had one GRAS gene but chr 09 had no HuGRAS genes (Figure 4).

Identification of Cisacting Elements in HuGRAS Promoter Sequences
To identify the biological functions (stress response, growth and development) of the HuGRAS genes, all 45 HuGRAS gene sequences (2000 bp upstream of start codon) were selected for cis-element analysis using the PlantCARE web tool (Supplementary File S3). In total, 17 cis-elements were recorded in this study ( Figure S3). The cis-regulatory elements of 45 GRAS proteins are shown in Figure S3. Nine genes that exhibited higher expression among the GRAS gene family in the pitaya plant exhibited various cis-acting regulatory elements, as shown in Figure 7. Fourteen cis-acting elements were categorized into Figure 6. Protein interaction network. Green and red genes are designated as hub genes because they interact with more than 10 genes. Different colors show the interactions of the genes as follows: green (16)(17)(18)(19)(20), red (11)(12)(13)(14)(15), yellow (6-10) and gray (2-5).

Identification of Cisacting Elements in HuGRAS Promoter Sequences
To identify the biological functions (stress response, growth and development) of the HuGRAS genes, all 45 HuGRAS gene sequences (2000 bp upstream of start codon) were selected for cis-element analysis using the PlantCARE web tool (Supplementary File S3). In total, 17 cis-elements were recorded in this study ( Figure S3). The cis-regulatory elements of 45 GRAS proteins are shown in Figure S3. Nine genes that exhibited higher expression among the GRAS gene family in the pitaya plant exhibited various cis-acting regulatory elements, as shown in Figure 7. Fourteen cis-acting elements were categorized into four groups: lightresponsive elements, growth and development elements, stress-and defense-responsive elements and hormone-responsive elements.

Expression of HuGRAS Genes at Developmental Stages of Pitaya
To confirm the expressions of the predicted genes from the transcriptome data, we conducted qRT-PCR for HuGRAS  and HuGRAS-37. We designed primers for the 12 candidate genes, categorized in six GRAS subfamilies (Supplementary File S4). The results thereof exhibited that HuGRAS-1 ,  HuGRAS-7, HuGRAS-12, HuGRAS-18, HuGRAS-25, HuGRAS-34, HuGRAS-35, HuGRAS-41 and HuGRAS-37 showed higher levels of expression across the tissues (Figure 8). The expression levels of the HuGRAS members varied widely in different tissues. The HuGRAS-1 gene, categorized in the DELLA subfamily, was significantly expressed across the tissues, including the stems, FBs and pericarp. However, relatively weaker expression was observed in the pulp of the fruit. Among the PAT1 subfamily members, HuGRAS-34, HuGRAS-35 and HuGARS-41 exhibited strong expression in the plant tissues as compared to HuGRAS-6, which exhibited weak expression in the pericarp and the pulp. HuGRAS-7, which belongs to the SCL-3 subfamily, was expressed at a low level in the one-month-old stem cells but abundant in other tissues. HuGRAS-12, a gene categorized in the SCR subfamily, was expressed at a higher level in other tissues than the flower buds. HuGRAS-21 and HuGRAS-29, members of the LISCL subfamily, were expressed at lower levels than the HuGRAS-18 and HuGRAS-25 grouped in the same subfamily, which were expressed at higher levels in the flower buds, the pericarp and the pulp of the pitaya plant. HuGRAS-37, grouped into the HAM subfamily, was highly expressed in the flower buds but weakly expressed in the one-month-old stems. Nine genes, which were categorized into six subfamilies, exhibited higher expression levels and might play key roles in the growth and development of the pitaya plant.

Expression of HuGRAS Genes at Developmental Stages of Pitaya
To confirm the expressions of the predicted genes from the transcriptome data, we conducted qRT-PCR for HuGRAS-1, HuGRAS-6, HuGRAS-7, HuGRAS-12, HuGRAS-18, HuGRAS-21, HuGRAS-25, HuGRAS-29, HuGRAS-34, HuGRAS-35, HuGRAS-41 and HuGRAS-37. We designed primers for the 12 candidate genes, categorized in six GRAS subfamilies (Supplementary File S4). The results thereof exhibited that HuGRAS-1, HuGRAS-7, HuGRAS-12, HuGRAS-18, HuGRAS-25, HuGRAS-34, HuGRAS-35, HuGRAS-41 and HuGRAS-37 showed higher levels of expression across the tissues (Figure 8). The expression levels of the HuGRAS members varied widely in different tissues. The HuGRAS-1 gene, categorized in the DELLA subfamily, was significantly expressed across the tissues, including the stems, FBs and pericarp. However, relatively weaker expression was observed in the pulp of the fruit. Among the PAT1 subfamily members, HuGRAS-34, HuGRAS-35 and HuGARS-41 exhibited strong expression in the plant tissues as compared to HuGRAS-6, which exhibited weak expression in the pericarp and the pulp. HuGRAS-7, which belongs to the SCL-3 subfamily, was expressed at a low level in the one-month-old stem cells but abundant in other tissues. HuGRAS-12, a gene categorized in the SCR subfamily, was expressed at a higher level in other tissues than the flower buds. HuGRAS-21 and HuGRAS-29, members of the LISCL subfamily, were expressed at lower levels than the HuGRAS-18 and HuGRAS-25 grouped in the same subfamily, which were expressed at higher levels in the flower buds, the pericarp and the pulp of the pitaya plant. HuGRAS-37, grouped into the HAM subfamily, was highly expressed in the flower buds but weakly expressed in the one-month-old stems. Nine genes, which were categorized into six subfamilies, exhibited higher expression levels and might play key roles in the growth and development of the pitaya plant. The X-axis represents the plant tissues, including one-month-old stem (S1M), one-year-old stem (S1Y), two-year-old stem (S2Y), flower bud (FB), pericarp (PeriC) and pulp. Error bars represent the standard deviations for the three replicates.

Discussion
Pitaya (S. undatus L.) is a tropical fruit, typically cacti, evergreen, and consists of cladodes (a modified stem replaces the leaves for photosynthesis function) which perform its functioning as a leaf. The flowers and fruits are edible, and the pericarp and pulp of S. undatus are white in color. The fruits are highly enriched with polyphenols, tannis, betalains and nonbetalainic and antioxidant compounds [66]. Due to the importance of the pitaya tropical fruit, the present study was carried out to explore the growth and the developmental process of the plant. The GRAS TF is being explored in other crop species, such as Arabidopsis [55], Medicago truncatula [48], pepper [54], cotton [50], soybeans [27], tomatoes [67], Chinese cabbage [16] and tropical fruit such as litchi [28], but we could not find any research studies about the GRAS gene family in S. undatus L. GRAS proteins have been recognized as important TF, playing different functions in plant growth and development, including patterning of roots and shoots, responses to various kinds of stresses, stem-cell initiation and maintenance [35], light signaling and the gibberellic-acid signaltransduction pathway [21,68].
With the availability of the pitaya reference genome [2] and pitaya tissue expression data via the pitaya genome and multiomics database [56], we performed a genome-wide identification of the GRAS gene family members in the pitaya genome. In the current study, we found 45 GRAS gene family members in this genome, named HuGRAS-1 to HuGRAS-45; they were widely distributed on 11 chromosomes (Figure 4). Most of the HuGRAS genes were found were on the ends of these chromosomes, which is in accordance with other plant species, such as watermelon, potatoes, rice and Arabidopsis [69]. The conserved motif structures ( Figure 3) and HuGRAS gene sequences ( Figure S2) exhibited the same pattern of conserved motifs and exon-intron sequences, respectively, The X-axis represents the plant tissues, including one-month-old stem (S1M), one-year-old stem (S1Y), two-year-old stem (S2Y), flower bud (FB), pericarp (PeriC) and pulp. Error bars represent the standard deviations for the three replicates.

Discussion
Pitaya (S. undatus L.) is a tropical fruit, typically cacti, evergreen, and consists of cladodes (a modified stem replaces the leaves for photosynthesis function) which perform its functioning as a leaf. The flowers and fruits are edible, and the pericarp and pulp of S. undatus are white in color. The fruits are highly enriched with polyphenols, tannis, betalains and nonbetalainic and antioxidant compounds [66]. Due to the importance of the pitaya tropical fruit, the present study was carried out to explore the growth and the developmental process of the plant. The GRAS TF is being explored in other crop species, such as Arabidopsis [55], Medicago truncatula [48], pepper [54], cotton [50], soybeans [27], tomatoes [67], Chinese cabbage [16] and tropical fruit such as litchi [28], but we could not find any research studies about the GRAS gene family in S. undatus L. GRAS proteins have been recognized as important TF, playing different functions in plant growth and development, including patterning of roots and shoots, responses to various kinds of stresses, stem-cell initiation and maintenance [35], light signaling and the gibberellic-acid signal-transduction pathway [21,68].
With the availability of the pitaya reference genome [2] and pitaya tissue expression data via the pitaya genome and multiomics database [56], we performed a genome-wide identification of the GRAS gene family members in the pitaya genome. In the current study, we found 45 GRAS gene family members in this genome, named HuGRAS-1 to HuGRAS-45; they were widely distributed on 11 chromosomes (Figure 4). Most of the HuGRAS genes were found were on the ends of these chromosomes, which is in accordance with other plant species, such as watermelon, potatoes, rice and Arabidopsis [69]. The conserved motif structures ( Figure 3) and HuGRAS gene sequences ( Figure S2) exhibited the same pattern of conserved motifs and exon-intron sequences, respectively, suggesting that these genes may have similar functions to those reported in previous studies [21].
In accordance with phylogenetic analysis, we compared the 45 HuGRAS gene sequences with 335 sequences of GRAS proteins from maize, soybeans, Medicago truncatula, rice, Arabidopsis and tomatoes. HuGRAS genes were divided into nine subfamilies based on clade support values: PAT1, SHR, LISCL, HAM, SCR, RGL, LAS, DELLA and SCL3 ( Figure 2). Each subfamily carried varying numbers of HuGRAS genes, and the PAT1 subfamily contained the largest number of HuGRAS genes. The protein sequences and differential expression profiles of pitaya tissues aid in the identification of the genes that play key roles in growth and development. Expression and network analysis provide a clue to locating genes that exhibit high levels of expression ( Figure 5) and interact with many other genes ( Figure 6). With the help of expression analysis and network analysis, 12 selected genes were categorized into their respective subfamilies, as predicted in the phylogenetic tree ( Figure 2). All genes were placed in their respective GRAS families: HuGRAS-1 in the DELLA subfamily; HuGRAS-6, HuGRAS-34, HuGRAS-35 and HuGRAS-41 in the PAT1 subfamily; HuGRAS-7 in the SCL-3 subfamily; HuGRAS-12 in the SCR subfamily; HuGRAS-18, HuGRAS-21, HuGRAS-25 and HuGRAS-29 in the LISCL subfamily; and HuGRAS-37 in the HAM subfamily. qRT-PCR was carried out for the predicted gene subfamilies ( Figure 8) to confirm their expressions in different stages of the plant. The HuGRAS-1 gene, categorized in the DELLA subfamily, was significantly expressed across the tissues, as the DELLA subfamily is involved in the growth and development of the plant [41]. In the absence of gibberellic acid, DELLA proteins interact with light-responsive TFs, including phytochrome-interacting factors (PIFs), to form inactive complexes [70], while higher expression of gibberellic acid degrades DELLA proteins and initiates the growth rate [71]. Our network analysis revealed that the HuGRAS-1 gene interacts with almost 20 other proteins, so our results are consistent with the prediction of Hirsch and Oldroyd, 2009 [41] that DELLA proteins interact with other PIF families and make complexes with them. This DELLA-PIF TF complex is possibly competitive but a common mechanism for DELLAs to make complexes for light-and gibberellic-acid-signaling to alter environmental conditions [41]. DELLA proteins also regulate immune responses by regulating the jasmonic-and salicylic-acid pathways. The PAT1 subfamily members, HuGRAS-34, HuGRAS-35 and HuGARS-41, exhibited strong expressions in plant tissues as compared to HuGRAS-6, which was weakly expressed. PAT1 is a specific member of the GRAS family that interacts with light signaling via phytochrome A to regulate the plant developmental process, including de-etiolation and hypocotyl elongation [45]. HuGRAS-7, which belongs to the SCL-3 subfamily, was expressed at a low level in one-month-old stem cells but exhibited higher expression abundantly in other tissues. The SCL-3 subfamily acts antagonistically, downstream to gibberellic-acid DELLA responses and upstream to gibberellic-acid-biosynthesis pathways, during plant growth and development [47]. HuGRAS-12 was categorized in the SCL-3 subfamily and was expressed at a higher level in other tissues than flower buds. HuGRAS-37 is grouped into the HAM subfamily and was highly expressed in flower buds but weakly expressed in one-month-old stems. The SCR and HAM subfamilies play key roles in root/shoot patterning and cell differentiation in shoot meristem maintenance [72]. HuGRAS-21 and HuGRAS-29, members of the LISCL subfamily, were expressed at lower levels than the HuGRAS-18 and HuGRAS-25 genes of the same subfamily, which were expressed at a higher level in the flower buds, the pericarp and the pulp of the pitaya plant. Higher levels of expression of LISCL genes (HuGRAS-18 and HuGRAS-25) may predict their role in flower development and fruit ripening. The LISCL subfamily of the GRAS protein has been reported to play a role in another development, of Lilium longiflorum L. [73]. Previous studies have revealed that redundancy of relative expression and phytohormones in different parts of the reproductive tissue (panicle) can lead to defects in growth. Similarly, varying expression levels of GRAS-family genes might also have different functions in different pitaya tissues [74,75].
In this study, we analyzed the GRAS TF family in the pitaya plant (S. undatus L.) and six other species, including maize, soybeans, Medicago truncatula, rice, Arabidopsis and tomatoes. A total of 380 GRAS genes were analyzed in this research, in addition to 45 genes that were predicted from the pitaya genome. We categorized these genes into nine subfamilies based on phylogenetics and previous studies of other crops. Among the nine subfamilies of GRAS, few genes showed higher expression in different tissues of pitaya plant. These genes were categorized into six sub-families including DELLA (HuGRAS-1), SCL-3 (HuGRAS-7), PAT1 (HuGRAS-34, HuGRAS-35, HuGRAS-41), HAM (HuGRAS-37), SCR (HuGRAS-12) and LISCL (HuGRAS-18, HuGRAS-25) which may have potential key role in the growth and development of the pitaya plant. Their roles were also confirmed using in silico cis-acting analysis ( Figure S3). As we could see, cis-acting elements, including gibberellin, auxin, ABA, jasmonic-acid and salicylic-acid-responsive elements were abundantly present in the HuGRAS promoters. These genes can be used to study the regulatory pathways of specific plant traits. Positive and negative regulators can be identified from the pathways; then the CRISPR system can be used to produce a transgene-free pitaya plant. Previously, many crops have been improved using the latest genome editing technique [76,77]. Collectively, our results lay a theoretical foundation for the role of GRAS genes in pitaya growth and development. It provides valuable information to improve the pitaya breeding program.

Conclusions
This study is the first comprehensive genome-wide identification of the GRAS gene family in pitaya (S. undatus L.). This research might aid in the interpretation of the GRAS genes function, protein interactions, signaling-pathway regulations and expression patterns in different tissues. The comparative study between the GRAS families of six species, the phylogenetic tree, the expression pattern and the gene network analysis will lay a foundation for the functional characterization of the genes in pitaya. Understanding the possible roles of nine predicted genes (HuGRAS-1, HuGRAS-7, HuGRAS-12, HuGRAS-18, HuGRAS-25, HuGRAS-34, HuGRAS-35, HuGRAS-37, HuGRAS-41) from the six subfamilies of GRAS gene and their expression patterns in different tissues provides insightful information for the development of pitaya fruit's economic, agronomic and ecological benefits. Altogether, the current study is the first report on the GRAS gene family in pitaya tropical fruit. The identification of the genes will assist in clarifying the molecular genetic basis and aid in improving the genotypes in the breeding program.