Genome-Wide Identiﬁcation and Characterization of G2-Like Transcription Factor Genes in Moso Bamboo ( Phyllostachys edulis )

: G2-like (GLK) transcription factors contribute signiﬁcantly and extensively in regulating chloroplast growth and development in plants. This study investigated the genome-wide identiﬁcation, phylogenetic relationships, conserved motifs, promoter cis-elements, MCScanX, divergence times, and expression proﬁle analysis of PeGLK genes in moso bamboo ( Phyllostachys edulis ). Overall, 78 putative PeGLKs ( PeGLK1 – PeGLK78 ) were identiﬁed and divided into 13 distinct subfamilies. Each subfamily contains members displaying similar gene structure and motif composition. By synteny analysis, 42 orthologous pairs and highly conserved microsynteny between regions of GLK genes across moso bamboo and maize were found. Furthermore, an analysis of the divergence times indicated that PeGLK genes had a duplication event around 15 million years ago (MYA) and a divergence happened around 38 MYA between PeGLK and ZmGLK . Tissue-speciﬁc expression analysis showed that PeGLK genes presented distinct expression proﬁles in various tissues, and many members were highly expressed in leaves. Additionally, several PeGLKs were signiﬁcantly up-regulated under cold stress, osmotic stress, and MeJA and GA treatment, implying that they have a likelihood of affecting abiotic stress and phytohormone responses in plants. The results of this study provide a comprehensive understanding of the moso bamboo GLK gene family, as well as elucidating the potential functional characterization of PeGLK genes.


