The latest version of the sequence of genome DNA, CDS (coding sequence), and the protein of L. chinense was downloaded from CNGBdb (CNA0002404, China National GeneBank DataBase). We downloaded GRAS family HMM (hidden Markov model) profile PF03514 from the Pfam database (https://pfam.xfam.org) to search candidate GRAS genes by using HMMER 3.3. Furthermore, we used the de novo transcriptome data of somatic embryos in hybrid Liriodendron (data not published) to screen candidate GRAS genes and map them with the results. The protein sequences of candidate GRAS genes were uploaded to the NCBI protein database, which was used as initial queries by BLASTp for further verification. Finally, we checked the sequence and annotation by using the online Softberry software (https://www.softberry.com/) and CD-search (https://www.ncbi.nlm.nih.govStructure/cdd/cdd.shtml). In addition, we used the online EXPASY software website (https://www.expasy.org/) to analyze the theoretical pI (isoelectric point of protein), molecular weight, and other basic analyses of the GRAS genes’ information.

We used the MEME (https://meme-suite.org/meme/tools/meme) (version 5.3.3) program to identify the conserved motifs of LcGRAS using the classic mode. We increased the minimum width to 10 and increased the number of motifs to 10. We extracted the 2,000-bp length of the DNA sequence, which is upstream from every LcGRAS gene, as the promoter sequence to analyze the cis element by using PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). Then, repeated motifs or motifs without explicit functional annotations in the results were manually deleted and sorted by function.

We aligned the protein sequences of LcGRAS using Clustal Omega. A phylogenetic tree was constructed via the NJ (neighbor-joining) method with 1,000 bootstrap replicates using MEGA7. To further understand the classification of the GRAS gene family, we downloaded the GRAS protein sequences of other several species from PlantTFDB (http://planttfdb.gao-lab.org/), including A. thaliana, A. trichopoda, O. sativa, Picea abies, and Populus trichocarpa. The genomic data of Euryale ferox were downloaded from CoGe (https://genomevolution.org/coge/). EfGRAS was identified by HMM, and the blast was the same as LcGRAS. After alignment, the GRAS data of Magnolia biondii were available in the Dryad database: https://doi.org/10.5061/dryad.s4mw6m947. RaXml (Stamatakis, 2014) was used to construct the ML (maximum likelihood) tree with 1,000 bootstraps, and it was drawn using the online iTOL tool (https://itol.embl.de/). The Ka/Ks of the GRAS gene was calculated using KaKs_calculator 2.0.

We obtained genomic data at the chromosome level in L. chinense internally. Because some contigs are still not assembled into chromosomes, we merged all remaining contigs that were not assembled according to the serial number. Each contig was segmented into a 50-bp N base. The final sequence was named chrF (chr-Fake). MCScanX was used to find out the collinearity and duplication relationship among genes in L. chinense, and the results are shown using Circos (Krzywinski et al., 2009). RIdeogram (Hao et al., 2020) was used to create a figure for the collinearity analysis carried out on four species.

The materials of somatic embryogenesis were collected sequentially, starting with collecting pre-embryonic mass (PEM) and then collecting materials from the small plant stage (PL); the materials were collected from a somatic embryogenesis system (Li et al., 2012). In detail, the embryogenic callus with pro-embryogenic cell mass (PEM) was cultured in the induced media after the preculture for 10 days (ES1). The cells after sieving were named ES2. We collected the samples that were cultured in the induced media at the first and third day named ES3 and ES4, respectively. Globular embryos could be observed approximately 7 days that we collected the globular embryos as the third sample named ES5. And then, we collected heart-shaped embryos (ES6), torpedo-shaped embryos (ES7), immature cotyledonary embryos (ES8), mature cotyledonary embryos (ES9), and little plant (PL) in proper order. Each sample has three independent biological replications. The plants of hybrid Liriodendron, which regenerate from somatic embryos, were cultivated for 2 months in 3/4 MS medium. All plants were cultivated in a room containing tissue culture under white light (light for 16 h; darkness for 8 h). The transient genetic transformation method used is the same as the method used in Tamarix hispida (Zang et al., 2017). We collected whole plants 48 h after transformation. Based on the limitation of the genetic transformation system, only the first mature leaves have high transgenic efficiency. Therefore, only apical leaves were collected after the transformation. Additionally, three independent biological replicates were prepared for later expression validation. To simulate cold stress, plants were held at 4°C for 12 h, 24 h, and 48 h. To simulate heat and drought stress, plantlets were held at 40°C or subjected to a 15% PEG 6000 treatment for 1 h, 3 h, 6 h, 12 h, 1 day, and 3 days. We collected leaves as samples with three independent biological replicates. All samples were frozen in liquid nitrogen immediately and then stored at −80°C for RNA extraction by RNA-seq analysis and quantitative real-time PCR (qRT-PCR) experiments.

We downloaded transcriptome data from NCBI (PRJNA679101 and PRJNA679089). We obtained the expression data of LcGRAS in treatments with cold (raw data unpublished), hot, and drought ( Supplementary Table S3 ) conditions and during somatic embryogenesis (raw data unpublished) ( Supplementary Table S4 ). The transcriptome data of different tissues are shown in Supplementary Table S5 . We used R package (pheatmap) to make the heatmap. Total RNA was extracted from each sample by using a FastPure Plant Total RNA Isolation Kit (Vazyme). HiScript 1st Strand cDNA Synthesis Kit (Vazyme) was used for first-strand cDNA synthesis. We used the AceQ qPCR SYBR Green Master Mix (Vazyme) to carry out qRT-PCR analysis in the LC480 platform. The primers used in this study are listed in Supplementary Table S4 . We performed three biological replicates for each gene. The gene named Tubulin was used as the reference gene. The relative expression levels were calculated using 2-ct.

We performed a genome-wide HMM search and homology analysis to obtain putative GRAS genes by using all protein sequences from the reference genome of Liriodendron. We kept the putative targets when the E-value was less than 1e-5. The de novo transcriptome of Liriodendron was used to confirm that no possible genes have been missed. We further confirmed that all candidates have the right domains or motifs and reannotated the genes by using the online FGENESH software and Conserved Domain Search Service (CD Search) in Softberry and NCBI, respectively. In total, 49 genes were considered as GRAS genes in the L. chinense with a typical GRAS domain and complete gene structure ( Supplementary Table S1 ). The length of the majority of proteins was 400–700 amino acids, but a few of them were longer than 800 amino acids. Protein Lc34435, the longest gene in all GRAS genes in L. chinense, is 995 amino acid-long.

The GRAS genes do not have a complex domain structure. All have a representative domain named GRAS (pfam03514) or GRAS superfamily (cl15987). Notably, the DELLA subfamily is named based on its extra DELLA (pfam12041) domain in the N-terminal ( Figure 1A ). We searched the conserved motifs of all GRAS proteins by using MEME via the classic mode. Almost all motifs were located on the GRAS domain or GRAS superfamily ( Figures 1A, B ). Overall, 10 motifs were present for most genes. However, based on the results of subfamily classification ( Figure 2 ), we found that some motifs help in the identification and functional characterization of GRAS subfamilies. For example, only Lc09089 missed Motif 1. Motif 1 may be an element for identifying SCL4/7. Motif 8 was missing in SCR, HAM-a, and HAM-b. HAM-a also missed Motif 5, and SCR also missed Motif 6 ( Figure 1B ).

It is worth noticing that only Lc34435 has a domain named the Retrotrans_gag superfamily (cl29674) ( Figure 1A ), which does not exist in any GRAS gene. However, Lc34435 has a standard and representative domain named the GRAS superfamily, which is unique to the GRAS genes. Moreover, Lc34435 is a highly homologous relative to other LcGRAS proteins. For these reasons, we kept this gene as a member of the LcGRAS family. Lc23231 is a structurally incomplete GRAS protein with a shorter gene length and an incomplete domain in the C terminus ( Figure 1 ). We found that there are similar structurally incomplete genes in the GRAS family in A. thaliana and other species as well. However, when aligning the protein sequences from these several species, we found that they do not share any motifs ( Supplementary Figure S1 ). Since this is a common phenomenon for the GRAS family in different species, we kept this gene for the following analysis.

LcGRAS was distributed on the chromosomes except for Chr6, Chr12, and Chr15 ( Figure 3 ). We found that some genes are highly identical in terms of motif and domain distribution, such as Lc29387 and Lc35375. They have the same length with respect to coding sequences and only a two-nucleotide difference. It might be derived from a duplication event. To find out all possible duplication genes, we used MCScanX to calculate the collinearity relationship in LcGRAS ( Figure 3 ). There are a total of 16 genes belonging to WGD (whole-genome duplications), which form eight pairs. In addition, three couples of tandem repeats and six genes of proximal repeats were found in Chr7, Chr9, and ChrF. Duplication may be a cause for the increased number of LcGRAS family members.

Cis-acting regulatory elements (CREs) are closely related to gene function regulation, and their divergence is a common cause of evolutionary change (Wittkopp and Kalay, 2012). To address what could be the upstream regulators, we used a 2,000-bp-length sequence upstream of each gene to search for CREs by using PlantCARE software. Except for common CREs such as A-box and TATA-box, CREs are related to regulated hormones, stress responsiveness, light responsiveness, development, and MYB ( Supplementary Figure S2 ). Light responsiveness and MYB-binding sites have the greatest number of all CREs. Except for the GT1 motif and G-box, there were 24 different CREs related to light responsiveness. This suggests that LcGRAS may be intensively regulated by the light signal. Furthermore, LcGRAS may be under hormone regulation. The CREs, such as AuxRR-core, ABRE, P-box, O2-site, TCA-element, as-1, and ERE, were related to auxin, ABA, GA, zein, salicylic acid, MeJA, and ethylene, respectively. W-box, WRE3, WUN motif, and STRE are CREs that are involved in wound or stress responsiveness. MBS is an MYB-binding site involved in drought inducibility, and LTR is involved in low-temperature responsiveness. This suggests that LcGRAS may have functions when plants are under stress, such as cold or drought conditions, and MYB may potentially be a regulator of LcGRAS.

To explore the phylogenetic relationship of GRAS and classify the subfamily of LcGRAS, we include the GRAS gene family from A. thaliana, A. trichopoda, O. sativa, P. abies, P. trichocarpa, and E. ferox to construct the phylogenetic tree. Firstly, we identified 67 GRAS genes in E. ferox by using the blastp of NCBI and HMM, which is the same as the work carried out on L. chinense ( Supplementary Figure S3 ). E. ferox, which belongs to magnoliids, is used to support the phylogenetic status of L. chinense.

In total, 385 genes were divided into 14 clusters ( Figure 2 ). One of 14 clusters, which contain Os_04g35250, At_00017.27, Pt_019G085600.1, and Pt_013G114900.1, was named as unclassified. This is because these four genes independently form a branch and cannot be properly clustered into any known clade. HAM genes were separated into three subclades, and this separation was based on whether they have Motif 5 and Motif 8 or not ( Figure 1B ). Notably, it is interesting that the members of cluster HAM-t all came from woody plants.

Next, based on subfamilies, we calculated Ka and Ks for each LcGRAS gene with genes in the same subfamily. The Ka/Ks value of gene pair Lc29387 and Lc35375 is expected to be 49.4, and the other gene pairs belong to neutral evolution or purifying selection ( Supplementary Figure S4 ). The resulting p-value of Lc29387 and Lc35375 is approximately 0.375. Since these two genes are extremely similar to each other and only have two SNPs (single nucleotide polymorphism), evidence that these two genes belong to positive selection is insufficient. In general, the gene pairs in Liriodendron have a lower non-synonymous substitution rate, suggesting that the genes in the same subfamily may not have strong functional divisions.

To further address the phylogenetic relationship with other species, we constructed four synteny maps containing L. chinense, A. thaliana, E. ferox, and P. trichocarpa. A total of 34 LcGRAS genes showed a syntenic relationship with these species ( Figure 4 ). It is obvious that LcGRAS had the most syntenic gene pair with respect to P. trichocarpa (63 pairs) and the least syntenic gene pair with respect to A. thaliana (17 pairs). We found that the unique syntenic gene pairs of L. chinense and P. trichocarpa have an independent clade of the phylogenetic tree, which only consisted of genes from woody plants, such as Lc00207a (HAM-t), Lc06445 (HAM-t), Lc14185 (SHR), and Lc14801 (SHR). The function of these genes is closely related to the growth, development, and morphogenesis of trees. Moreover, the genes, such as Lc01491 (HAM-a), Lc29299 (HAM-a), Lc18765 (HAM-a), Lc11923 (DLT), Lc02091 (PAT), Lc09089 (SCL4/7), Lc16825 (LAS), Lc19053 (LAS), Lc26517 (SHR), Lc12790 (LiSCL), and Lc29387 (LiSCL), had syntenic genes with GRAS genes in all three species, suggesting an important role of these genes in evolution, and they exhibit more basic functions in plants.

The GRAS gene family not only has diverse functions in plant growth but also is responsive to most biotic or abiotic stresses. Heat, drought, and cold are the common stresses for plants in nature. We analyzed the transcriptome data of Liriodendron under heat, drought, and cold stress treatments to study the expression pattern of LcGRAS. We found that the expression pattern of LcGRAS is variable ( Figure 5 ). Most genes of LcGRAS have changed expression levels within 24 h and 48 h of these three treatments. Some genes, such as Lc02461 and Lc34435, were barely expressed ( Supplementary Table S3 ). In heat stress, the genes in Group 1 and Group 4 are sensitive to heat stress, for which their expression levels could change within 1 h ( Figure 5A ). Under drought treatment, LcGRAS genes in Group 1, Group 2, and Group 3 were transiently upregulated, and the genes in Group 4 were gradually downregulated ( Figure 5B ). In the cold treatment, most LcGRAS genes respond after 12 h, and the expression patterns could be divided into four groups ( Figure 5C ). In Group 1, genes were downregulated and showed persistent low expression over 48 h. The genes in the second group were upregulated at 12 h and were immediately downregulated at subsequent time points. The two other groups were upregulated at 24 h or 48 h. With longer stress treatments (72 h), the expression levels show fluctuation. In contrast, all genes in the PAT subfamily were all upregulated after cold stress, suggesting that PAT genes may not only have a function in phytochrome A signal transduction but can also help plants resist stressful environments. In addition, we selected a subset of genes to verify that the expression level changes were consistent with the transcriptome. The relative expression levels of PAT3, PAT4, and PAT5 by qRT-PCR were similar with transcriptome data ( Supplementary Table S7 and Figure 6 ). With the published papers (Wu et al., 2021; Wu et al., 2022), the studies demonstrate the reliability of transcriptome data.

As mentioned above, PAT genes are upregulated after cold treatment. We reasoned that this subfamily may have functions in responding to cold stress. It has been reported that there is complex crosstalk between light and cold signaling pathways (Kidokoro et al., 2022). For example, phytochrome B could enhance cold resistance by promoting the degradation of PIF (phytochrome-interacting factor) proteins, which in turn release COR (cold-regulated) from the inhibition of PIF (Xu and Deng, 2020). Phytochrome A has also been reported to play a role in acquiring freezing tolerance in Arabidopsis (Prerostova et al., 2021). Here, we identified two PIF and four COR genes in L. chinense. Consistently, we noticed that one PIF (Lc20729) gene and one COR (Lc23651) gene were strongly regulated under cold stress in our transcriptome data ( Supplementary Table S7 ). Additionally, the result of qRT-PCR demonstrated that the expression levels of PIF and COR changed in accordance with transcriptome data under cold stress ( Supplementary Figure S5 ). In order to investigate the function of PAT in cold resistance, we transiently overexpressed six LcPATs (PAT1-PAT6) in Liriodendron based on the method developed in T. hispida (Zang et al., 2017) and quantified the expression of COR and PIF. First, we used 35S:GUS as the control to test whether the transient expression system works in Liriodendron, and we found that the leaves exhibit strong GUS staining. By contrast, we could not detect clear GUS staining in the root or stem ( Figure 7 ). Therefore, we only collected the transformed leaves to analyze the expression levels of COR and PIF via real-time quantitative PCR (qPCR). PAT1–PAT6 were overexpressed after transformation, although the overexpression level varied among biological replicates ( Supplementary Figure S5 ). The sample that was transformed with an empty vector was used as the control. We found that the COR was upregulated, and PIF was downregulated after overexpressing PATs. COR and PIF responded to the overexpression of PAT1 and PAT4 more actively, but not to PAT3 and PAT5. We could not rule out the possibility that PAT3 and PAT5 were lowly expressed after transformations. These results suggest that the cold may stimulate the expression of PAT via PIF and COR to help plants acquire cold tolerance.

It is reported that some GRAS genes are highly expressed during embryogenesis and in non-differentiated proliferating embryogenic tissue (Abarca et al., 2014). However, whether GRAS genes are expressed during somatic embryogenesis is not known yet. We analyzed the expression level of LcGRAS during somatic embryogenesis in Liriodendron hybrid. We found that not all LcGRASs are expressed at a high level during somatic embryogenesis. They could be divided into four groups ( Figure 8 ; Supplementary Table S4 ). In Group 1 and Group 2, the genes were not expressed until after a heart-shaped embryo developed, which is similar to GRAS genes in Cucumis sativus (Wiśniewska et al., 2013). The genes in Group 3 and Group 4 were expressed in PEMs or ES1. The majority of LcGRAS genes are not expressed in ES3 and ES4. Therefore, LcGRAS may not be involved in the early patterning stages of somatic embryogenesis, but they have a function in late embryo development.

Next, we checked the transcriptome data of different tissues in mature plants ( Figure 8B ; Supplementary Table S5 ). Most LcGRAS genes show a strongly specific expression pattern in different tissues. Two of the DELLAs were expressed in all eight tissues, and this indicated that DELLAs are closely related to organ development and produce a long-term effect on plant growth ( Supplementary Table S5 ). There are a few genes expressed in the phloem and xylem, but many LcGRASs are expressed in the bark, leaves, bud, and stigma. Interestingly, some LcGRASs have higher expression levels in the stigma than in the stamen, such as Lc11923, Lc01491, Lc18765, Lc02091, Lc09089, and Lc29007. These genes may have key roles in the formation of floral organs.

In L. chinense, we identified 49 genes that belong to the GRAS gene family. The gene number of this gene family is large. There could be more than a hundred members in G. hirsutum (Zhang et al., 2018), or there could only be 34 members in A. thaliana. The core function of GRAS is closely related to plant development and stress resistance.

Among 49 LcGRASs, we found that Lc23231 is a short gene with a GRAS domain. The short genes were also identified in A. thaliana, A. trichopoda, O. sativa, and other several species. The alignment result shows that there is no obvious motif that could be found among these short genes ( Supplementary Figure S1 ). We found that Lc23231 has a unique absence of Motif 3 ( Figure 1 ), which leads to the sequence terminating prematurely. Similarly, AT5G67411 (SCL16) is believed to be a putative expressed pseudogene and a unique member of the GRAS family, lacking signature motifs in A. thaliana (Lee et al., 2008). However, Lc23231 still responds to stress treatments and has different expression levels between different tissues in Liriodendron. VviSCL3B also presents a premature stop codon and lacks motifs PFYRE and SAW, but it was found to be expressed (Grimplet et al., 2016). Therefore, we speculate that Lc23231 may not be a pseudogene and has functions, which may be due to replication events or structure variation. The delay of Lc34435 termination was also unexpected. We hypothesized that this was related to retrotransposon insertion. Since Lc23231 and Lc34435 are barely expressed in our existing transcripts, determining whether changes in the nucleotide sequence affect their gene function is difficult. This problem requires performing more molecular experiments for proof. In addition, the majority of LcGRAS genes (76.3%) were intron-less. This is similar to most results from analyses of other plants, such as C. annuum (84%) (Liu et al., 2018), Hippophae rhamnoides (58.1%) (Yu et al., 2021), and Hibiscus hamabo (62.71%) (Ni et al., 2022). One possible reason is that GRAS genes originate from the bacterial prokaryotic genome (Zhang et al., 2012).

Tandem and segmental duplications are the main ways that plants can expand the gene family (Cannon et al., 2004). A large number of duplications was identified in G. hirsutum (Li et al., 2018). However, WGD is the main form of duplication in Liriodendron, which contributed significantly to gene family expansion for Liriodendron.

In general, 49 genes of LcGRAS were divided into 13 subgroups, including LiSCL, SHR, PAT, SCL3, RAM, DELLA, DLT, SCR, SCL4/7, LAS, HAM-a, HAM-t, and HAM-b. We distinguished the three subclusters from the cluster of HAM. The members of cluster HAM-t are from woody plants. The following question still requires more data and experiments in order to obtain answers: Is this cluster unique for woody plants or do the genes of HAM-t have a function that is unique for woody plants and does this indicate the evolution traits of woody and herbaceous plants? Collinearity analyses with L. chinense, A. thaliana, P. trichocarpa, and E. ferox show that the number and homology of GRAS are species-specific, other than evolutionary status. L. chinense has more collinearity blocks with P. trichocarpa than E. ferox. The evolution of the GRAS gene family was hypothesized to occur after the divergence event (Grosche et al., 2018), and GRAS proteins may play significant roles in complex morphogenetic pattern formation in higher plants such as rootless, leafless thallophytes (Laskar et al., 2021).

We detected the cis-acting regulatory elements in the promotors of LcGRAS. Except for gibberellin responsiveness and light responsiveness, we also found cis-acting regulatory elements related to auxin, MeJA, ABA, and stress responsiveness, suggesting that GRAS has complex biological functions. GRAS not only may be involved in signaling pathways for GAs (Peng et al., 1999), phytochrome (Tong and Chu, 2009), and brassinosteroids (BRs) (Silverstone et al., 1998) but also may be involved in other plant hormones such as auxin, MeJA (Hou et al., 2010), and ABA. Their functions of responding to stress and growth may be inextricably linked with the regulator of ABA or auxin. GRAS genes play significant roles in plant growth, production, and defense against abiotic and biotic stresses (Khan et al., 2021). The development stage, tissue, environmental conditions, and genotype affect expression levels (Hirsch and Oldroyd, 2009), and this is consistent with the results of our qRT-PCR and transcriptome data. Moreover, it has been reported that GRAS is involved in maintaining the normal development of seeds (Chen et al., 2021). The transcriptome data show that LcGRAS is expressed during somatic embryogenesis from the globular embryo stage. This suggests that LcGRAS may have similar functions in maintaining embryo development during somatic embryogenesis in Liriodendron. Furthermore, nine genes (Lc01789, Lc02091, Lc07168, Lc16536, Lc17878, Lc22978, Lc31454, Lc33380, and Lc35375) were highly expressed and regulated during somatic embryogenesis and under abiotic stress. These genes may respond to stress and participate in development control at the same time. We found that the expression level of PAT is upregulated after cold and drought treatment ( Figure 8 ). There was a similar discovery with respect to Vitis vinifera (Grimplet et al., 2016), which indicated that PAT may have a function during cold and drought stress. We used a quick method to confirm its expression levels in a cold-regulated network. The results show that PAT interacts with the cold regulatory network, but the specific regulatory mode is still unclear. The method of transient genetic transformation is not stable in Liriodendron, rendering the results unrepeatable. Stable transformations may be a better method for exploring the gene function of PATs under cold stress in Liriodendron.

Overall, differences in the number of GRAS gene family members among different species may lead to differences in response to stress or hormones. It is worth further studying the diverse functions of the members of the GRAS family.