Comprehensive Analysis of Glutamate Receptor-like Genes in Rice (Oryza sativa L.): Genome-Wide Identification, Characteristics, Evolution, Chromatin Accessibility, gcHap Diversity, Population Variation and Expression Analysis

Glutamate receptors (GLR) are widely present in animals and plants, playing essential roles in regulating plant growth, development and stress response. At present, most studies of GLRs in plants are focused on Arabidopsis thaliana, while there have been few studies on rice. In this study, we identified 26 OsGLR genes in rice (Oryza sativa L.). Then, we analyzed the chromosomal location, physical and chemical properties, subcellular location, transmembrane (TM) helices, signal peptides, three-dimensional (3D) structure, cis-acting elements, evolution, chromatin accessibility, population variation, gene-coding sequence haplotype (gcHap) and gene expression under multiple abiotic stress and hormone treatments. The results showed that out of the 26 OsGLR genes, ten genes had the TM domain, signal peptides and similar 3D structures. Most OsGLRs exhibited high tissue specificity in expression under drought stress. In addition, several OsGLR genes were specifically responsive to certain hormones. The favorable gcHap of many OsGLR genes in modern varieties showed obvious differentiation between Xian/indica and Geng/japonica subspecies. This study, for the first time, comprehensively analyzes the OsGLR genes in rice, and provides an important reference for further research on their molecular function.


