Functional characterization and allelic mining of OsGLR genes for potential uses in rice improvement

Glutamate-like receptor (GLR) genes are a group of regulatory genes involved in many physiological processes of plants. With 26 members in the rice genome, the functionalities of most rice GLR genes remain unknown. To facilitate their potential uses in rice improvement, an integrated strategy involving CRISPR-Cas9 mediated knockouts, deep mining and analyses of transcriptomic responses to different abiotic stresses/hormone treatments and gene CDS haplotype (gcHap) diversity in 3,010 rice genomes was taken to understand the functionalities of the 26 rice GLR genes, which led us to two conclusions. First, the expansion of rice GLR genes into a large gene family during evolution had gone through repeated gene duplication events occurred primarily in two large GLR gene clusters on rice chromosomes 9 and 6, which was accompanied with considerable functional differentiation. Secondly, except for two extremely conserved ones (OsGLR6.2 and OsGLR6.3), rich gcHap diversity exists at the remaining GLR genes which played important roles in rice population differentiation and rice improvement, evidenced by their very strong sub-specific and population differentiation, by their differentiated responses to day-length and different abiotic stresses, by the large phenotypic effects of five GLR gene knockout mutants on rice yield traits, by the significant association of major gcHaps at most GLR loci with yield traits, and by the strong genetic bottleneck effects and artificial selection on the gcHap diversity in populations Xian (indica) and Geng (japonica) during modern breeding. Our results suggest the potential values of the natural variation at most rice GLR loci for improving the productivity and tolerances to abiotic stresses. Additional efforts are needed to determine the phenotypic effects of major gcHaps at these GLR loci in order to identify ‘favorable’ alleles at specific GLR loci specific target traits in specific environments to facilitate their application to rice improvement in future.