Introduction
Chloroplasts contain the green pigment chlorophyll, which can carry out photosynthesis and is thus essential in plant life processes [1][2][3]. Recent research has supported the idea that chloroplasts are derived from protoendosymbiosis events associated with these cyanobacteria [4][5][6]. Furthermore, photosynthetic organs are assembled under the coordinated regulation of chloroplasts and nucleus. Plastids function not only in photosynthesis and nutrient storage, but also in the compartmentalization of various metabolic intermediates [7]; for example, most amino acids, all fatty acids, purines and pyrimidine bases, terpenes and various pigments, and hormones are synthesized in plastids. In addition, proplastids in plant meristem cells (leaf sheaths in dark cotyledons) are converted to chloroplasts of mesophyll cells under light, while GLK genes participate in the regulation and control of chloroplast formation in transition and early maturity and are indispensable in the development of angiosperm chloroplasts [2,[7][8][9][10].
GLK genes encode members of the GARP superfamily of nuclear transcription factors [11] defined by Golden2 in maize (Zea mays L.), Response regulator-B (ARR-B) proteins in Arabidopsis thaliana [12], and phosphate starvation response1 (PSR1) protein in Chlamydomonas. Most GLK proteins contain two conserved domains: an Myb-DNA binding domain (DBD, containing an HLH motif) and a C-terminal domain (containing a conserved GCT box) [13,14].
GLK genes promote the production of plant chloroplasts, and optimise photosynthesis under different biotic and abiotic environmental stress conditions [15,16]. For instance, At-GLK1 and AtGLK2 function in the formation and development of chloroplasts redundantly in Arabidopsis [17,18], and AtGLK1 OE increases resistance to Fusarium graminearum, a wide range of host pathogens that cause significant losses in cereal crops. In contrast, the glk1 glk2 double mutant provides resistance to Hyaloperonospora arabidopsidis (Hpa), implying a potential role of GLKs in plant defense and disease resistance [19]. In addition, GLK1 and GLK2 participate in chloroplast growth and development of C3 photosynthetic plants redundantly [15]. In tomato (Solanum lycopersicum), SlGLK2 OE enhances fruit photosynthesis and chloroplast development gene expression, resulting in elevated carbohydrates and carotenoids in ripe fruit, which can lead to strong resistance to cold, drought, heat and other abiotic stresses [20]. In maize, the spatial compartmentalization of ZmGLK1 transcripts and G2 might function as a specialization that was essential to the chloroplast development of mesophyll cells and distinct bundle sheath in plant tissues [14,17].
Moso bamboo (Phyllostachys edulis) is a rapidly growing forest product that is widely used to manufacture paper, art ware, and food resources, with multiple economic, ecological, and cultural values in China and many other countries, sees a possibility for gene identification and function analysis that will help with it contributing on an even wider spectrum with the completion of the genome sequence [6,20,21]. The GLK gene family, previously identified and described within certain plants, for example, maize [22], Arabidopsis [6], tomato [23], and tobacco [24]. However, there are scarcely any reports on the characteristics and functions of GLKs in Phyllostachys edulis. Thus, this study identified 78 putative PeGLK genes and conducted a systematic analysis, including phylogenetic relationships, gene and protein structures, conserved motifs, promoter cis-elements, chromosomal localization, MCScanX, evolutionary patterns, and tissue-specific expression. The expression patterns of PeGLKs in response to abiotic stress (cold and osmotic) and phytohormone (MeJA and GA) treatment were analyzed as well. These results are valuable insightful for further understanding of potential roles of the moso bamboo GLK gene family.

Plant Material and Growth Conditions
For this experiment, moso bamboo (Phyllostachys edulis) seeds, which were provided by Gongcheng Yao Autonomous County, Guangxi Zhuang Autonomous Region, China, were germinated and grown under daylight conditions of 16 h light/8 h dark and maintained at 28 • C and 80% relative humidity in an artificial growth chamber for 100-day-old. The 100day-old plants were subjected to cold stress, osmotic stress, and MeJA and GA treatments. Specifically, Phyllostachys edulis plants were placed in a 4 • C growth chamber and leaves were gathered under cold stress at 0, 1, 3, 6, 12, and 24 h, and under osmotic stress at 0, 1, 2, 3, 6, and 12 h after treatment with pure water containing 20% polyethylene glycol (PEG) 6000. For hormone treatments, a solution of 100 µM Methyl jasmonic (MeJA) and 300 mg/L gibberellic acid (GA) was sprayed onto seedlings based on the requirements, and seedlings were sampled randomly at 0, 1, 3, 6, and 12 h. Negative control plants were treated in a 28 • C artificial growth chamber and sprayed with distilled water.

Identification of PeGLK Genes
Phyllostachys edulis protein and nucleic acid sequences were acquired from the moso bamboo genome database (http://parrot.genomics.cn, accessed on 10 March 2022 members in moso bamboo (Phyllostachys edulis) by analyzing conserved domains of these proteins. Detailed information of the PeGLK genes, including gene IDs, physical positions, gene and protein sequences, and coding sequences (CDSs), were retrieved from the bamboo database. The ExPASy online tool (http://www.expasy.ch/tools/pi_tool.html, accessed on 10 May 2022) was used to determine the physical parameters of PeGLK genes, including molecular weight (MW), isoelectric point (pI), and open reading frame length.

PeGLK Sequence Alignment and Phylogenetic Analysis
To identify the domain organization of 78 PeGLK proteins, multiple alignments of PeGLK domain-containing sequences were performed using ClustalW [25], and a phylogenetic tree based on the complete PeGLK sequences was constructed using the neighborjoining method by MEGA 7.0 with 1000 bootstrap replications (https://megasoftware.net/, accessed on 10 April 2022). In addition, the combined phylogenetic tree of PeGLK and ZmGLK protein sequences was also generated with the N-J method.

Gene and Protein Structure Analysis
The exon-intron structures were obtained by contrasting the CDSs and relevant genomic DNA sequences of PeGLK genes from GSDS (http://gsds.cbi.pku.edu.cn/, accessed on 14 April 2022). Conserved motifs present in the PeGLK proteins were examined by the online MEME tool (http://meme-suite.org/tools/meme, accessed on 18 April 2022) [26]. The optimum width of motifs was limited to 5-50, and the maximal number of motifs was set to 10 residues. Then, we used Pfam tools to determine the motif annotations. The protein homology models were predicted by aligning GLK protein sequences with the HMM-HMM search in intensive mode on the Phyre2 website (http://www.sbg.bio.ic.ac.uk/phyre2/html/page.cgi?id=index, accessed on 26 April 2022) [27].
To perform a collinearity analysis of PeGLK protein in Phyllostachys edulis, BLAST was used to compare the entire sequence, with a cutoff of truncated E-value of 1 × 10 −20 [29]. The collinearity analysis of moso bamboo and maize also used the same truncated E-value. MCScanX software was used to produce collinearity blocks for the entire genome by calculating the BLASTP results [28].

Identification of PeGLK Genes in Phyllostachys edulis
To identify the moso bamboo GLK gene family members, previously reported Arabidopsis AtGLK1 (AT2G20570) and AtGLK2 (AT5G44190) protein sequences [17] were used as BLASTP queries to search for available protein databases in the moso bamboo genome bank. We identified 78 PeGLK candidate genes, most of which were used to confirm the presence of the Myb DNA-binding domain (DBD) and C-terminal domain (containing a conserved GCT box) using the Pfam database. Except for the conserved GLK Myb DNA-binding domain, the members of each subfamily have suspected functional diversity, represented by unique motifs. For instance, among groups 1-7, each has an Myb-CC-LHEQLE domain, which is a type of Myb-like domain. Group 10 has an REC-typeB-ARR-like domain ( Figure 1). To discern the similarities between moso bamboo GLK domains, the domain sequences of 78 PeGLK proteins were blasted with the DNAMAN 8 platform ( Figure S1). The results showed that PeGLKs were conserved in two regions of a putative Myb DNA-binding domain, with the HLH structure of the first helix containing the initial sequence, PELHRR, and the second helix containing NI/VASHLQ, which was consistent with the GLK members in Arabidopsis, maize [22], tobacco, and tomato [40].
Basic information about PeGLK genes, such as the accession number, position, and physicochemical parameters, is presented in Table 1. The molecular weight (MW) varies from 18.08 to 75.00 kDa and the length of CDSs ranged from 486 to 2079 bp, and encoded sequences ranged from 161 to 692 aa. In addition, the theoretical isoelectric point (pI) ranged from 5.51 to 10.43.

Phylogenetic and Exon-Intron Structure Analysis of GLK Genes
From a phylogenetic tree formulated on the aligning of 137 sequences of GLK proteins from moso bamboo (78) and maize (59), the phylogenetic relationships of GLK proteins among different species were derived. The characteristics of ZmGLK genes from maize are listed in detail in Table S1. In this phylogenetic tree, GLK family members were classified into 16 groups based on evolutionary relationships and motif analysis, and poplar GLK family members were allocated into 13 groups (G1-G13, but not G14-G16). The numbers of PeGLKs in different groups were uneven: groups 1 to 13 contained 9, 4, 7, 5, 3, 4, 2, 15, 14, 7, 3, 6, and 2 genes, respectively ( Figure 2). To further explore the evolutionary relationships of PeGLK genes, a single phylogenetic tree was constructed using the moso bamboo PeGLK protein sequences. PeGLK proteins were also classified into 13 distinct subfamilies, in accordance with the phylogenetic tree constructed by moso bamboo and maize. The exon-intron predictions were identified by the online Gene Structure Display Server (GSDS) tool. As shown in Figure 3, the number of exons in different subfamilies varied from 1 to 12, and most PeGLK genes (91%) had five or more exons. In addition, most genomic sequences contained upstream and downstream sequences, except for PeGLK5, PeGLK13, PeGLK48, PeGLK50, PeGLK67, PeGLK68, PeGLK69, and PeGLK77, whereas PeGLK12 and PeGLK41 had only upstream sequences, and PeGLK19, PeGLK29, PeGLK36, PeGLK51, PeGLK70, and PeGLK74 had only downstream sequences. In general, the vast majority of putative paralogous had significant changes in structural organization and intron/exon numbers, except for three paralogous (PeGLK5/PeGLK6, PeGLK62/PeGLK65, PeGLK48/PeGLK67), which had the same intron/exon numbers but different intron lengths. To further explore the evolutionary relationships of PeGLK genes, a single phylogenetic tree was constructed using the moso bamboo PeGLK protein sequences. PeGLK proteins were also classified into 13 distinct subfamilies, in accordance with the phylogenetic tree constructed by moso bamboo and maize. The exon-intron predictions were identified by the online Gene Structure Display Server (GSDS) tool. As shown in Figure 3, the number of exons in different subfamilies varied from 1 to 12, and most PeGLK genes (91%) had five or more exons. In addition, most genomic sequences contained upstream and downstream sequences, except for PeGLK5, PeGLK13, PeGLK48, PeGLK50, PeGLK67, PeGLK68, PeGLK69, and PeGLK77, whereas PeGLK12 and PeGLK41 had only upstream sequences, and PeGLK19, PeGLK29, PeGLK36, PeGLK51, PeGLK70, and PeGLK74 had only downstream sequences. In general, the vast majority of putative paralogous had significant changes in structural organization and intron/exon numbers, except for three paralogous (PeGLK5/PeGLK6, PeGLK62/PeGLK65, PeGLK48/PeGLK67), which had the same intron/exon numbers but different intron lengths.

Identification of Conserved Sequence Motifs and Homology Modeling in PeGLK Genes
To evaluate the phylogenetic kinship, we examined the conserved motifs of 78 PeGLK proteins using the MEME software ( Figure 4). In Table S2 the eight distinct motifs found are shown, with their specific characteristics noted. Functional annotation was made for six putative motifs with the Conserved Domain Database (CDD), in which motifs 1 and 2 were defined as Myb-SHAQKYF, motif 3 was Myb-CC-LHEQLE, and motifs 5, 6, and 8 were REC-type B-ARR-like. However, the remaining two putative motifs did not have functional annotations. Protein family members within the same group shared similar or identical components and spatial distributions, which implies that these proteins have similar functions. For instance, all of the PeGLK proteins were characterized by motif 1 in the Myb DNA-binding domain, with an HLH structure (the first helix contains initial sequence PELHRR and invariably comprises 14 amino acids, and the second helix contains an initial NI/VASHLQ motif. These helices are separated by a 22 amino acid

Identification of Conserved Sequence Motifs and Homology Modeling in PeGLK Genes
To evaluate the phylogenetic kinship, we examined the conserved motifs of 78 PeGLK proteins using the MEME software ( Figure 4). In Table S2 the eight distinct motifs found are shown, with their specific characteristics noted. Functional annotation was made for six putative motifs with the Conserved Domain Database (CDD), in which motifs 1 and 2 were defined as Myb-SHAQKYF, motif 3 was Myb-CC-LHEQLE, and motifs 5, 6, and 8 were REC-type B-ARR-like. However, the remaining two putative motifs did not have functional annotations. Protein family members within the same group shared similar or identical components and spatial distributions, which implies that these proteins have similar functions. For instance, all of the PeGLK proteins were characterized by motif 1 in the Myb DNA-binding domain, with an HLH structure (the first helix contains initial sequence PELHRR and invariably comprises 14 amino acids, and the second helix contains an initial NI/VASHLQ motif. These helices are separated by a 22 amino acid loop.). Additionally, there was specificity within different subfamilies. For example, motif  The structure and homology modeling of PeGLK proteins were determined by Phyre2, and alignment of the protein sequences with HMM-HMM search was conducted in intensive mode [41]. As a result, all 78 PeGLK proteins could be modeled with certainty. The structure and homology modeling of PeGLK proteins were determined by Phyre2, and alignment of the protein sequences with HMM-HMM search was conducted in intensive mode [41]. As a result, all 78 PeGLK proteins could be modeled with certainty. As shown in Figure 5, 10 PeGLKs (PeGLK5, PeGLK12, PeGLK13, PeGLK29, PeGLK49, PeGLK51, PeGLK53, PeGLK70, PeGLK71, and PeGLK74) had 100% of their predicted length modeled with >98% confidence.

Physical Locations and Duplication Events of PeGLK Genes in Phyllostachys edulis
The 78 PeGLKs were unevenly distributed on 23 moso bamboo chromosomes (Chr2-Chr24), except Chr1 ( Figure 6). Chromosomes 6 and 8 contained eight PeGLK genes, the highest number; chromosomes 15, 17, and 21 had seven, and 20 had five, whereas there was only one PeGLK gene on chromosomes 2, 5, 12, 14, and 18. Interestingly, Chr13 was the longest chromosome with only four PeGLK genes. These results indicate that the chromosome length was not proportional to the number of genes. To explore the potential mechanism of the PeGLK gene family, the MCScanX program was used to investigate potential segmental and tandem duplication events. As shown in Figure 7A, 31 segmental duplicated pairs and one tandem duplicated pair of PeGLK genes were identified in a synteny map ( Table 2). The 31 segmental duplicated pairs presented a biased distribution pattern, and no pairs were distributed on chromosome 1. These results suggest that moso bamboo PeGLKs may have arisen mainly from segmental duplication events.

Physical Locations and Duplication Events of Pelf Genes in Phyllostachys edulis
The 78 PeGLKs were unevenly distributed on 23 moso bamboo chromosomes (Chr2-Chr24), except Chr1 ( Figure 6). Chromosomes 6 and 8 contained eight PeGLK genes, the highest number; chromosomes 15, 17, and 21 had seven, and 20 had five, whereas there was only one PeGLK gene on chromosomes 2, 5, 12, 14, and 18. Interestingly, Chr13 was the longest chromosome with only four PeGLK genes. These results indicate that the chromosome length was not proportional to the number of genes. To explore the potential mechanism of the PeGLK gene family, the MCScanX program was used to investigate potential segmental and tandem duplication events. As shown in Figure 7A, 31 segmental duplicated pairs and one tandem duplicated pair of PeGLK genes were identified in a synteny map ( Table 2). The 31 segmental duplicated pairs presented a biased distribution pattern, and no pairs were distributed on chromosome 1. These results suggest that moso bamboo PeGLKs may have arisen mainly from segmental duplication events.      To get additional insight into PeGLK's evolutionary orthologous relationships, a synteny map was plotted between moso bamboo and maize. As shown in Figure 7B, 42 orthologous pairs (Pe-Zm) of moso bamboo and maize were determined (Table 3). Across the two, highly conserved microsynteny was observed in the regions of PeGLK genes, especially in Pe3 and Zm5, Pe7 and Zm6, Pe15 and Zm1, Pe17 and Zm5, and Pe21 and Zm1, all with three synteny genes.

Evolutionary Patterns and Divergence Times of GLK Genes in Moso Bamboo and Maize
To explore the selective constrains for the duplicated PeGLK gene pairs, a comprehensive analysis of the Ka/Ks ratio and Ks values was conducted using full-length sequences to estimate divergence times. As shown in Figure 8 and Table 2, the distribution of calculated Ks values of paralogous pairs (Pe-Pe) showed an average value of~0.2, indicating that PeGLK genes underwent a large-scale duplication event circa 15 MYA. A previous study estimated that the timing for whole-genome duplication of moso bamboo was 7-12 MYA [42], implying that large-scale duplication of PeGLK genes occurred earlier [33]. In addition, for the orthologous maize-moso bamboo pairs, the averaged Ks value was~0.5 ( Figure 8 and Table 3), indicating GLK genes differentiated 38 MYA. A comparison with a previous study revealed 42-52 divergence times between moso bamboo and maize, implying that GLK genes encountered gene evolution before the isolation of maize. In principle, Ka/Ks ratios greater than, equal to, and less than 1 indicate accelerated evolution with positive, neutral, and negative or stabilizing selection, respectively [26,43]. The Ka/Ks ratios of the Pe-Pe (Table 2) and Pe-Zm (Table 3) genomes were less than 1, which indicates that GLK genes experienced highly positive purifying selection among moso bamboo-maize and moso bamboo. 022, 27, x FOR PEER REVIEW 16 of 26

PeGLK Expression Levels in Different Tissues
Tissue-specific gene expression profiles provide critical clues for exploring gene function. In order to characterize the expression levels of PeGLKs, we analyzed the transcription levels of 36 PeGLK genes selected as representatives of each subfamily in four different tissues (leaf, stem, rhizome, and root) of moso bamboo (Table S3). As shown in Figure 9, the vast majority of PeGLK genes (except for PeGLK6, PeGLK7, PeGLK52, PeGLK61, PeGLK63, PeGLK68, PeGLK71) in the leaf, nine genes (PeGLK47, PeGLK51, PeGLK52, PeGLK60, PeGLK63, PeGLK65, PeGLK66, PeGLK69, PeGLK71) in stem, four genes (PeGLK44, PeGLK52, PeGLK53, PeGLK68) in rhizome, and 12 genes (PeGLK6, PeGLK7, PeGLK20, PeGLK21, PeGLK28, PeGLK36, PeGLK37, PeGLK39, PeGLK53, PeGLK61, PeGLK72, PeGLK78) in root presented high expression levels. Interestingly, the expression level of most PeGLKs was significantly higher in leaf than in the other tissues. We also defined that some genes showed a tissue-specific expression pattern. For example, PeGLK20, PeGLK21, PeGLK28, PeGLK36, PeGLK37, PeGLK39, PeGLK72, and PeGLK78 displayed high expression levels in leaves and roots, but low levels in stems and rhizomes. PeGLK52 exhibited high expression levels in stems and rhizomes, but low levels in leaves and roots. PeGLK47, PeGLK51, PeGLK60, PeGLK65, PeGLK66, and PeGLK67 presented high expression levels in leaves and stems, but low levels in rhizomes and roots. PeGLK44 and PeGLK70 displayed high expression levels in leaves and rhizomes, but low levels in stems and roots. In addition, only PeGLK63 and PeGLK67 were prominently expressed in stems, and only PeGLK68 was prominently expressed in rhizomes, whereas PeGLK6, PeGLK7, PeGLK21, PeGLK36, PeGLK61, and PeGLK78 were relatively higher in roots, PeGLK63 and PeGLK71 in stems, and PeGLK52 and PeGLK68 in rhizomes than in the other tissues. Among our previously identified paralogous genes, five pairs (PeGLK60 and PeGLK63, PeGLK62 and PeGLK65, PeGLK67 and PeGLK68, PeGLK70 and PeGLK72) shared the same or similar expression patterns in the four tissues, revealing the same evolutionary fate of

PeGLK Expression Levels in Different Tissues
Tissue-specific gene expression profiles provide critical clues for exploring gene function. In order to characterize the expression levels of PeGLKs, we analyzed the transcription levels of 36 PeGLK genes selected as representatives of each subfamily in four different tissues (leaf, stem, rhizome, and root) of moso bamboo (Table S3). As shown in Figure 9, the vast majority of PeGLK genes (except for PeGLK6, PeGLK7, PeGLK52, PeGLK61, PeGLK63, PeGLK68, PeGLK71) in the leaf, nine genes (PeGLK47, PeGLK51, PeGLK52, PeGLK60, PeGLK63, PeGLK65, PeGLK66, PeGLK69, PeGLK71) in stem, four genes (PeGLK44, PeGLK52, PeGLK53, PeGLK68) in rhizome, and 12 genes (PeGLK6, PeGLK7, PeGLK20, PeGLK21, PeGLK28, PeGLK36, PeGLK37, PeGLK39, PeGLK53, PeGLK61, PeGLK72, PeGLK78) in root presented high expression levels. Interestingly, the expression level of most PeGLKs was significantly higher in leaf than in the other tissues. We also defined that some genes showed a tissuespecific expression pattern. For example, PeGLK20, PeGLK21, PeGLK28, PeGLK36, PeGLK37, PeGLK39, PeGLK72, and PeGLK78 displayed high expression levels in leaves and roots, but low levels in stems and rhizomes. PeGLK52 exhibited high expression levels in stems and rhizomes, but low levels in leaves and roots. PeGLK47, PeGLK51, PeGLK60, PeGLK65, PeGLK66, and PeGLK67 presented high expression levels in leaves and stems, but low levels in rhizomes and roots. PeGLK44 and PeGLK70 displayed high expression levels in leaves and rhizomes, but low levels in stems and roots. In addition, only PeGLK63 and PeGLK67 were prominently expressed in stems, and only PeGLK68 was prominently expressed in rhizomes, whereas PeGLK6, PeGLK7, PeGLK21, PeGLK36, PeGLK61, and PeGLK78 were relatively higher in roots, PeGLK63 and PeGLK71 in stems, and PeGLK52 and PeGLK68 in rhizomes than in the other tissues. Among our previously identified paralogous genes, five pairs (PeGLK60 and PeGLK63, PeGLK62 and PeGLK65, PeGLK67 and PeGLK68, PeGLK70 and PeGLK72) shared the same or similar expression patterns in the four tissues, revealing the same evolutionary fate of duplicated genes. Above all, PeGLKs displayed diverse expression profiles in different tissues, which provides insight into the role of PeGLKs in multiple growth and development of Phyllostachys edulis.

Expression Profiles of GLK Genes under Abiotic Stress and Phytohormone Treatment
To determine the cis-regulatory elements in the putative promoter regions, we used PlantCARE to detect the 2000 bp upstream sequences of PeGLK genes [44]. Results portray that many cis-regulatory elements corresponding to LTRE (cold-responsive element), DRE (drought-responsive element), MeJA, and GA were detected, indicating that PeGLKs participate in plant development and stress responses ( Figure 10).

Expression Profiles of GLK Genes under Abiotic Stress and Phytohormone Treatment
To determine the cis-regulatory elements in the putative promoter regions, we used PlantCARE to detect the 2000 bp upstream sequences of PeGLK genes [44]. Results portray that many cis-regulatory elements corresponding to LTRE (cold-responsive element), DRE (drought-responsive element), MeJA, and GA were detected, indicating that PeGLKs participate in plant development and stress responses ( Figure 10). Several GLK genes have been reported to respond to cold and osmotic stresses in maize [22], tobacco [45], and tomato [24]. To investigate whether the expression of PeGLK Several GLK genes have been reported to respond to cold and osmotic stresses in maize [22], tobacco [45], and tomato [24]. To investigate whether the expression of PeGLK genes was affected by abiotic stress and phytohormone treatment, we detected the dynamic expression levels of 13 genes (PeGLK1, PeGLK7, PeGLK21, PeGLK36, PeGLK40, PeGLK48, PeGLK50, PeGLK53, PeGLK60, PeGLK61, PeGLK67, PeGLK70, and PeGLK72) as representatives of each subfamily under cold and osmotic stress and MeJA and GA treatment (Figures 11 and 12, and Table S5) (Figures 11 and 12, and Table S5). Figure 11. Expression patterns of 13 representative PeGLK genes in response to abiotic stress treatments. X-axis and Y-axis indicate time points after cold and osmotic treatments, and relative expression levels are standardized to reference gene TIP41. In the cold stress treatment, eight genes (PeGLK1, PeGLK21, PeGLK36, PeGLK40, PeGLK53, PeGLK60, PeGLK67, and PeGLK72) were obviously up-regulated compared with 0 h. For example, the expression of two genes (PeGLK40 and PeGLK60) peaked at 12 h, and PeGLK60 was the most highly expressed (>80-fold) after 12 h of treatment. Four genes (PeGLK7, PeGLK36, PeGLK53, and PeGLK72) gradually increased over time and peaked at 24 h. Seven genes (PeGLK1, PeGLK21, PeGLK48, PeGLK50, PeGLK61, PeGLK67, and PeGLK70) showed slight (<5-fold) changes in response to cold stress treatment.
In the osmotic stress treatment, five genes (PeGLK7, PeGLK36, PeGLK53, PeGLK67, and PeGLK70) were obviously up-regulated compared with 0 h. For instance, two genes (PeGLK53 and PeGLK70) gradually accumulated and peaked at 6 h, especially PeGLK70, which showed an expression level more than 50-fold higher at 6 h, whereas the expression levels of PeGLK7, PeGLK36, and PeGLK67 presented parabolic trends and were highest at 12 h.
In the MeJA treatment, four genes (PeGLK1, PeGLK36, PeGLK53, and PeGLK70) showed the highest expression at 12 h, and then a gradually increasing trend. The expression levels of four genes (PeGLK7, PeGLK40, PeGLK50, and PeGLK61) were highest at 6 h, especially PeGLK7, which presented more than 60-fold higher expression at 6 h. In the GA treatment, two genes (PeGLK50 and PeGLK53) were up-regulated, and PeGLK53 exhibited an expression level more than 35-fold higher at 12 h.
In total, 46.2, 38.5, 61.5, and 15.4% of genes, including six cold stress-related genes, five osmotic stress-related genes, eight MeJA-related genes, and two GA-related genes, were differently expressed; the detailed gene lists are shown in Figure S2. However, only the expression of PeGLK53 showed significant changes in response to the two abiotic stresses ( Figure 11) and two phytohormone treatments ( Figure 12). Additionally, the above results coincide with the promoter analysis of PeGLK members, indicating the widespread presence of several osmotic, cold, MeJA, and GA responsive cis-elements. Overall, these results imply that some PeGLK genes have potential functions in response to abiotic stress and phytohormone treatments.

Discussion
The function of GLK genes was first analyzed in maize and GLK genes were served as transcription factors. The GLK gene sequences were only found in photosynthetic eukaryotes such as green algae and higher plants, not present in the genome of cyanobacteria, which indicated that GLK genes were involved in the formation and development of chloroplast [45]. Members of the GLK gene family were characteristic groups in nuclear coding genes, and general feature of organisms with GLK genes were associated with chloroplasts and chlorophyll. Chloroplasts convert light energy into chemical energy and play a significant role in plant growth and development; thus, the study of GLK genes has become an important prospect of investigation [46].

GLKs in Phyllostachys edulis
The GLK transcription factor, which is a member of the newly classified GARP superfamily [14], plays an essential role in the formation and development of chloroplasts [11,47]. As reported in previous studies, detailed features and functions of GLK genes have been uncovered in Arabidopsis [6], maize [22], tobacco [24], and tomato [23]. However, little was known about the members of the moso bamboo GLK family until now. Here, 78 putative PeGLK genes were identified, which is 22 more than maize (59), and there are 30 and 42 times as many genes as tomato (54) and sorghum (45), respectively. The higher number of PeGLK genes was in accord with the genome duplication event that occurred in Phyllostachys edulis [33,48]. According to the phylogenetic analysis of grouping with maize sequences, the predicted moso bamboo GLK family members were divided into 13 groups, G1-G13. However, there were no orthologous genes in maize G14 and G15, which suggests divergence between maize and moso bamboo. In addition, the PeGLKs belonging to the same subfamilies displayed high similarity based on their domain structures and number of exons/introns, and these identical results indicate that the PeGLK groupings are relatively reliable.

Physical Locations and Divergence of GLKs in Moso Bamboo and Maize
According to the chromosome location analysis, PeGLKs were extensively and unevenly spread across 31 chromosomes of Phyllostachys edulis, which may be because of insertions, deletions, duplications, and reversions. Previous reports have shown that gene duplication is an important means of gene expansion in the process of plant genome evolution and is also essential for organisms to adapt to various environments during development and growth [49,50]. Thus, we performed genomic collinearity analysis, and found that the moso bamboo genome has 32 pairs of duplicated PeGLK genes, including 31 segmental and one tandem duplicated pair. Segmental duplication was obviously the major method of expansion of the moso bamboo GLK gene family. Segmental duplication has been shown to be more universal than tandem replication and to play a significant role in long-term evolution [22,[51][52][53]. Otherwise, syntheny analysis of the moso bamboo genome with maize sequenced plant genomes indicated significant collinearity of family members between bamboo and monocot maize, which is in accord with the evolutionary relationship between dicotyledons and monocotyledons.
To get more insight into the Phyllostachys edulis macroevolutionary profile and divergence times, Ks and Ka for paralogous and orthologous gene pairs were estimated and Ks and Ka/Ks values were calculated for each gene pair. The Ks values showed that Phyllostachys edulis had a large-scale duplication event about 15 MYA, and the divergence time for orthologous gene pairs (Pe-Zm) was approximate 38 MYA. It has also been shown that a whole genome duplication event occurred in Phyllostachys edulis 7-12 MYA and the divergence time between moso bamboo and maize was 35-42 MYA [33]. Compared with the previous reports, our results suggest that the moso bamboo GLK gene family underwent an earlier large-scale duplication event and diversified earlier than maize. In general, the Ka/Ks ratio was available to determine the selective pressure of choice acts on the coding sequences [37]. In the current study, the Ka/Ks ratios for Pe-Pe and Pe-Zm gene pairs were both <1, which imply that the homologous PeGLK gene pairs experienced strong purifying selection during the evolutionary process [26,43].

Potential Functions of PeGLKs in Moso Bamboo Development and Stress Responses
GLK genes in maize play different roles in the growth and development of chloroplasts in mesophyll cells and bundle sheaths, and they also guide chloroplast development into monomorphic form in Arabidopsis [14,17]. Here, we examined the expression levels of 36 putative PeGLK genes in four different tissues of Phyllostachys edulis and found that great majority genes exhibited high expression levels in the leaf. In view of the important function of chloroplasts in the conversion of light energy into chemical energy in plant growth and development [54], these results imply that PeGLK might be involved in the development of chloroplasts, and the high expression level in leaves might satisfy the active metabolism or gene activity. Moreover, ZmGLK47 was reported to have high expression in all maize tissues, related to the formation and development of chloroplasts [22]. PeGLK67, as an ortholog to ZmGLK47 in Phyllostachys edulis, exhibits the same protein structure and expression patterns.
Plant genomes include various stress-related genes which enable them to respond to diverse environmental conditions [54]. Additionally, the cis-elements of promoter regions largely determine the stress-response gene expression profiles that help plants adapt to adverse conditions, and are also closely linked with multiple stimulus-response genes [55,56]. In view of the function of cis-elements in abiotic stress resistance, we examined the cis-elements of PeGLKs and discovered a high number of DRE, LTRE, MeJA, and GA responsive sequences in PeGLK promoters, which indicates the significant role of PeGLKs in stress responses [22]. Plants produce many signaling molecules in response to abiotic stresses and phytohormone, treatments such as low temperature, drought, salinity, ABA, MeJA, GA, and SA. To date, an association of plant GLK with abiotic and phytohormone stress responses has been discovered in many species [6,[22][23][24]. Thus, we examined the expression of 13 selected PeGLK genes in response to cold, osmotic, MeJA, and GA stress in Phyllostachys edulis. Previous studies have shown that orthologous genes in different species were conservative in their gene functions, and paralogous genes present different functions owing to gene duplication [57]. Compared with previous studies, we found that the expression of PeGLK36 and its ortholog in maize, ZmGLK1, exhibited similar patterns in their cold and osmotic stress responses [22]. Conversely, the expression of PeGLK67 and its ortholog in maize, ZmGLK50, showed opposite patterns, indicating that PeGLKs may have lost or acquired functions in the process of evolution. Additionally, previous reports have also shown that AtGLK1 may be involved in the SA and JA signaling pathway in disease prevention [19], and OsGLK1 is related to resistance to pathogen invasion [58]. Further research on PeGLKs might provide more insight into stress-resistant and disease-resistant Phyllostachys edulis varieties.

Conclusions
In this study, 78 moso bamboo GLK family genes were identified and could be divided into 11 subfamilies. Members of the same subfamily displayed similar gene structures and motif compositions, implying that PeGLK genes were conserved in the process of evolution. Furthermore, 78 PeGLK genes were unevenly distributed across the 23 moso bamboo chromosomes. Evolutionary analysis showed that segmental duplication was the main reason for the expansion of moso bamboo GLK gene family. The expression profiles of moso bamboo GLK genes indicated that PeGLKs play significant roles in different tissues of moso bamboo growth and development, and PeGLKs was highly expressed in leaf than in other tissues. The expression patterns of PeGLK genes under cold stress, osmotic stress, and MeJA and GA treatments offer a different viewpoint about the function of moso bamboo GLKs in responses to different abiotic stress and phytohormones treatments. Taken together, these results may be helpful to select suitable candidate genes for further cloning and provide valuable information for functional analysis of PeGLK genes.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27175491/s1, Figure S1: Multiple sequence alignment of the moso bamboo GLK conserved domain. Identical residues are shown in black shadow and similar residues in pink; Figure S2: Percentage of PeGLKs response to abiotic stress and hormone treatments. Up and Down represent up and down regulated after cold stress, osmotic stress, MeJA, and GA treatments, respectively. No represents no expression change after treatments; Table S1: Detailed information about ZmGLK genes in maize; Table S2: The MEME motif sequences and lengths of PeGLK proteins in moso bamboo; Table S3: The microarray data of 36 PeGLK genes in Phyllostachys edulis; Table S4: The primers of qRT-PCR of 37 PeGLKs and PeTIP41; Table S5: Expression level of PeGLK genes in response to cold stress, osmotic stress, MeJA and GA treatments.