Introduction
In 280 BC, Aristotle asserted that plants also have minds and can perceive the world as animals do. Recent studies have demonstrated that plants can indeed use their "wisdom" to adapt to changes in the external environment. With the rapid development of electrophysiology as well as cell and molecular biology, numerous molecules in plants have been found to perform neurological-like functions, which paves the way to explore the mystery of "plant neurology". Glutamate acts as an excitatory neurotransmitter in the vertebrate central nervous system, facilitating long-range information exchange via activation of glutamate receptor channels [1]. Glutamate receptor-like (GLR) genes are widely found in animals and plants. There are various GLR genes in animal neurons and glial cells, which can be mainly divided into ionotropic glutamate receptors and metabotropic glutamate receptors.
Ionotropic glutamate receptors (iGluRs), such as kainate, AMPA and NMDA receptors, are glutamate-gated ion channels mainly expressed in the brain [2]. The iGluRs are complete membrane proteins consisting of four large subunits (>900 residues), including the extracellular amino-terminal domain (ATD) with a signal peptide, the extracellular ligand-binding domain (LBD), the transmembrane domain (TMD) with M1-M4 helixes, and the intracellular carboxyl-terminal domain (CTD), which form a central ion channel hole [2]. Plant glutamate receptor-like (GLRs) are specific membrane proteins with ligand-gated ion channel activities. It has been demonstrated that plant GLRs have a high degree of homology with non-NMDA receptors [3]. Fully functional plant GLRs have similar structures to iGluRs in animals, which also consist of four subunits [4]. Since it is difficult to study these GLRs in humans due to the complexity and ethical concerns, studies of them in plants may facilitate a better understanding of the more complex functions such as memory, learning or neurodegenerative diseases from the perspective of fundamental cellular biology.
Rice (Oryza sativa L.) is an important model species for monocot plants and cereals, such as maize (Zea mays), wheat (Triticum aestivum L.), barley (Hordeum vulgare L.) and sorghum (Sorghum bicolor). Therefore, it is of great significance to carry out basic research on GLR genes in rice. However, there have been few studies of GLR-related studies in rice. Li et al. [5], for the first time, screened an osglr3.1 mutant from a rice T-DNA insertion library, and osglr3.1 was found to be critical for the division and survival of individual cells in root apical meristems. In 2016, Ni et al. [6] identified 13 OsGLR genes solely through the Phytozome website (http://www.phytozome.net/, accessed on 29 July 2022). In this study, we identified 26 OsGLR genes in Oryza sativa by BLAST, HMM search and literature analysis [5][6][7], which genes were then divided into four subgroups. Subsequently, a comprehensive analysis of OsGLR genes was carried out, including the chromosomal location, physical and chemical properties, subcellular location, phylogenetic relationship, conserved motif, transmembrane (TM) helices, signal peptides, three-dimensional (3D) structure, cis-acting elements, evolution, chromatin accessibility, population variation, gene-coding sequence haplotype (gcHap) and gene expression analysis under multiple abiotic stress and hormone treatments. The data provide an important reference for further exploring the molecular functions of OsGLR genes in regulating rice growth and stress response.

Identification and Chromosomal Location of OsGLR Genes
To identify all GLR genes in the whole genome of rice, we first downloaded the whole genome data of rice (Oryza sativa Geng/japonica), tomato (Solanum lycopersicum) and Arabidopsis thaliana from the Ensembl plant database (http://plants.ensembl.org/ index.html, accessed on 29 July 2021), and constructed a local database with protein sequence files. Then, we downloaded the GLR protein (ID: PF00060), searched PF00060 with hidden Markov models (HMM) from the Pfam database (http://pfam.xfam.org/, accessed on 29 July 2021) and used the HMM-search program of the HMMER software package (European Bioinformatics Institute, Cambridge, UK) to screen protein sequences containing this domain in rice. By using Pfam, SMART (http://smart.emblheidelberg.de/, accessed on 29 July 2021) and the NCBI-CDD database (https://www.ncbi.nlm.nih.gov/ structure/cdd/cdd.shtml, accessed on 29 July 2021) to detect candidate proteins, and after removal of protein sequences containing incomplete domains, and analysis of the GLRrelated literature, we finally identified all OsGLR genes [6], renaming them according to their positions on the chromosome.

Structural Characteristics and Phylogenetic Analysis of OsGLR Genes
We first used the software TBtoolsv0.66833 (South China Agricultural University, Guangzhou, China) to extract the protein sequences of the OsGLR genes from the rice genome data [8], and we predicted the physical and chemical properties of all OsGLRs on the website Expasy (https://web.expasy.org/protparam/, accessed on 22 July 2021). After extracting the position of the UTRs, introns, and exons of OsGLRs from the RAP-BD (the Rice Annotation Project Database) website, the structural information for OsGLRs was obtained and visualized in GSDS [9] (http://gsds.gaolab.org/, accessed on 22 July 2021). WoLF PSORT can predict the subcellular localization of proteins based on their amino acid sequences. Therefore, the subcellular localization of all OsGLR genes was completed on the online website (https://wolfpsort.hgc.jp/, accessed on 22 July 2021). For phylogenetic tree analysis, we downloaded the genome sequences of GLRs in rice (Oryza sativa Japonica Group), tomato (Solanum lycopersicum), and A. thaliana from the Ensembl plant database (http://plants.ensembl.org/index.html, accessed on 22 July 2021). MEGA [10] (Molecular Evolutionary Genetics Analysis,) is a very powerful software for molecular evolutionary genetic analysis, which can be used for sequence alignment, evolutionary tree inference, estimation of the molecular evolution rate and validation of evolutionary hypotheses. The phylogenetic tree of OsGLR, SlGLR and AtGLR sequences was constructed using the ML (maximum likelihood) method in MEGA7, and beautified on the online website (iTOL https://itol.embl.de/, accessed on 22 July 2021).

Prediction of Transmembrane Helices and Signal Peptide of the OsGLR Genes
Prediction of transmembrane (TM) helices was performed with reference to the guide for the online website TMHMM Server v.2.0 (https://services.healthtech.dtu.dk/ service.php?TMHMM2.0, accessed on 23 July 2021). Protein sequences of the 26 Os-GLRs were converted into FASTA format and submitted using the default parameters on the online website. The SignalP 5.0 server can predict the presence of signal peptides and the locations of their cleavage sites in proteins [11]; it is based on a deep convolutional and recurrent neural network architecture including a conditional random field. The deep recurrent neural network architecture is more powerful in recognizing sequence motifs with varying lengths, such as signal peptides, than traditional feedforward neural networks. The prediction method was similar to TM prediction. Protein sequences of the 26 OsGLRs were converted into FASTA format and input into the online website (https://services.healthtech.dtu.dk/service.php?SignalP5.0, accessed on 23 July 2021). Eukarya was selected as the Organism Group and the other parameters were the default values.

Three-Dimensional (3D) Structure and Conserved Motifs of OsGLR Genes
The 3D structure can provide a wealth of information on the biological function and evolutionary history of macromolecules. It can be used to examine the relationship between the sequence structure and function, protein interactions and active sites. The protein sequence of each OsGLR was used to predict their 3D structure with the Phyre2 (Structural Bioinformatics Group, Imperial College, London, UK) [12] (Protein Homology/analogY Recognition Engine V 2.0) online website (http://www.sbg.bio.ic.ac.uk/phyre2/html/ page.cgi?id=index, accessed on 16 July 2021). The pipeline mainly includes the following steps: 1. detection of sequence homologs with PSIBlast; 2. prediction of secondary structure and disorder with Psipred and Disopred; 3. construction of an HMM of the sequence based on the homologs detected above; 4. scanning of this HMM against a weekly updated library of HMMs of proteins with experimentally solved structures; 5. construction of 3D models of the protein based on the alignment between the HMM of the sequence and the HMMs with known structures; 6. modelling of insertions and deletions using a loop library, a fitting procedure (cyclic coordinate descent) and a set of empirical energy terms; 7. modelling of amino acid side chains using a rotamer library from Roland Dunbrack's laboratory and our own implementation of a fast graph-based approach (R3) to optimize the choice of rotamer for each sidechain whilst avoiding steric clashes; 8. submission of the top model (if sufficiently confident) for binding site prediction by 3DLigandSite(Structural Bioinformatics Group, Imperial College, London, UK); 9. prediction of transmembrane helix and topology by memsatsvm. Finally, the predicted 3D structure results were visualized in PyMOL 4.5.0 software (Warren Lyford DeLano, Philadelphia, PA, USA) [13]. MEME (Multiple Em for Motif Elicitation) (National Institutes of Health, Bethesda, MD, USA) can discover novel, ungapped motifs (recurring, fixed length patterns) and split variable-length patterns into two or more separate motifs [14]. For the motif distribution in OsGLR sequences, we selected Zero or One Occurrence per Sequence (zoops) and 10 conserved motifs. All other parameters were set as their default values (https://memesuite.org/meme/tools/meme, accessed on 16 July 2021).

Cis-Acting Regulatory Elements, Functional Interaction Network and Chromatin Accessibility of OsGLR Genes
Prediction of cis-acting regulatory sites was carried out through the online website (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/, Last accessed on 26 July 2021). The whole genome sequence of rice was used to extract the 2000 bp sequence before the target gene promoter in TB tools v0.66833 Toolkit [8] (South China Agricultural University, Guangzhou, China). Four genes were not in the annotation file, which were downloaded separately from the Ensembl Plant Database. On the PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/, accessed on 26 July 2021) website, the cis-acting regulatory elements of OsGLRs were predicted, and the results were visualized in Gene Structure Display Server (GSDS) [9] (http://gsds.gaolab.org/, accessed on 26 July 2021). The functional interaction network analysis was completed according to the instructions on the online website (https://www.stringdb.org/, accessed on 26 July 2021). STRING used a spring model to generate the network images. Nodes were modeled as masses and edges as springs. The final positions of the nodes in the image were computed by minimizing the 'energy' of the system. The amino acid sequences of the 26 OsGLRs were used to predict the functional interaction protein on the website STRING (https://www.stringdb.org/, accessed on 26 July 2021). Prediction of variant effects of chromatin accessibility was carried out through the online website (http://ricevarmap.ncpgr.cn/effect_predict/, accessed on 26 July 2021). RiceVarMap v2.0 [15] (Chinese Academy of Agricultural Sciences, Peking, China) is a comprehensive database for rice genomic variation and its functional annotation, providing curated information for 17,397,026 genomic variations (including 14,541,446 SNPs and 2,855,580 small INDELs) from the sequencing data of 4726 rice accessions. For chromatin accessibility, we collected six tissues (root (RT), young leaf (YL), flag leaf (FL), young panicle (YP), lemma and palea (LP), and stamen and pistil (SP)) of Zhenshan 97 (a Xian/indica variety) for the ATACseq experiment, with at least two replicates for each tissue. After mapping to the reference genome, an average of 39.9 million qualified fragments were obtained per sample. The chromatin accessibility results of OsGLRs were downloaded and plotted.

Population Variation, Gene-Coding Sequence Haplotype (gcHap) and KAKS of OsGLR Genes
Population variation analysis was carried out based on resequencing of 3010 rice accessions [16], which were divided into five subspecies, including GJ (Geng/japonica, 801 accessions), XI (Xian/indica, 1764 accessions), AUS (aus/boro, 221 accessions), ARO (aromatic basmati/sadri, 101 accessions) and ADM (admixed, 123 accessions). All these five subspecies were further grouped into 12 groups according to the classification of their corresponding rice accessions. These groups included five subgroups (XI1-5) of the XI subspecies AUSG6, four subgroups (GJ7-10) of the GJ subspecies AROG11, and admixtures (ADM). All genes were categorized according to their presence/absence in the 453 highquality accessions. Core genes were defined as genes existing in all high-quality rice accessions. Distributed genes were defined as genes existing in significantly less than 99% of accessions (binomial tests, p-value < 0.05, null hypothesis is "loss rate < 1%"). Candidate core genes were taken as those existing in > 99% (not all) of high-quality rice accessions (binomial test, fdr < 0.05). Random genes were defined as genes without differences among groups and subgroups. By using the identified OsGLRs for gene distribution and presence frequency analysis, population distribution results were obtained from the full dataset of 32 million Nipponbarebased SNPs from 3010 rice accessions (http://cgm.sjtu.edu.cn/3kricedb/search.php, accessed on 26 April 2021). Haplotype information [17] was obtained from the 3K database website (https://www.rmbreeding.cn/Index/, accessed on 26 April 2021), and modern varieties were defined in the literature [18]. Each gcHap of the OsGLRs genes was constructed by the concatenating SNPs in the CDS region, in which synonymous SNPs were ignored. KAKS was calculated in TBtools Toolkit (South China Agricultural University, Guangzhou, China) using the protein sequence and CDS sequence of OsGLRs, and the divergence time (Mya) was calculated using the formula: Mya = Ks/2λ × 10 −6 , where λ = 6.5 × 10 −9 .

Expression Difference of OsGLR Genes under Different Hormone Treatments
Gene expression data generated by the Affymetrix ATH1 (Thermo Fisher Scientific, Massachusetts, USA) array were normalized by the GCOS method, with a trimmed mean target intensity (TGT) value of 100. Most tissues were sampled in triplicate. The tissue expression data of all the OsGLR genes were downloaded from the BAR (The Bio-Analytic Resource for Plant Biology) database (http://www.bar.utoronto.ca/, accessed on 8 October 2021), and processed and visualized in TBtools Toolkit using the default settings. Gene expression data under abiotic stress conditions were obtained from http://rice.uga.edu/expression.shtml, accessed on 8 October 2021. The values of presence/absence variation were assigned for digital gene expression (DGE) libraries, and genes were regarded as 'expressed' if at least one sequence read was mapped uniquely within an exon.
The seeds of rice (Nipponbare) were sterilized with 70% ETOH and 1% hypochlorous acid. After disinfection, the seeds were imbibed and germinated in distilled water. The germinated seeds were transplanted into a 96-well plate and grown in a growth chamber at 28 • C, 24 h of light and 60% humidity. After three days of germination, the seedlings were transferred to the medium diluted five times with Yoshida nutrient solution [19] and adjusted to pH 5.5 with 1 mMMES + 1 M NaOH. The seedlings were allowed to grow for another four days. The medium was updated every two days. The 7-day-old seedlings were transferred to a medium containing plant hormones (50 µM ABA + 0.1% Ethanol, 10 µM GA + 0.1% Ethanol, 10 µM IAA + 0.1% Ethanol, 1 µM BR + 0.1% Ethanol, 1 µM CK + 0.02% DMSO and 100 µM JA + 0.02% DMSO). Shoot samples were collected at 0, 1, 3, 6 and 12 h after hormone treatment. Root samples were collected at 0 min, 15 min, 30 min, 1 h, 3 h and 6 h after hormone treatment. TRIzol (TIANGEN, W9330, Peking, China) was used to extract RNA from a total of 138 samples, and a Takara reverse transcription kit (SparkJade, AG0304, Shandong, China) was used to reverse message RNA into cDNA, which was diluted by five times and amplified by PCR MIX (SparkJade, SMQFK, Shandong, China). These RNA were labeled with Cy3 (mock treatment) and Cy5 (hormone treatment), and used for hybridization using the Agilent two-color microarray analysis system. The time-course expression profile for each gene is shown as the log-ratio of the signal intensity (log2 Cy5/Cy3).

Semi-Quantitative PCR Analysis of OsGLR Genes under Low-Temperature Stress
The dormancy of Nipponbare (Oryza sativa L. japonica) seeds was broken in a drying oven (Boxun Ltd., GZX9240MBE, Hefei, China) at 43 • C for 5 days, and then 250 round and plump seeds of Nipponbare were selected and disinfected with 5% hypochlorous acid for 1 h. After disinfection, the sample was rinsed with ddH 2 O four to five times, and then put in 10 × 9 cm glass petri dishes with two layers of germination paper, followed by the addition of 10 mL distilled water for imbibing and germination. The germinated rice seeds were transplanted into a 12 × 8 cm germination box containing Hoagland nutrient solution. The nutrient solution was changed every four days, and cultivation was carried out under natural conditions. When the seedlings are cultivated to the three-leaf stage, low-temperature (4 • C) stress treatment was performed. Leaf samples were collected in liquid nitrogen at 0 h and 3 h, 6 h, 9 h, 12 h, 24 h and 48 h after low-temperature stress treatment, with sampling of seven seedlings at a time, and the experiment was repeated three times. The collected rice leaf samples were stored in a refrigerator (Haier, DW86L486, Qingdao, China) at -80 • C for later use. The primers of OsGLR genes are shown in Table 1. TRIzol (TIANGEN, W9330, Beijing, China) was used to extract RNA from rice leaf tissue, and a Takara reverse transcription kit (SparkJade, AG0304, Shandong, China) was used to reverse message RNA into cDNA, which was diluted five times and amplified by PCR MIX (SparkJade, SMQFK, Shandong, China). Table 1. Semi-quantitative primers of the glutamate receptor-like (GLR) genes in rice.

Genome-Wide Identification and Characterization of the OsGLR Genes in Rice
The OsGLR genes in the whole genome of Oryza sativa were identified with BLAST, HMM search and literature analysis [5,6]. Finally, a total of 26 OsGLR genes were identified in the genome, which were renamed according to their location information on chromosomes ( Table 2). To characterize these putative OsGLR genes, we further identified their physical and chemical characteristics, including their start and stop locations, protein lengths, theoretical PI and subcellular localization. The coding sequence length of the 26 OsGLR genes ranged from 361 bp to 9868 bp, and the protein length ranged from 95 to 988 amino acids ( Table 1). The molecular weight of the OsGLR genes was predicted to be between 10497.73 and 107252.1 Da ( Table 1). The PI (theoretical PI) ranged from 4.05 to 9.52, with an average of 6.52 (Table 1). Subcellular localization revealed that the 26 OsGLR genes were located in cytoplasm, chloroplast, endoplasmic reticulum, mitochondrion, plasma membrane, vacuole and nucleus. Almost all genes were predicted to be located on the plasma membrane or chloroplast (Table 1). Different physical and chemical characteristics of the OsGLR genes or proteins indicated their different biological functions.
The 26 OsGLR genes were unevenly located on chromosomes 2, 4, 6, 7 and 9 ( Figure 1a). There was only one gene on chromosomes 2 and 4, two genes on chromosomes 7 and 12, and nine genes on chromosomes 6 and 9 (Figure 1a), most of which were clustered at the end of chromosomes. To further reveal the evolution and structural characteristics of the OsGLR genes, we constructed a phylogenetic tree of 56 GLR genes using the ML (maximum likelihood) method, including 26 GLR genes in rice (Oryza sativa Japonica Group), 20 GLR genes in A. thaliana and 10 GLR genes in tomato (Solanum lycopersicum) (Figure 1b).  12, and nine genes on chromosomes 6 and 9 (Figure 1a), most of which were clustered at the end of chromosomes. To further reveal the evolution and structural characteristics of the OsGLR genes, we constructed a phylogenetic tree of 56 GLR genes using the ML (maximum likelihood) method, including 26 GLR genes in rice (Oryza sativa Japonica Group), 20 GLR genes in A. thaliana and 10 GLR genes in tomato (Solanum lycopersicum) ( Figure  1b). The conserved motif can provide insights into how patterns of residue conservation and divergence in a protein family relate to functional properties, and can provide useful links to more detailed information that may be helpful in understanding that those sequences have all four subfamily/structure/function relationships. There are four subfamilies, two large subfamilies and two small subfamilies. Arabidopsis thaliana has only three subfamilies (groups 1, 2 and 4).
The 56 GLR genes were divided into two large subfamilies and two small subfamilies. The OsGLR, AtGLR and SlGLR genes were present in four, three and two subfamilies, respectively. The third subfamily was unique to the OsGLR genes, including eight OsGLR genes, namely, OsGLR6.3, OsGLR6.4, OsGLR6.5, OsGLR6.6, OsGLR6.7, OsGLR6.8, OsGLR6.9 and OsGLR6.10 ( Figure 1b and Table 2), which were mostly located on chromosome 6. There were four GLR genes for rice (OsGLR2. 2 Figure 1. Characteristics of the glutamate receptor-like (GLR) genes in rice. (a) Chromosome location of the OsGLR genes; (b) phylogenetic tree of GLR genes in rice, tomato and Arabidopsis thaliana; (c) intron-exon distribution of the OsGLR genes; (d) conserved motif prediction of the OsGLR genes. The conserved motif can provide insights into how patterns of residue conservation and divergence in a protein family relate to functional properties, and can provide useful links to more detailed information that may be helpful in understanding that those sequences have all four subfamily/structure/function relationships. There are four subfamilies, two large subfamilies and two small subfamilies. Arabidopsis thaliana has only three subfamilies (groups 1, 2 and 4).
To understand the structural characteristics of the OsGLR genes, the intron-exon structure of the 26 OsGLR genes was analyzed using the GSDS toolkit [9]. The OsGLR genes had 2-7 exons and 1-6 introns (Figure 1c). Different OsGLR genes have different structures, but it is worth noting that the genes from the same subfamily had the same or similar intron-exon structure. For example, OsGLR9.1, OsGLR9.3, OsGLR9.4, OsGLR9.5 and OsGLR9.6 in subfamily II had five exons and four introns. OsGLR6.10, OsGLR6.4, OsGLR6.8 and OsGLR6.9 in subfamily III also had five exons and four introns. OsGLR4.1, OsGLR6.11 and OsGLR9.9 located in subfamily IV harbored two exons and one intron (Figure 1c). We also predicted the important motifs of the OsGLR genes (Figure 1d and Figure S1). MEME (Multiple Em for Motif Elicitation) analysis showed that the OsGLR genes had 10 relatively conserved motifs ranging from 21 to 50 amino acids ( Figure S1). Most of the OsGLR genes had ten, nine or six motifs, which was similar to the intronexon structure. The OsGLR genes in the same subfamily had the same or similar motifs. OsGLR6.3, OsGLR6.4, OsGLR6.8 and OsGLR6.9 in subfamily III had six identical or similar motifs. OsGLR9.3, OsGLR9.4, OsGLR9.5, OsGLR9.6 and OsGLR9.8 had ten identical or similar motifs. OsGLR6.12 and OsGLR7.1 had nine identical or similar motifs. According to the phylogenetic tree, structure and motif analysis of OsGLRs, OsGLR9.1, OsGLR9.3, OsGLR9.4, OsGLR9.5 and OsGLR9.6 in subfamily II and OsGLR6.4, OsGLR6.8 and OsGLR6.9 in subfamily III had the same or similar intron-exon structure and motifs. These genes might have the same or similar biological functions, while different intron-exon structures and motifs within the same subfamily might also be involved in different functions.

Transmembrane (TM) Helices and Signal Peptides of OsGLR Genes
To understand the structural characteristics of the OsGLR genes, we further predicted their transmembrane (TM) helices and signal peptides (Figure 2 and Figure S2). The OsGLR genes could be divided into two types according to TM helices: TM-related Os-GLR genes and non-TM OsGLR genes. There were six genes without the TM domain (OsGLR4. 1

Three-Dimensional (3D) Structure of OsGLR Genes
The 3D structure can provide important information such as sequence-structurefunction relationships, interactions and active sites. By using the Phyre2 web portal [12], we could identify the distinct biological units, interactions among molecular components and secondary structures (alpha helices and beta strands) and 3D domains within individual protein molecules. The results showed that there were significant differences in the 3D structure of different OsGLRs

Three-Dimensional (3D) Structure of OsGLR Genes
The 3D structure can provide important information such as sequence-structure function relationships, interactions and active sites. By using the Phyre2 web portal [12 we could identify the distinct biological units, interactions among molecular component and secondary structures (alpha helices and beta strands) and 3D domains within ind vidual protein molecules. The results showed that there were significant differences in th 3D structure of different OsGLRs (Figure 3). Interestingly, 10 OsGLR genes with both a TM domain and signal peptides had similar 3D structures, including OsGLR2. 1 Figure 3).

Cis-Acting Regulatory Elements and Functional Interaction Network of the OsGLR Genes
To predict the function of the OsGLR genes, first we extracted their 2000-bp upstream promoter regions from the rice GJ genome to identify the cis-acting regulatory elements. A total of 23 cis-regulatory elements were found (Figure 4), including ABRE, ACE, AT-rich, AUXRR-core, Box, C-box, CAAT-box, G-box, TC-rich, Gare-motif, GC-motif, GCN4_motif, GT1-motif, GTGGC-motif, I-B, Ox, LTR, MBS, MRE, MSA-like, Sp1, TatA-box and TATC-Box, which could be divided into four categories, namely, light responsive elements, growth and development related elements, plant-hormone-related elements and stressrelated elements. ABRE is a transcription factor, as well as the main cis-element of ABA pathway participating in drought, salt and osmotic-stress response [22,23]. The CAATbox is mainly involved in photosynthesis [24]. The G-box plays an important role in the early senescence of rice leaves [25]. The GCN4_motif reduces protein synthesis [26]. The ten OsGLR genes with both the TM domain and signal peptides and similar 3D protein structures were further analyzed. The ten genes had similar structures but different cisacting elements. The tata-box was a cis-acting regulatory element shared by all 10 of these OsGLR genes. Except for OsGLR9.1, the other nine genes all had the CAAT-box, and all genes had the G-box except for OsGLR9. 6 Figure 3. Three-dimensional (3D) structure prediction of the glutamate receptor-like (GLR) genes in rice. Green helices represent α-helices, blue arrows indicate β-strands and faint lines indicate coil Confidence is colored from high (red) to low (blue) using a rainbow spectrum.

Cis-Acting Regulatory Elements and Functional Interaction Network of the OsGLR Genes
To predict the function of the OsGLR genes, first we extracted their 2000-bp upstream promoter regions from the rice GJ genome to identify the cis-acting regulatory elements. A total of 23 cis-regulatory elements were found (Figure 4), including ABRE, ACE, ATrich, AUXRR-core, Box, C-box, CAAT-box, G-box, TC-rich, Gare-motif, GC-motif, GCN4_motif, GT1-motif, GTGGC-motif, I-B, Ox, LTR, MBS, MRE, MSA-like, Sp1, TatAbox and TATC-Box, which could be divided into four categories, namely, light responsive elements, growth and development related elements, plant-hormone-related elements and stress-related elements. ABRE is a transcription factor, as well as the main cis-element of ABA pathway participating in drought, salt and osmotic-stress response [22,23]. The CAAT-box is mainly involved in photosynthesis [24]. The G-box plays an important role in the early senescence of rice leaves [25]. The GCN4_motif reduces protein synthesis [26]. The ten OsGLR genes with both the TM domain and signal peptides and similar 3D protein structures were further analyzed. The ten genes had similar structures but different cis-acting elements. The tata-box was a cis-acting regulatory element shared by all 10 of these OsGLR genes. Except for OsGLR9.1, the other nine genes all had the CAAT-box, and all genes had the G-box except for OsGLR9.6. OsGLR2.1, OsGLR2.2, OsGLR7.1, OsGLR9.1  Confidence is colored from high (red) to low (blue) using a rainbow spectrum.

Chromatin Accessibility of OsGLR Genes
Since the OsGLR genes with the same structure had differences in function, we analyzed the chromatin accessibility of OsGLR genes to determine whether epigenetics affects their functions. We predicted the variant effects of these genes on chromatin through the RiceVarMap v2.0 toolkit [15]. The results showed that 19 OsGLR genes had significant differences in chromatin accessibility ( Figure 6). OsGLR6.3, OsGLR6.4, OsGLR6.5, OsGLR6.9, OsGLR6.10 and OsGLR9.7 had the highest chromatin accessibility mainly in roots; Os-GLR2.2, OsGLR4.1 and OsGLR9.4 exhibited the highest chromatin accessibility mainly in flag leaves; OsGLR6.1, OsGLR6.12 and OsGLR7.1 had the highest chromatin accessibility mainly in young leaves; and OsGLR2.1, OsGLR7.2 and OsGLR9.8 had the highest chromatin accessibility mainly in lemma (palea). The highest chromatin accessibility of OsGLR6.2, OsGLR6.11, OsGLR9.6 and OsGLR9.9 was mainly found in young ears. In particular, the chromatin accessibility of OsGLR2.2 had a large number of peaks, which was significantly more than other OsGLR genes ( Figure 6). Moreover, we further analyzed the ten OsGLR genes with both the TM domain and signal peptides and similar 3D protein structures, and they had the same interacting proteins except for OsGLR2.2. Their common interacting proteins were AMT3-2, AMT3-1, AMT2-3, AMT2-2, AMT2-1, AMT1-1, OsJ_07562, OsJ_07561, TPC1 and OEP21, which probably work through involvement in ammonium transport ( Figure 5). Prediction of the cis-acting regulatory elements and interacting proteins of the OsGLR genes revealed that the same structure might represent the same function. However, there were sometimes exceptions. For example, OsGLR2.2 had a similar structure but different interacting proteins compared with the other nine proteins.

Chromatin Accessibility of OsGLR Genes
Since the OsGLR genes with the same structure had differences in function, we analyzed the chromatin accessibility of OsGLR genes to determine whether epigenetics affects their functions. We predicted the variant effects of these genes on chromatin through the RiceVarMap v2.0 toolkit [15]. The results showed that 19 OsGLR genes had significant differences in chromatin accessibility ( Figure 6). OsGLR6. 3 OsGLR6.11, OsGLR9.6 and OsGLR9.9 was mainly found in young ears. In particular, the chromatin accessibility of OsGLR2.2 had a large number of peaks, which was significantly more than other OsGLR genes ( Figure 6).

Evolutionary, Population Variation and Gene-Coding Sequence Haplotype (gcHap) Analysis of the GLR Genes in Rice
Evolutionary and population variation analysis were carried out with RPAN [28] (Rice Pan-genome Browser, (Chinese Academy of Agricultural Sciences, Beijing, China)) and the RFGB [17] (Rice Functional Genomics and Breeding) database. According to the pangenomic analysis of OsGLR genes, we obtained the distribution and presence frequency of the OsGLR genes in 453 high-quality accessions (Figure 7a,b and Figure S3). Moreover, according to the heat map of the OsGLR gene frequency in different subspecies and subgroups of 453 high-quality accessions, we counted the core, candidate core and distributed OsGLR genes. OsGLR6.4, OsGLR6.5, OsGLR6.11 and OsGLR9.7 were present in all 453 high-quality accessions ( Figure 7a). As for the distribution of genes, the OsGLR genes accounted for a higher proportion in japonica subspecies (Figure 7b and Figure S3) (Figure 7b and Figure S3), and were absent in other subspecies except for OsGLR9. 5

Expression Analysis of OsGLR Genes under Multiple Abiotic Stress and Hormone Treatments
Firstly, we analyzed the expression levels of 16 OsGLR genes in dry seeds, flowers, roots and young leaves (Figure 9). OsGLR2.1, OsGLR2.2, OsGLR4.1, OsGLR6.8 and Os-GLR6.9 had the highest expression level in young leaves. The expression levels of OsGLR6.3, OsGLR6.4, OsGLR6.11 and OsGLR9.6 were the highest in dry seeds. OsGLR6.10, OsGLR6.12 and OsGLR7.1 exhibited the highest expression levels in roots (Figure 9). The expression of OsGLR9.7 was the highest in flowers (Figure 9). The different expression levels of these genes in different tissues implied that they have different functions. In other plants, GLR genes are mainly involved in regulating the response to several abiotic stresses [29][30][31][32]. To identify the function of OsGLR genes, expression analysis of the OsGLR genes was carried out under multiple abiotic stress and hormone treatments ( Figure 10).    Subsequently, we analyzed the expression variation of OsGLR genes in the roots and shoots under a variety of hormone treatments (Figure 10a,b). For the shoots, the expression of OsGLR6.2 was down-regulated by treatment with ABA, BR, CK and JA ( Figure  10a). In contrast, the expression of OsGLR6.3 was induced by JA treatment and reached the highest level after 1 h. The expression of OsGLR9.8 showed slow and steady increase after ABA treatment (Figure 10a). The expression of OsGLR9.8 reached the maximum after 3 h of IAA treatment. OsGLR9.7 expression was down-regulated by IAA and GA3 treatment (Figure 10a). The expression level of OsGLR9.7 reached the highest after 6 h of ABA treatment or JA treatment, and it was almost not detected at other indicated time points (Figure 10a). The expression of OsGLR6.5 reached the highest level after 1 h of ABA, GA3, IAA and BR treatment. The expression of OsGLR6.8 was the highest after 1 h of ABA, GA3 and IAA treatment (Figure 10a). The expression level of OsGLR6.10 showed an increase immediately after IAA and BR treatment, while it decreased in the first hour and then increased to reach the highest at 6 h under ABA and GA3 treatment (Figure 10a). JA treatment up-regulated the expression of OsGLR7.1. The expression of OsGLR2.1 reached the highest after 12 h of BR treatment. OsGLR6.12 was down-regulated after CK treatment (Figure 10a). The expression of OsGLR4.1 was up-regulated after BR and CK treatment, and that of OsGLR2.2 decreased after ABA, IAA and GA3 treatment. OsGLR6.11 had low expression under normal conditions and various treatments, but was up-regulated after ABA treatment. The expression level of OsGLR9.5 was very low under various treatments and normal conditions (Figure 10a).
For the roots, OsGLR6.2 was significantly up-regulated after 1 h of ABA treatment, 6 h of GA3 treatment, after BR treatment, 3-6 h after CK treatment, and 0.5 h after JA treatment. The expression levels of OsGLR6.3 and OsGLR6.11 were very low under normal Subsequently, we analyzed the expression variation of OsGLR genes in the roots and shoots under a variety of hormone treatments (Figure 10a,b). For the shoots, the expression of OsGLR6.2 was down-regulated by treatment with ABA, BR, CK and JA (Figure 10a). In contrast, the expression of OsGLR6.3 was induced by JA treatment and reached the highest level after 1 h. The expression of OsGLR9.8 showed slow and steady increase after ABA treatment (Figure 10a). The expression of OsGLR9.8 reached the maximum after 3 h of IAA treatment. OsGLR9.7 expression was down-regulated by IAA and GA3 treatment (Figure 10a). The expression level of OsGLR9.7 reached the highest after 6 h of ABA treatment or JA treatment, and it was almost not detected at other indicated time points (Figure 10a). The expression of OsGLR6.5 reached the highest level after 1 h of ABA, GA3, IAA and BR treatment. The expression of OsGLR6.8 was the highest after 1 h of ABA, GA3 and IAA treatment (Figure 10a). The expression level of OsGLR6.10 showed an increase immediately after IAA and BR treatment, while it decreased in the first hour and then increased to reach the highest at 6 h under ABA and GA3 treatment (Figure 10a). JA treatment up-regulated the expression of OsGLR7.1. The expression of OsGLR2.1 reached the highest after 12 h of BR treatment. OsGLR6.12 was down-regulated after CK treatment (Figure 10a). The expression of OsGLR4.1 was up-regulated after BR and CK treatment, and that of OsGLR2.2 decreased after ABA, IAA and GA3 treatment. OsGLR6.11 had low expression under normal conditions and various treatments, but was up-regulated after ABA treatment. The expression level of OsGLR9.5 was very low under various treatments and normal conditions (Figure 10a).
For the roots, OsGLR6.2 was significantly up-regulated after 1 h of ABA treatment, 6 h of GA3 treatment, after BR treatment, 3-6 h after CK treatment, and 0.5 h after JA treatment. The expression levels of OsGLR6.3 and OsGLR6.11 were very low under normal conditions (Figure 10b). The expression of OsGLR6.5, OsGLR6.10 and OsGLR7.2 was limited under various hormone treatments. The expression of OsGLR2.1 was significantly down-regulated after 15 min and up-regulated after 30 min of ABA treatment (Figure 10b). OsGLR7.1 was significantly down-regulated after IAA treatment, but significantly up-regulated after ABA, JA and CK treatment (Figure 10b). The expression of OsGLR6.8 was slowly but steadily upregulated after BR treatment, and significantly up-regulated after 1 h of JA treatment. The expression of OsGLR9.8 was significantly up-regulated after ABA treatment, but showed little change under other hormone treatments. The expression level of OsGLR9.7 decreased significantly after 1 h but increased significantly after 3 h under CK treatment. OsGLR6.12 was down-regulated after ABA, IAA, BR and CK treatment (Figure 10b). The expression of  (Figure 10d). Under low-temperature stress, there were significant differences in the expression levels of OsGLR genes in roots and shoots (Figure 10d). The expression of OsGLR9.2 was the lowest at 3 h and the highest at 12 h in shoots (Figure 10d). OsGLR4.1 and OsGLR9.3 had the lowest expression levels in shoots before low-temperature treatment, and the highest expression at 12 h after treatment (Figure 10d). The expression level of OsGLR2.2 in shoots remained unchanged at 12 h and before treatment, and reached the highest at 24 h (Figure 10d). However, the expression of OsGLR9.2, OsGLR4.1, OsGLR9.3 and OsGLR2.2 in roots showed little change. The expression levels of OsGLR9.9, OsGLR6.4, OsGLR6.9, OsGLR7.1, OsGLR6.5, OsGLR9.4, OsGLR6.3 and OsGLR6.12 were the highest before low-temperature treatment (Figure 10d). The expression level of OsGLR9.9 was the highest at 1 h of low-temperature treatment. OsGLR9.4, OsGLR6.3 and OsGLR6.12 showed the highest expression before low-temperature treatment. The expression levels of OsGLR6.10, OsGLR9.8, OsGLR9.5, OsGLR9.6, OsGLR9.1 and OsGLR6.2 showed no significant change in shoots and roots before and after treatment (Figure 10d). The expression of OsGLR6.10 and OsGLR9.8 reached the highest at 6 h, while that of OsGLR9.6 and Os-GLR6.2 was the highest at 24 h, and that of OsGLR9.1 reached the highest at 3 h under low-temperature treatment (Figure 10d).

Semi-Quantitative PCR Analysis of OsGLR Genes under Low-Temperature Stress
To verify whether the OsGLR genes in Xian/indica (XI) and Geng/japonica (GJ) are different in response to low-temperature (4°C) stress, we further performed semi-quantitative PCR of four genes (OsGLR9.7, OsGLR9.3, OsGLR6.12 and OsGLR4.1) mainly distributed in Xian subspecies, and another four genes (OsGLR9.1, OsGLR7.1, OsGLR6.4 and OsGLR2.2) mainly distributed in Geng subspecies ( Figure S4). As a result, there were differences between these genes in expression under low-temperature stress. The four GJ OsGLR genes showed more obvious expression variation than the four XI OsGLR genes. For the four XI OsGLR genes, the expression level of OsGLR4.1 increased at 24 h and reached the highest level at 48 h ( Figure S4). OsGLR9.7 had a low expression level before lowtemperature stress, but its expression reached the highest level at 6 h of low-temperature stress. For OsGLR6.12 and OsGLR9.3, there was little change in expression before and after low-temperature stress. In addition, for the four GJ OsGLR genes, the expression level of OsGLR6.4 was the lowest at 6 h of low-temperature treatment, and reached the highest at 48 h. OsGLR9.1 had the highest expression before low-temperature treatment, and was almost not expressed after low-temperature treatment. The expression of OsGLR2.2 and OsGLR7.1 showed little difference before and after low-temperature treatment ( Figure S4).

Discussion
Plant GLR genes are involved in numerous physiological and biochemical processes, such as signal transduction [32,33], stress response [21,30,34], growth and development [5,29,[35][36][37], and stomata closure [38,39] in A. thaliana, tomato, radish and other species. However, this GLR gene family has not been systematically and comprehensively characterized in rice. This study for the first time systematically analyzed the OsGLR genes and explored their roles in the growth and development and stress response of rice.
In addition, the cis-acting elements of the 26 OsGLR genes could be divided into four categories, including light-responsive elements, elements related to growth and development, elements related to plant hormones, and other components (Figure 4). The light-responsive elements mainly include CAAT-box [24], ACE (angiotensin-converting enzymes), G-box, MRE (MYB-recognizing elements), and Box 4. The elements related to growth and development mainly comprise MSA-like (mitosis-specific activator) and circadian rhythm, and the elements related to plant hormone ABRE (ABA-responsive element) [23] and TATC-box [40]. There are many cis-acting elements in the promoter region of the OsGLR genes. The functions of some elements have been studied in other plant species, but the majority of them remain to be studied. The cis-acting element analysis of the OsGLR genes would help further clarify the future direction of research on this gene family.
Next, in order to explore the breeding potential of OsGLR genes, we conducted population variation and gcHap analysis based on the 3K project [16,18,28]. We found that the favorable gcHap of some OsGLR genes is only present in GJ subspecies (Figure 7 and Figure S3). The favorable gcHap of many OsGLR genes in modern varieties showed obvious differentiation between indica and japonica subspecies. For example, the favorable gcHap of the six genes, OsGLR9.9, OsGLR9.4, OsGLR9.1, OsGLR7.1, OsGLR6.4 and OsGLR2.2, was only found in GJ modern varieties (Figure 8). In 3010 rice accessions, the favorable gcHap of the genes OsGLR6.2, OsGLR6.3, OsGLR6.8, OsGLR9.1, OsGLR9.2 and OsGLR9.5 was mainly distributed in GJ subspecies (Figure 8). These genes may be more resistant to lowtemperature stress. Previous research has reported that SIGLR2.1, SIGLR4.2 and AtGLR3.4 from subfamily IV are responsive to low-temperature stress by regulating apoplastic H 2 O 2 production and redox homeostasis [21,30]. In addition, AtGLR1.2 and AtGLR1.3 in subfamily I enhance cold tolerance by increasing the endogenous JA level under cold stress [41]. Therefore, OsGLR9.9 and OsGLR7.1 in subfamily IV and OsGLR2.2 and OsGLR6.2 in subfamily I distributed in GJ subspecies are likely involved in regulating the low-temperature stress response.
Ultimately, we carried out expression analysis of the OsGLR genes under multiple abiotic stress and hormone treatments in rice (Figure 10). Under a variety of hormone treatments (ABA, GA, IAA, BR, CTK and JA), the OsGLR genes showed extremely complex expression variations (Figure 10a,b). Notably, the OsGLR6.2 gene was highly expressed in root and seedling tissues under almost all hormone treatments (Figure 10a,b). We also found that the OsGLR genes had a strong tissue specificity of expression under drought stress (Figure 10c). For example, the expression levels of OsGLR9.9, OsGLR6.10, OsGLR6.1, OsGLR6.11, OsGLR9.3, OsGLR9.1, OsGLR6.2, OsGLR9.8, OsGLR6.12 and OsGLR6.3 in roots varied with the time of drought treatment (Figure 10c), while there was no significant difference in the expression of these OsGLR genes in shoots. To test whether the OsGLR genes in indica and japonica are different in response to low-temperature stress, we carried out semi-quantitative PCR of four genes mainly distributed in indica subspecies and another four genes mainly distributed in japonica subspecies ( Figure S4). As a result, there were differences in expression between these genes under low-temperature stress, and those genes in japonica showed wider variations in expression. In short, this study carried out a systematic and comprehensive analysis of OsGLR genes, which has important reference value for follow-up, in-depth molecular function research.
Additionally, a recent paper found that plant GLRs participate in regeneration by dynamically altering chromatin and transcription in reprogrammed cells near callus, precisely regulating regeneration and defense [42]. Feijo and his collaborators studied the function of glutamate receptor-like proteins (GLR) in the basal land plant Physcomitrella patens, revealing a previously unknown and important role: controlling the movement of sperm in order to find the female reproductive organs and regulation of zygote development [43]. They found that this glutamate receptor-like protein functions as a channel that allows calcium ions to flow through. Many human glutamate receptors function in the same way, suggesting that the plant and human versions of this receptor maintained this conserved function during parallel evolution. We are passionate about using plants to study the function of animal neurons. Although the road ahead will be long and steep, we believe this to be a worthy pursuit.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cimb44120437/s1, Figure S1: The 10 predicted conserved motifs of the glutamate receptor-like (GLR) genes in rice. The length of 10 predicted conserved motifs ranged from 21 to 50 amino acids. Figure S2: Prediction of signal peptides and the location of their cleavage sites in the 26 glutamate receptor-like (GLR) proteins in rice. On the plot, three marginal probabilities are reported; the protein can have a Sec signal peptide (Sec/SPI), a Lipoprotein signal peptide (Sec/SPII), a Tat signal peptide (Tat/SPI) (depending on what type of signal peptide is predicted), CS (the cleavage site) or No signal peptide at all (Other, the probability that the sequence does not have any kind of signal peptide). Figure S3: Gene distribution of the glutamate receptor-like (GLR) genes in 3010 rice accessions. The 3010 rice accessions were divided into five subspecies, including JAP (japonica, 801 accessions), IND (indica, 1764 accessions), AUS (aus/boro, 221 accessions), ARO (aromatic basmati/sadri, 101 accessions) and ADM (admixed, 123 accessions). Figure S4: Semi-quantitative PCR analysis of the glutamate receptor-like (GLR) genes under low-temperature (4 • C) stress condition. After low-temperature stress, seeding leaves were sampled at seven different time points: 0 h, 3 h, 6 h, 9 h, 12 h, 24 h and 48 h. Table S1: Summary of favorable haplotypes of glutamate receptor genes in rice.