Introduction
Glutamate receptors (GLRs) are an important gene family that are known to play important roles in a variety of physiological processes of eukaryotes.In animals, GLRs are known to act as neurotransmitters that allow animals to respond external stimuli through ion channels (Lam et al., 1998).In plants, GLRs are known to play a role in a variety of physiological processes, including growth and development (Kang and Turano, 2003;Li et al., 2006;Michard et al., 2011;Vincill et al., 2013;Kong et al., 2015;Singh et al., 2016;Wudick et al., 2018), stress responses (Lu et al., 2014;Cheng et al., 2016;Cheng et al., 2018;Li et al., 2019;Philippe et al., 2019;Wang et al., 2019;Zheng et al., 2021), and pathogen defense (Mousavi et al., 2013;Shao et al., 2020;Liu et al., 2021;Xue et al., 2022;Yu et al., 2022).However, since GLR genes represent a large gene family that encodes membrane proteins embedded within the cell membrane, studying their functionalities in plants through traditional biochemical methods has been proven difficult.
In recent years, there has been growing interest to understand the molecular structure, expression patterns and functional roles of GLR genes in plants.Specially, technical advances in structural biology have enabled researchers to gain insights into the molecular structure of plant GLRs.One breakthrough study in this regard was the high-resolution structure of a glutamate receptor-like AtGLR3.4resolved from the model plant Arabidopsis thaliana using a cryoelectron microscopy (cryo-EM) in which the tetrameric assembly of AtGLR3.4 subunits was revealed in a three-layer domain architecture that implies how GLRs transmit signals in plants (Green et al., 2021).In parallel, the expression patterns of plant GLR genes in response to different environmental cues such as abiotic stresses, pathogen infections and developmental transitions have recently been investigated.GLR genes have been shown to be essential for proper root growth and branching in Arabidopsis thaliana and rice, while plants with reduced GLR activity exhibit abnormal root growth (Li et al., 2006;Vincill et al., 2013;Singh et al., 2016).Plant GLR genes are also involved in plant responses to biotic and abiotic stresses.In particular, they have been shown to be vital for plant responses to salt, drought and chilling (cold) stresses (Lu et al., 2014;Cheng et al., 2016;Cheng et al., 2018;Zheng et al., 2021).In addition to their roles in abiotic stress tolerances, GLR genes also play an important role in plant defense systems against their pathogens.Cotton lines with CRISPR/Cas9-mediated GhGLR4.8 knockout become highly susceptible to the fungal pathogen, Fusarium oxysporum f. sp.vasinfectum (Fov) and lost the ability to induce the calcium influx in response to the secreted proteins of Fov (Liu et al., 2021).AtGLR3.3 and AtGLR3.6 were shown to trigger leaf-to-leaf systemic signaling and induce systemic defenses against its piercing-sucking insect, myzus persicae (green peach aphid) or its chewing insect, spodoptera litura (cotton leafworm) when under insect attack (Xue et al., 2022).These findings highlight the importance of GLR genes in sensing and responding to pest and disease erosions and suggest their potential uses in developing new strategies for enhancing plant resistance to pests and diseases.While the importance of GLR genes in plant growth/development, stress responses and pathogen defense are becoming increasingly clear, how different GLR genes are involved in these complex processes remains unclear.Thus, it remains an interesting subject regarding how these GLR receptors can be used as targets for improving crop production and sustainability.
As the most important cereal crop grown worldwide, rice is also known to be the model species of cereals for functional and population genomic research.Today, more than 4,600 rice genes have been cloned and functionally characterized at the molecular level (https://funricegenes.github.io/).Meanwhile, the 3010 diverse rice accessions (3KRG) from 89 countries worldwide have been sequenced, revealing extremely rich genomic diversity in the primary gene pool of this important species (Project, 2014;Wang et al., 2018;Zhang et al., 2021).Unfortunately, the rapid progresses in rice functional and population genomic research have not yet been widely applied to development of more efficient breeding technologies.This is because, three pieces of information are essential to most plant breeders, but generally lacking for virtually all cloned rice genes.First, the current information on phenotypic effect(s) of a cloned gene is incomplete because this information was obtained from typical molecular biology experiments using near isogenic lines in specific genetic backgrounds and environments and thus is irrelevant to target environments of most breeding programs.In other words, the genotype by environment interaction and genetic background effects on important agronomic traits are largely uncharacterized for most cloned rice genes.Secondly, although there is rich natural allelic variation at most gene loci in rice populations (Zhang et al., 2021), but it remains a huge challenge regarding how to determine and mine desirable allele(s) at the cloned rice gene loci from rice germplasm accessions for improving specific target trait(s).Third, the functionalities of most rice genes, such as most GLR genes, remain unknown.Thus, given the fact that the functionalities of most rice genes remain unknown, one important challenge is to obtain the required information on these uncharacterized genes for breeding application without going through the typical time-consuming gene cloning processes.
In this study, we employed an integrated approach that involved phylogenetic analysis of rice GLR genes alongside their counterparts in closely related plants, phenotypic and functional characterization of important rice GLR genes using CRISPR technology and transcriptomic analyses, as well as comprehensive population genetic analyses of their gene-CDS-haplotype (gcHap) diversity in rice populations.Collectively, we demonstrated a more efficient strategy to generate important information for large numbers of gene loci required for developing novel breeding technologies in future.

The phylogenetic relations between plant GLRs and mammals iGluRs (ionotropic glutamate receptors)
To understand the relationships between the rice GLR genes and those in animals and other plants, we compared all GLR genes in the Os-Nipponbare-Reference-IRGSP-1.0 reference genome with those in other flowering plants (Simon et al., 2023), including Arabidopsis thaliana (At), Zea mays (Zm), Solanum lycopersicum (Sl), Gossypium hirsutum (Gh), and Raphanus sativus (Rs), as well as basal land plants such as the moss Physcomitrium patens (Pp), the liverwort Marchantia polymorpha (Mp), and the lycophyte Selaginella moellendorffii (Sm), compared to the invertebrate Caenorhabditis elegans (Ce) and a-amino-3-hydroxy-5-methyl-4isoxazolepropionic acid receptors (AMPARs), N-methyl-daspartate receptors (NMDARs), kainate receptors (KARs), and dreceptors from mammals.Also included is the bacterial GluR0 from Synechocystis PCC 6803 and AvGluR1 from the freshwater rotifer Adineta vaga (Av).Sequences are available in the Supplementary Data 1.Phylogenetic analysis of all plant GLRs and mammals iGluRs was conducted using the MEGA7.0software (Shi et al., 2022).Muscle sequence comparison, which is a logarithmic expected multi-sequence comparison method, was performed to compare the protein sequences of GLR proteins of the relevant species.Then, a Bootstrap method with 1,000 replicates using the neighbor-joining (NJ) method (with p-distance correction) was used to assess the robustness of all nodes in the tree.The resulting phylogenetic tree was visualized using the website iTOL (https://itol.embl.de/),which allowed for the interactive exploration and customization of phylogenetic trees.

The gene-CDS-haplotype (gcHap) diversity analyses
The construction method of gcHap dataset in the 3,010 rice genomes (3KRG, Wang et al., 2018) was referred to the Molecular Plant article of Zhang (Zhang et al., 2021).A total of 913,360 single nucleotide polymorphisms (SNPs) were identified within the coding sequence (CDS) regions of 42,497 of the 45,963 annotated genes in the Os-Nipponbare-Reference-IRGSP-1.0 reference genome.In the rice genome, no SNPs were found within the CDS regions of 3,466 genes, which were defined as house-keeping (HK) genes.Shannon's equitability (E H ) was used to evaluate the level of gcHap diversity at all rice GLR loci in specific rice populations or the whole species (Sheldon, 1969).Nei's genetic identity ((I Nei ) was used to measure the genetic similarity between two populations or individuals based on their gcHap frequencies at different GLR loci (Nei, 1972).The formula for calculating Shannon's equability (E H ) and Nei's genetic identity (I Nei ) were referred to the method used in this article by Zhang (Zhang et al., 2021).The GraphPad Prism 9 software was used for statistical analysis and picture plotting of the obtained data.

The gcHaps diversity of rice GLR genes in modern varieties (MVs) and landraces (LANs)
To understand how modern breeding during past decades has affected gcHap diversity of rice GLR genes, detailed information was gathered for the total of 3,010 rice accessions of 3KRG.Of these, 732 were identified as Xian (indica) landraces (LANs-Xian), 358 were identified as Geng (japonica) landraces (LANs-Geng), 328 were identified as modern Xian varieties (MVs-Xian), and 139 were identified as modern Geng varieties (MVs-Geng).The first step was to download the predominant gcHaps (with the highest frequency in 3KRG) of each GLR gene from the RFGB webpage (https://www.rmbreeding.cn/Index).The frequency shift of the predominant gcHap for each GLR gene between modern varieties (MVs) and landraces (LANs) was calculated based on an R script (available at website https://github.com/isaac-Tsang/gcHap_diversity_LANs_MVs-).Then, A comparison was made between the distribution of gcHaps in modern Xian and Geng varieties and their respective landraces.At the same time, the loss and emergence of gcHap diversity were also analyzed in modern Xian/Geng varieties.Finally, the above data was plotted using GraphPad Prism 9 software.

Expression profiles analysis of rice GLR genes under normal and abiotic stress conditions
Firstly, we downloaded a rice expression database (RED), a repository of gene expression profiles for rice (Xia et al., 2017).The database includes rice gene expression data derived from RNA-Seq analysis of eight different tissues (aleurone, anthers, callus, leaves, panicles, pistils, roots, seeds and shoot) at different growth stages, as well as under a variety of biotic and abiotic treatments, including ABA, JA, Cd, drought, and cold treatment conditions.The expression profiles of all GLR genes are available in the Supplementary Data 2. The expression patterns of all rice GLR genes from the dataset were analyzed using the heatmap.2function of the R language and the layout of the image was done in the Adobe I l l u s t r a t o r s o f t w a r e ( h t t p s : / / w w w .a d o b e .c o m / c n / products/illustrator.html).

Knockout and field phenotypic experiments of rice GLR genes
Based on their strong expression in more rice tissues and/or their strong responses to abiotic stresses, five rice GLR genes, OsGLR2.2,OsGLR9.8,OsGLR6.8,OsGLR4.1 and OsGLR7.1, were chosen for functional analysis using the CRISPR-Cas9 technology in three steps.First, we designed two target sequences and corresponding detection primers for each of the five GLR genes (Supplementary Data 3).Cloning and transformation of these GLR genes were commissioned by Wuhan Boyuan Biotechnology Co., Ltd, using a Geng variety, Nipponbare as a background material for genetic transformation.Screening of the transformants via PCR amplification and sequencing to determine those carrying the desired mutations for the five target GLR genes were described previously (Zeng et al., 2020).2×Spark Taq PCR Master Mix (with dye) (Shandong Sparkjade Biotechnology Co., Ltd.) is used for PCR amplification.The PCR products were sent to Sangon Biotech (Shanghai, China) for Sanger sequencing.DNASTAR Lasergene software was used for sequence analysis and alignment (https:// www.dnastar.com/software/lasergene/).
Field phenotypic experiments of these GLR gene knockout mutants were conducted in the short-day environment of Hainan (18°15' N and 109°30' E) from Nov. 25 of 2021 -April 10 of 2022 and in the long-day environment of Hefei (31°52' N and 117°17' E) from May 27 -Oct.7 of 2022.Each homozygous GLR knockout mutant was planted in a 3-row plot with 8 plants per row at a spacing of 15 cm between individual plants and 20 cm between rows, with three replications for each of the knockout mutant lines and Nipponbare as the wild type check in the field experiments.Four agronomic traits, 1000-grain weight (TGW), seed setting rate (SR), panicle number (PN), panicle length (PL), and plant height (PH) were measured on the five middle plants in each plot.Statistical analyses of the phenotypic data of the experiments were performed using the GraphPad Prism 9 software, while the layout of the obtained picture plotting data was done using the Adobe Illustrator software (https://www.adobe.com/cn/products/illustrator.html).

Construction of gene CDS haplotype (gcHap) networks of rice GLR genes and association analyses of the major gcHaps at the GLR loci with yield traits in the 3KRG accessions
The construction of the gcHap networks containing major gcHaps with frequencies ≥1% in the 3KRG accessions or specific populations for each GLR gene was performed using the R package pegas (Paradis, 2010).The network of gcHaps for each GLR gene was produced using statistical parsimony algorithm, which is a method that connects the most closely related haplotypes first via the smallest number of mutations (Templeton et al., 1992).For more detailed steps, please refer to this Molecular Plant article of Zhang (Zhang et al., 2021).Finally, the layout of the image is done in the Adobe Illustrator software (https://www.adobe.com/cn/products/illustrator.html).

Results
The phylogenetic relations between plant GLRs and mammals iGluRs (ionotropic glutamate receptors) We identified 26 GLR genes in the Os-Nipponbare-Reference-IRGSP-1.0 reference genome clustered primarily on chromosomes 2, 6, and 9. To understand the evolution and functional diversification of GLR genes, we performed a phylogenetic analysis to compare the rice GLR genes with their corresponding ones from various flowering plants, including Arabidopsis thaliana (At), Zea mays (Zm), Solanum lycopersicum (Sl), Gossypium hirsutum (Gh), and Raphanus sativus (Rs), three basal land plants, moss Physcomitrium patens (Pp), liverwort Marchantia polymorpha (Mp), and lycophyte Selaginella moellendorffii (Sm), one invertebrate Caenorhabditis elegans (Ce) and four mammals, AMPARs, NMDARs, KARs, and d-receptors (Figure 1).Obviously, there is a clear distinction between GLRs in plants and animals, indicating a distant genetic relationship between the two groups.This is not unexpected as plants and animals represent distinct evolutionary lineages that diverged from a common ancestor more than a billion years ago.Consistent with our previous results (Shi et al., 2022), the 26 GLR genes could be roughly classified into four groups (Figure 1).

Functional differentiation and tissue specificity of rice GLR genes
To understand the functionalities of different rice GLR genes, we analyzed the expression levels of all GLR genes in eight different tissues of rice, including aleurone, anther, callus, leaf, panicle, pistil, root, seed and shoot tissues.Under the normal conditions, different GLR genes showed varied expression levels in different tissues of rice (Figure 2).Notably, OsGLR2.1,OsGLR4.1,OsGLR6.5, OsGLR6.8,OsGLR6.10,OsGLR6.12,OsGLR7.1,OsGLR7.2, and OsGLR9.8 were each expressed strongly (mean FPKM > 5) in more than one tissue, indicating their importance in rice growth and development.Among these genes, OsGLR4.1 had the highest expression levels in leaf and seed tissues (Figure 2).OsGLR2.2,OsGLR6.3 and OsGLR6.9showed moderate levels of expression in one or more tissues, while the remaining ones (OsGLR6.1,OsGLR6.2,OsGLR6.4,OsGLR6.6,OsGLR6.9,OsGLR6.11,OsGLR9.7 and OsGLR9.9)showed low levels of expression (<1) in all tissues.Interestingly, four GLRs (OsGLR9.3-OsGLR9.6)showed a moderate expression specifically in the pistil tissue but extremely low expression in other tissues.These results indicated that one important aspect in the functional differentiation of rice GLR genes was reflected in their tissue specificity, though it remains a mystery regarding how the expression pattern of each rice GLR gene related its functions in its specifically expressed tissues.

Expression profiles of rice GLR genes in response to abiotic stresses
In order to gain insights into the roles of rice GLR genes under different abiotic stress conditions, we conducted an integrative analysis of their expression profiles under abscisic acid (ABA), JA, cadmium (Cd), drought, and cold treatments.When treated with ABA, rice OsGLR6.5,OsGLR6.8,OsGLR6.10,OsGLR7.2 and OsGLR9.8showed increased expression in both shoots and roots (Figure 3A), indicating that they might be involved in the ABA regulated pathways.Interestingly, OsGLR4.1 and OsGLR9.8showed The phylogenetic relations between plant GLRs and mammals iGluRs (ionotropic glutamate receptors).The phylogenetic relations of glutamate receptors from Oryza sativa Geng (Os) based on the Os-Nipponbare-Reference-IRGSP-1.0 reference genome and other flowering plants that are described to show a conserved phenotype, including Arabidopsis thaliana (At), Zea mays (Zm), Solanum lycopersicum (Sl), Gossypium hirsutum (Gh), and Raphanus sativus (Rs), as well as basal land plants such as the moss Physcomitrium patens (Pp), the liverwort Marchantia polymorpha (Mp), and the lycophyte Selaginella moellendorffii (Sm), compared to the invertebrate Caenorhabditis elegans (Ce) and AMPARs, NMDARs, KARs, and dreceptors from mammals (without prefix).Also included is the bacterial GluR0 from Synechocystis PCC 6803 and AvGluR1 from the freshwater rotifer Adineta vaga (Av).OsGLRs are shown in colored text for clarity.Sequences were aligned using MUSCLE software, and the phylogenetic tree was constructed using the neighbor-joining method.Sequences are available in the Supplementary Data 1.AMPAR, a-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptor; KAR, kainate receptor; NMDAR, N-methyl-d-aspartate receptor; d, d-receptor.
ten rice GLR genes showed significantly increased expression in response to the five external treatments.Of these, OsGLR7.1,OsGLR2.1 and OsGLR9.8 were the most noteworthy ones.The expression of OsGLR7.1 in both roots and shoots and expression of OsGLR2.1 in roots were strongly upregulated in response to JA, Cd, drought and cold, suggesting their important roles in the JA signaling pathways involved in abiotic stresses.In contrast, the expression of OsGLR9.8 in shoots was strongly upregulated under the Cd, drought and cold stresses, indicating its important role in rice responses to abiotic stresses.

Functional determination of five rice GLR genes by knockout and phenotypic assessment
To assess the functional importance of rice GLR genes, five important OsGLR genes, OsGLR2.2,OsGLR9.8,OsGLR6.8,OsGLR4.1 and OsGLR7.1 were selected based on previous transcriptomic analyses for further gene knockout experiments using the CRISPR-Cas9 technology.Specifically, three independent homozygous mutants carrying different mutation sites for each of the GLR genes were obtained (Supplementary Figure 1) and assessed for 1000-grain weight (TGW), seed setting rate (SR), panicle number (PN), panicle length (PL) and plant height (PH) under the short-day and long-day environments of Hainan and Hefei.Three interesting observations were noted regarding the phenotypic differences for the measured traits between the GLR knockout mutants and the wild type (Figure 4).First, significant differences were observed between the knockout mutants of each GLR gene and the wild type for all the measured traits, but these differences varied greatly depending on the environments where these materials were assessed (Figures 4A-E).In other words, all five GLR genes showed huge genotype x environment interactions for their phenotypic effects on the measured traits, particularly for SR, PN and PH.When compared with the wild type, the mutant lines of all five GLR genes showed greatly reduced SR under the short-day environment of Hainan, but their SR was either significantly increased (OsGLR7.(OsGLR2.2) in the long-day environment of Hefei (Figure 4B).The mutant lines of all five GLR genes showed greatly increased PN under the short-day environment of Hainan, but not in the longday environment of Hefei (Figure 4C).More interestingly, the mutant lines of all five GLR genes showed significantly reduced PH under the short-day environment of Hainan, but significantly increased PH in the long-day environment of Hefei (Figure 4E).Secondly, one or more mutant lines of the GLR genes, except for OsGLR6.8 in Hainan, showed significant increased TGW in both environments (Figure 4A).Thirdly, the phenotypic differences between the mutant lines and wild type varied among different mutant lines of the same GLR genes, particularly for PL (Figure 4D), indicating that mutation sites in different GLR genes did make differences in generating phenotypically detectable mutants.These results led us to an important conclusion, i.e. obvious functional differentiation and functional redundancy regarding their phenotypic effects on the measured traits, the five rice GLR genes appeared to play an important role in contributing to the significant G x E interaction of agronomic traits in rice across the short-day and long-day environments.

The rich genetic and allelic diversity at rice GLR loci in rice populations
Given the functional importance of the rice GLR genes, it is important to understand their genetic diversity in rice populations in order to facilitate their applications in rice improvement.Using the gene CDS haplotype (gcHap) data from the 3,010 rice genome project (3KRG), we obtained the Shannon's equitability, E H , and number of gcHaps and major gcHaps (with frequency ≥1% in 3KRG) of the 26 GLR genes in four major rice populations (Supplementary Data 4, Figure 5).Of the 26 GLR genes, OsGLR6.2 and OsGLR6.3 are extremely conserved house-keeping (HK) genes with E H = 0, and gcHapN = 1 (Figure 5A), indicating all rice varieties share the same alleles at the two loci.For the remaining 24 polymorphic GLR genes, the average numbers of gcHaps (major gcHaps) and E H were 476.4 (8.4) and 0.368 in the 3KRG accessions.However, different GLR genes vary considerably in their genetic diversity.OsGLR9.4 has the highest genetic diversity with E H of 0.821 and 2,182 gcHaps (Figures 5C, D), followed by OsGLR9.2 with E H of 0.691 and 1,681 gcHaps, and OsGLR9.5 with E H of 0.649 and 1,374 gcHaps.OsGLR6.1 has the lowest diversity with E H of 0.050 and six gcHaps (Supplementary Data 4).
Different rice populations also vary considerably in their diversity of the GLR genes.The mean E H of the 24 GLR genes were 0.367 in population Xian, 0.213 in population Geng, 0.409 in population Aus, and 0.352 in population Bas.The numbers of detected gcHaps and major gcHaps were 342.1 and 7.7 in population Xian, 76.2 and 4.8 in population Geng, 47.8 and 5.2 in population Aus, and 17.8 and 5.0 in population Bas (Supplementary Data 4).Clearly, the numbers of gaHaps and major gcHaps in different rice populations were largely determined by their population sizes, implying that most gcHaps are rare ones with very low frequencies in different rice populations (Figures 5C, D).
To understand how different GLR genes were differentiated among major rice populations, we estimated the Nei's genetic identity (I Nei ) of all GLR genes using the gcHap data of the 24 polymorphic OsGLR genes, except for OsGLR6.2,OsGLR6.3, between all pairwise populations.Except for OsGLR6.1,OsGLR6.5 and OsGLR6.11,strong Xian-Geng sub-specific differentiation (I Nei <0.35) was observed at 21 of the GLR genes, 20, 17 and 18 of which also showed strong Xian-Bas, Aus-Geng and Aus-Bas differentiation (Figure 6, Supplementary Data 5), indicating that the allelic differentiation at most GLR loci had contributed very strongly to the differentiation of major rice populations and expectedly to the adaptations of different rice populations to their environments.
Additionally, the proportion of newly emerged gcHaps observed in modern varieties is significantly lower than that of lost gcHaps (Figures 7B, C).MVs-Xian gained average 62.8 new gcHaps/locus and 2.2 new major gcHaps that were absent in the LANs-Xian.The new gcHaps were apparently generated by intragenic recombination during breeding, while all new major gcHaps were rare ones in the LANs-Xian but became major ones with significantly increased frequencies in the MVs-Xian (Table 1).We noted four GLR loci (OsGLR9.2-OsGLR9.5) in a large cluster on chromosome 9 (Figure 7A), each of which had a huge number of gcHaps and showed more dramatically changes in gcHapN, suggesting genetic drift and unequal crossing-over occurred more frequently among duplicated GLR genes in the GLR gene clusters.The increased diversity and reduced gcHapN observed at many of the GLR loci suggested that significant frequency shifts in major gcHaps had occurred at these GLR loci in population Xian during breeding.Indeed, significant shifts in frequencies of the predominant gcHaps, F (P) , were observed at 13 of the 24 GLR loci, including significantly reduced F (P) at nine of the GLR loci and significantly increased F (P) at four of the GLR loci (Supplementary Data 7).Notably, at OsGLR7.1 locus, the predominant allele (Hap1) had an F (P) of 0.60 and 0.34 in populations LANs-Xian and MVs-Xian, or a major reduction in F (P) by 0.25 in MVs-Xian (Supplementary Figure 3), accompanied by significantly increased F (P) of five newly major gcHaps, all of which were rare gcHaps with frequencies <0.01 in the LANs-Xian but became major gcHaps in the MVs-Xian.Similar observations were noted on four other loci (OsGLR6.11,OsGLR9.1,OsGLR6.5 and OsGLR6.4)where significantly reduced F (P) was accompanied by emergences of 1-2 new major gcHaps in the MVs-Xian from rare ones in the in the LANs-Xian.The remaining 11 loci showing no significant differences in F (P) between LANs-Xian and MVs-Xian included the four GLR loci (OsGLR9.2-OsGLR9.5)where no or few predominant gcHaps were present.
In population Geng, the average E H of the 24 GLR genes was similar in the MVs-Geng (0.262) to 0.257 of the LANs-Geng (Supplementary Data 6).However, significant reduced diversity in the MVs-Geng was observed at only two GLR loci (OsGLR6.4 and OsGLR6.6), while significantly increased diversity was observed at OsGLR2.2, OsGLR7.1 and OsGLR9.7 genes.Similarly, the MVs-Geng had an average of 20.6 gcHaps/locus, significantly less than 42.6 gcHaps/locus in the LANs-Geng (Figure 7; Table 1).In fact, the MVs-Geng lost an average 35.3 gcHaps/locus and 2.8 major gcHaps/locus that were present in the LANs-Geng but absent in the MVs-Geng, while the MVs-Geng gained an average 13.5 new gcHaps/locus and 3.0 new major gcHaps that were present in the MVs-Geng but absent in the LANs-Geng (Table 1).The new gcHaps were apparently generated by intra-genic recombination during breeding, while the new major gcHaps were all come from minor rare gcHaps in the population LANs-Geng.Significant shifts in frequencies of the predominant gcHaps, F (P) , were observed at 12 of the 24 OsGLR loci, including significantly reduced F (P) at OsGLR6.7,OsGLR6.9,OsGLR7.1 and OsGLR9.7 and significantly increased F (P) at OsGLR6.4,OsGLR6.5, OsGLR6.6,OsGLR6.8,OsGLR6.10,OsGLR9.2,OsGLR9.5 and OsGLR7.2(Supplementary Data 7).It should be pointed out that populations Xian and Geng have different predominant gcHaps at all polymorphic GLR loci, except for OsGLR6.1 and OsGLR9.4.This was because there is only a single predominant in both populations Xian and Geng at OsGLR6.1 while there is absence of any gcHaps at OsGLR9.4 whose frequency was higher than 1% in population Xian.

Associations of the major gcHaps at rice GLR genes with important agronomic traits
To provide additional evidence for the functional importance of the 24 polymorphic GLR genes, we constructed the gcHap networks of the major alleles at 24 GLR genes in the four major rice populations, and analyzed the associations between the major gcHaps at each of the 24 GLR gene loci with phenotypic values of four agronomic traits, plant height (PH), panicle length (PL), panicle number per plant (PN) and thousand grain weight (TGW) in 3KRG (Figure 8; Supplementary Figures 4-7).Strong (P < 10 -7 ) associations were observed in 61 (63.5%) of the 96 (24 x 4) cases and major alleles at many of the GLR gene were strongly associated with trait values of one or more traits.Notably, major alleles at OsGLR2.2, OsGLR4.1,OsGLR7.1,OsGLR6.8 and OsGLR9.8 were strongly associated with trait values of all four traits, except for PL in two cases (Figure 8), consistent with the results from above knock-out experiments (Figure 4).In particular, major alleles at ten GLR loci (OsGLR6.4,OsGLR6.5,OsGLR6.6,OsGLR6.7,OsGLR6.8,OsGLR6.9,OsGLR6.10,OsGLR6.12OsGLR7.1 and OsGLR7.2) were strongly associated with trait values of all four traits (Figure 8; Supplementary Figures 4-7).In contrast, no associations were detected for major alleles at OsGLR6.1 and OsGLR9.1 with the trait values of the four traits, largely because of too few major gcHaps present at the three loci (Supplementary Data 4).The Nei's genetic identity (I Nei ) of rice OsGLR genes between all pairwise populations calculated from the gcHap data.Comparison of trait values between the predominant and 'unfavorable' gcHaps of GLR genes in rice In contrast to the predominant gcHap (with highest frequency) at a GLR locus in a rice population would be considered to be favored by natural selection during evolution, the major gcHap with the lowest frequencies in the population would more likely to be the 'unfavorable' one.To test if this scenario was holding true, we compared the phenotypic differences between the predominant gcHap and 'unfavorable' one at each of the 24 GLR loci for 15 agronomic traits.Of the total 24 x 15 = 360 comparisons, significant differences were detected in 149 (41.4%) cases between the predominant gcHaps and unfavorable at all GLR gene loci except for OsGLR6.8,OsGLR9.7 and OsGLR9.8(Supplementary Figure 8-31; Supplementary Data 8).Of the 15 traits, DTH was detected at 16 (66.7%) of the GLR loci, followed by PN at 14 (58.3%) of the GLR loci, consistent with the result of the gene knockout experiments that the five knockout GLR genes had large and consistent effects on PN and their phenotypic effects were influenced strongly by the two environments of contrasting day-lengths that are known to great influence on DTH of different rice varieties (Figure 4).Across the 24 GLR loci, significant differences in mean trait values between the predominant gcHap and unfavorable one were detected for almost all traits for OsGLR9.1 (14 traits), OsGLR9.9(14 traits), and OsGLR9.2(13 traits), indicating these three GLR genes were more likely affecting more traits than other GLR loci (Supplementary Data 8).
Mining 'favorable' allele(s) of rice GLR loci for improved productivity As described above, there were significant differences between or among the major gcHaps at most rice GLR loci regarding their effects on traits of agronomic importance.Thus, it is reasonable to assume the presence of 'favorable' gcHap(s) at most GLR loci for traits leading to increased productivity, even though the 'favorable' gcHap(s) may differ across different target environments.Thus, there should be one or more 'favorable' alleles, defined as the gcHap (s) associated with trait values for high productivity.Then, the frequencies of these most 'favorable' gcHaps in the MVs and different rice populations would of great interest for rice breeders.Figure 9 and Supplementary Figure 32 show the frequencies of the 'favorable' gcHaps at 24 rice GLR genes that affects PN, TGW, PH and PL in populations MVs-Xian and MVs-Geng and the four major rice populations of 3KRG.The general observation was that the 'favorable' allele frequency varied considerably across different yield traits, different GLR loci, and across populations MVs-Xian and MVs-Geng or different rice populations.For example, the Hap3 at OsGLR9.2 and Hap5 at OsGLR6.12 were associated with the highest PN in 3KRG accessions, which reached fixation in MVs-Xian, but had low frequencies in MVs-Geng (Figure 9).The same was true for eight GLR loci (OsGLR9.2-OsGLR9.6,OsGLR6.7,OsGLR6.8 and OsGLR6.12)with the highest TGW.In contrast, Hap1 at OsGLR2.2, Hap12 at OsGLR6.5, Hap6 at OsGLR6.9, and Hap5 at OsGLR6.10 were associated with the highest PN in 3KRG accessions, which reached fixation in MVs-Geng, but had low frequencies in MVs-Xian.This was consistent with the empirical observation in rice breeding that the extremely high values for any specific yield components may not generally result in the highest productivity.This result also suggested that the 'favorable' alleles for different yield traits could be very different across the two rice subspecies or genetic backgrounds and possibly across different environments.

Discussion
Although glutamate receptors are specialized proteins that perform a crucial role in signal transduction and communication in both animals and plants.Present as a large gene family in the rice genome, the molecular properties and functional characteristics of most rice GLR genes remain largely unknown.The phylogenetic analysis based on sequence similarity classified the 26 rice GLR genes into four major groups with groups II and III closely related The lost, new and unchanged gcHap variation of rice 26 OsGLR genes during modern breeding.to maize GLR genes, while group IV containing eight GLR genes clustered on rice chromosome 6 was unique ones only present rice only (Figure 1).In this study, we demonstrated an integrated strategy by limited experimental efforts plus comprehensive data mining and analyses to understand the functionalities of the rice GLR genes, focusing on the phenotypic characterization and genetic diversity of rice GLR genes because of their importance for potential utilization in rice improvement.In this regard, our results revealed several important properties of the rice GLR genes that led us to several conclusions on the functionalities and evolutionary history of rice GLR genes.
Our first conclusion was that the expansion of the rice GLR genes into a large gene family during evolution had gone through repeated gene duplication events evidenced by the presence of two large GLR gene clusters on rice chromosomes 9 and 6, and their greater gcHap diversity and rapid evolutionary rates (Table 1), and this expansion of the rice GLR genes had been accompanied with considerable functional differentiation, evidenced by their different tissue expression specificities (Figure 2), by their differentiated responses to different combinations of abiotic stresses and hormone treatments (Figure 3), and by their associations with different combinations of agronomic traits (Supplementary Figures 8-31).
Secondly, rice GLR genes may have played an important roles in rice sub-specific and population differentiation during evolution, evidenced by their very strong sub-specific and population differentiation at most rice GLR loci (Figure 6), by the great differentiated responses to day-length, the major factor determining rice subspecific adaptations to their environments (Figure 4), and by their differentiated responses to different abiotic stresses (Figure 3).In fact, the observed population differentiation at most GLR loci were far stronger than most of the 3,000+ cloned rice genes or genomewide SNPs analyzed previously (Wang et al., 2018;Zhang et al., 2021).
Thirdly, many rice GLR genes may have played an important role in rice improvement, suggested by the strong and large phenotypic effects of the five GLR gene knockout mutants on rice yield traits (Figure 4), by the significant shifts of the predominant gcHaps and diversity changes at many rice GLR loci (Figures 5-7, Table 1), and by their strong associations of the major gcHaps with Frequencies of the favorable gcHaps of all rice GLR genes affecting important agronomic traits (PN and TGW) in Xian/indica (XI) and Geng/japonica (GJ) modern varieties and different rice subpopulations.The favorable gcHap of a gene is defined as the gcHap associated with the highest trait value."#accession" indicates the number of accessions that possess the favorable gcHap.Five subpopulations of XI (XI-1A, XI-1B, XI-2, XI-3, and XIadm) and four subpopulations of GJ (temperate GJ: GJ-tmp, subtropical GJ: GJ-sbtrp, tropical GJ: GJ-trp, and GJ-adm) were classified by Wang et al. (Wang et al., 2018).agronomic traits at many GLR loci (Figures 8-9; Supplementary data 8).This is not surprising since acting as signal receptors, GLR genes are expected to play important regulatory roles involved in many processes of rice growth and development and their responses to environmental cues.However, it remains to be elucidated in future how different rice GLR genes are involved in these important processes at the molecular level.
Fourthly, the observed diversity increases at some of the GLR loci during breeding were somehow unexpected since modern rice breeding is believed by many to cause reduced diversity in MVs.For example, the convergent phenotypic selection during breeding for similar phenotypes associated for high productivity in virtually all breeding programs, which include more compact plant type (erect leaves and smaller culm and leaf angles), large panicles (more), moderate sizes of TGW, high SR, more PN, and tolerance/resistance to the common abiotic/biotic stresses, etc.Also, the fact that only a very small portion of germplasm accessions maintained in genebanks worldwide have been utilized as parents in specific breeding programs would expected to result in severe genetic bottleneck and thus reduced genetic diversity at most gene loci.Indeed, we observed average losses of 39.3% (20.8%) of gcHaps (major gcHaps) at the 24 GLR loci during Xian breeding, which were accompanied with the emergences of 62.8 (18.4%) and new gcHaps per locus and 2.2 (28.6%) new major gcHaps/locus in MVs-Xian (Figure 7).The average losses of gcHaps (major gcHaps) at the 24 GLR loci reached were pronounced during Geng breeding, reaching 46.3% (58.3%), which were accompanied with an average emergences of 13.5 (17.7%) new gcHaps/locus and 3.0 (62.5%) new major gcHaps/locus during Geng breeding.In other words, the losses of gcHaps from the genetic bottleneck effects during modern breeding were more severe for MVs-Geng than MVs-Xian.This was expected since line hybridization is a common practice of modern breeding, which would result in greatly increased recombination as compared to the extremely low outcrossing rate in rice landraces.The much greater numbers of gcHap losses/gains observed in the two large GLR clusters on chromosomes 9 and 6 (Table 1) strongly suggest that more frequently occurred unequal crossing-over events in the large gene clusters would not only generate gene copy number variation but also new alleles of the duplicated GLR genes.More interestingly, significant reduced frequencies of the predominant gcHaps, losses of a few major gcHaps and the emergences of new major gcHaps in MVs-Xian and MVs-Geng from the rare ones in LANs-Xian and LANs-Geng at specific GLR loci would suggest cases that the predominant and/or major gcHaps favored by natural selection were against by the artificial selection during modern breeding.This also implies that many rare gcHaps might be of value in rice improvement.Taken together, modern breeding activities in the past had the greater impact on the GLR gene loci than the 3,000 + cloned genes analyzed previously (Zhang et al., 2021), implying again the greater importance of the GLR genes in rice improvement.
Finally, our results would suggest the potential values of the natural variation at most rice GLR loci for improving the productivity and tolerances to abiotic stresses.All rice GLR genes, except for OsGLR6.2,OsGLR6.3,OsGLR6.1 and OsGLR9.1, are potential values to be utilized in rice trait improvement.Specifically, OsGLR2.1,OsGLR2.2,OsGLR4.1,OsGLR6.12,OsGLR7.1 and OsGLR9.8 could potentially be used for improving rice tolerances to abiotic stresses, while others can potentially be used for improving rice yield traits.However, challenges remain in the identification of 'favorable' alleles at specific GLR loci for improving specific yield traits because of the greater G x E interaction of most rice GLR genes on yield traits observed in this study.Thus, additional efforts are needed to obtain the required information in order to facilitate their application to improving rice productivity in future.

FIGURE 2
FIGURE 2 FIGURE 3Integrative analysis of the rice OsGLR gene family expression profiles under ABA (A), Cd (B, C), JA (D), cold (E), and dry (F) treatment conditions.
FIGURE 8Haplotype networks of five GLR genes with CRSPR-9 mediated knocking-out and their associations with four agronomic traits in 3KRG.(A) Based on gcHap diversity and contribution to subspecific/population differentiation, the rice OsGLR2.2 is classified as the other gene.(B) The OsGLR4.1 is classified as a X-G divergent gene.(C) The OsGLR7.1 is classified as a X-G divergent gene.(D) The rice OsGLR6.8 is classified as the other gene.(E) The OsGLR9.8 is classified as a X-G divergent gene.Within each haplotype network, two adjacent gcHaps are separated by mutational changes with hatching line indicating differences between the two most related haplotypes.The right side of each gene haplotype network corresponds to the phenotypic variation among the haplotypes.Boxplots are shown for the following four traits: plant height (PH), panicle length (PL), panicle number (PN) and 1000-grain weight (TGW).The P values under trait names indicate differences between the haplotypes assessed by two-way ANOVA, with different letters on the boxplots indicating statistically significant differences at P < 0.05 based on Duncan's multiple range test.The bar charts on the right show the differences in frequency of the predominant gcHaps between landraces (LANs) and modern varieties (MVs) in Xian and Geng.Chisquare tests were used to determine significant differences in the proportions of the same gcHap between the different populations with ****P < 0.0001, ***P < 0.0001, **P < 0.01, *P < 0.05, and n.s., not significant.

TABLE 1
Comparisons in genetic diversity of the 26 GLR genes between the landraces and modern varieties.

TABLE 1 Continued
H and gcHapN (major gcHapN) are the Shannon's equitability, and number of gcHaps (≥ 1% varieties) identified.The bold GLR genes were those used in functional analyses by the CRISPR-Case9 knockout experiments. E