Transcriptome Analysis of Arabidopsis GCR1 Mutant Reveals Its Roles in Stress, Hormones, Secondary Metabolism and Phosphate Starvation

The controversy over the existence or the need for G-protein coupled receptors (GPCRs) in plant G-protein signalling has overshadowed a more fundamental quest for the role of AtGCR1, the most studied and often considered the best candidate for GPCR in plants. Our whole transcriptome microarray analysis of the GCR1-knock-out mutant (gcr1-5) in Arabidopsis thaliana revealed 350 differentially expressed genes spanning all chromosomes. Many of them were hitherto unknown in the context of GCR1 or G-protein signalling, such as in phosphate starvation, storage compound and fatty acid biosynthesis, cell fate, etc. We also found some GCR1-responsive genes/processes that are reported to be regulated by heterotrimeric G-proteins, such as biotic and abiotic stress, hormone response and secondary metabolism. Thus, GCR1 could have G-protein-mediated as well as independent roles and regardless of whether it works as a GPCR, further analysis of the organism-wide role of GCR1 has a significance of its own.

Introduction G-protein signalling pathways are implicated in a variety of plant processes, but their upstream receptors and the role of typical G-protein coupled receptors (GPCRs) in plants are not very well understood [1,2]. The canonical GPCR gene in Arabidopsis (GCR1) was isolated independently by two groups [3,4]. The latter group initially called it a cytokinin receptor but subsequently found that cytokinin response was due to an independent mutation [5]. Others implicated GCR1 in PIPLC mediated regulation of DNA synthesis [6], abolishing seed dormancy and reducing flowering time [7]. Later, knock-out mutants of GCR1 were used to implicate it in seed germination in response to brassinosteroids and gibberellins [8], reducing
Seeds corresponding to the identified super-pools and sub-pools were grown and DNA was extracted from leaves for genotyping and sequencing of the mutant lesions. The primers KK66 and KK68 were used to detect the GCR1 wild-type copy, and primers KK66 and JL202 were used to detect the presence of the T-DNA in the GCR1 gene. PCR conditions were as above. PCR products were separated on agarose gels, and individual segregating plants for the GCR1 mutant (gcr1-5) were genotyped based on the presence or absence of wild-type and T-DNA bands.

Phenotypic characterization of the mutants
Seeds of the wild type, Ws2 and GCR1 mutant (gcr1-5) were surface sterilized using 70% ethanol and washed thrice with sterile ultrapure water and stratified at 4°C for two days on halfstrength B5 plates. The plates were then kept in growth chamber maintained at 22±1°C with a light intensity of 150 μM sec -1 m -2 and a photoperiod of 16:8 (light:dark). 10 day old plantlets were then transferred to 3.5 cm pots containing 1:1 mixture of soilrite and vermiculite. The pots were watered using sub-irrigation. The plants were allowed to grow for full life cycle and various phenotypic characters were measured at appropriate durations.

Plant material and RNA isolation
Arabidopsis thaliana GCR1 mutant (gcr1-5) as well as the corresponding wild type, were grown on 1X B5 medium hydroponically in a growth chamber at 22±1°C with a light intensity of 150 μM sec -1 m -2 and a photoperiod of 16:8 hours of light:dark cycle. The seeds were vernalized prior to inoculation at 4°C for 2-3 days. Total RNA was isolated from 3-4 week old whole seedlings using a modified hot phenol-LiCl method optimized on Arabidopsis [29]and reported in the context of cyanobacteria [30]. RNA samples were analyzed by Nanodrop spectrophotometer and Bioanalyzer (Agilent technologies, Santa Clara, USA) to determine the quality, quantity and suitability for microarray. RNA samples having a RIN (RNA Integrity Number) value greater than 6.0 were used for microarray experiments. The isolated RNAs were also used for validating the mutants using qPCR with gene-specific primers.

Microarray experiments and data processing
Microarray experiments were performed using Agilent 8×60k Arabidopsis array (AMADID 037661) with independent biological duplicates both the wild type, Ws2 and GCR1 mutant. Total RNA was transcribed into Cy3 labelled cRNA using Agilent Quick-Amp labelling kit as per manufacturer's instructions. Labelled cRNA was purified using RNeasy minikit (Qiagen) and the specific activity of cRNA was determined as a quality control for all the samples. They were hybridized with the microarrays using Agilent in-situ hybridization kit as per manufacturer's instructions. The washed slides were scanned and the images were manually verified to ensure that they are devoid of uneven hybridization, streaks, blobs and other artifacts. Hybridization across the slide was analyzed based on the number of features that were positive and significantly above background, i.e g(r) is PosAndSignif. Overall the microarray images were clean, had uniform intensity and with very low background noise. The data was then extracted from images by using Feature Extraction 10.7 software (Agilent Technologies).

Data analysis
The data was then normalized using the recommended 'Per Chip and Per Gene Normalization' feature of the software GeneSpring GX Version 11.5. The correlation of replicates was checked using principal component analysis and correlation coefficients were obtained. The geometric mean (geomean) fold change values are represented as log2. The average data of biological replicates was taken for final calculations. Log2fold change value of 1.0 with p-value of 0.05 was taken cut-off for differential-regulation.

Functional classification of significant genes
The differentially regulated gene lists were assigned gene ontology terms according the Arabidopsis Information Resource (TAIR 10) [31]. The differentially regulated gene lists were subjected to enriched GO categorization using AgriGO with default settings. Pathway analysis of the DEGs to obtain the list of changed pathways was done using plant MetGenMAP, which takes AraCyc as the background. Further functional classification was also carried out using Mapman tool, where the DEGs were assigned to different biological processes (bins). This tool also takes into account the log2fold change and represents it as coloured boxes on the software generated biological process map.

Data validation using qPCR
Differentially expressed genes obtained from microarray analyses were verified by RT-qPCR using Stratagene Mx3000P (Agilent technologies). Typically, total RNA was digested by RNase-free DNase (Fermantas), repurified, quantified and 5 μg of RNA was used for cDNA preparation for each biological replicate using Oligo(dT) primers and RevertAid reverse transcriptase (Fermentas). The analyses were done using biological triplicates, out of which two were the same as used for microarray. Sequences for designing the primers were obtained from TAIR. PCR amplifications were performed in 20 μl by using the BrilliantIII Ultrafast SYBR Green QPCR mastermix (Agilent Technologies) with 1.0 μl of sample cDNA and 100 nmoles of each gene-specific primer. Primer efficiency was determined by serial dilution of the template and only primers that worked at 90-110% efficiency were used for all qPCR analyses. The specificity of primer pairs was obtained by melting curve analysis of the amplicons. Actin2 (ACT2) was used as an internal control for normalization. Quantification of the relative changes in gene expression was performed by using Pffafl method [32].

GCR1 mutant characterization
Sequencing of the PCR-amplified GCR1 gene confirmed the mutation in the second intron of GCR1 (Fig. 1A), similar to the mutant gcr1-3 described previously, except that it had two T-DNAs with their right borders facing each other, instead of the single T-DNA insertion reported earlier [9]. qPCR with gene-specific primers showed no expression of GCR1 transcript, confirming that this is a knock-out mutation (Fig. 1B). The mutant plants (gcr1-5) were phenotypically characterized for root length, plant height, leaf shape, etc. It was found that other than total number of siliques (S1 Fig.), gcr1-5 is barely distinguishable from the wild type as is the case with other known GCR1 mutants [8,9].

Microarray analysis and validation
The experiment was carried out in accordance with MIAME compliance and the high correlation coefficients between biological replicates (>0.95) clearly indicate the robustness and a high level of reproducibility of the data (S1 Table). A stringent cut-off value of 1.0 (geometric mean log 2 ) with a p-value of 0.05 was used for determining the up-or down-regulated genes in the mutant with respect to the wild type control. A total of 432 differentially regulated spots were obtained in the mutant (314 up-regulated and 118 down-regulated). These spots corresponded to 350 unique differentially expressed genes (DEGs) in the mutant (265 up-regulated and 85 down-regulated). A list of the top 20 up-and down-regulated genes are shown in Table 1. A heatmap of all the differentially regulated genes is shown in Fig. 2a. A set of 17 (8 up-regulated and 9 down regulated) genes representing all the affected biological processes were subjected to RT-qPCR using gene specific primers checked for efficiency (90-100%). The results of RT-qPCR matched with the microarray data in all the cases with a Pearson product moment correlation of 0.98 (p-value > 0.0001) (Fig. 3), validating the basic trends of regulation of gene expression on the microarray. The list of candidate genes along with the primer sequences are given in the S2 Table. Stress response genes among GCR1-regulated genes Mapman analysis mapped 119 DEGs as belonging to the biotic stress response category (Fig. 4), including some genes that are involved in both biotic and abiotic stresses and other related caterogies such as hormone signalling, transcription factors etc., with a putative role in biotic stress. This finding was substantiated by further analysis using GO slim classification. It mapped 310 genes, out of which 101 genes (32.6%) belonged to the stress responsive category, the largest  among the 'known' processes ( Fig. 2b), even though many more genes are categorized under 'unknown' or 'other' processes. Singular enrichment analysis (SEA) using AgriGO, revealed the very high significance levels of the gene clusters (S3 Table). Enriched GO classification revealed that all the classes can be broadly grouped into only four categories, i.e. response to stimulus (mainly stress), protein modification, response to hormones and response to nutrient levels or starvation (Fig. 5). The largest enriched class was that of response to stimulus with 75 genes, which included several stress-related classes like defense response, response to stress, defense response-incompatible interaction and response to stimulus as the most significant classes in their decreasing order of significance level. In terms of the number of genes per cluster, the largest clusters belong to plant defensins (PDF), mildew resistance locus O proteins (MLO), and the Toll Interleukin1 receptor-nucleotide binding site-leucine-rich repeat type R (TIR-NBS-LRR) class, which are involved in biotic stress response. In addition, there are genes involved in abiotic stress such as, Arabidopsis Zinc Finger (AZF) and Arabidopsis thaliana Phloem Protein 2 A5 (ATPP2-A5). Several peroxidases and transcription factors involved in biotic and abiotic stress response were also found to be differentially expressed in the mutant.

Hormone biosynthesis/response genes affected in GCR1 mutant
AgriGO analysis also revealed that some of the DEGs in the GCR1 mutant are involved in hormone biosynthesis as well as in response to various hormones (S3 Table). They include cytokinin oxidase/dehyrogenase (CKX4) which is involved in cytokinin biosynthesis as well in abiotic stress response, genes involved in salicylic acid response (PDF1.2, PDR12, WRKY18, GLIP1, etc), ethylene response (ERS2, ERF13, PDR12, etc) and ABA mediated pathway (LECRKA4.2, LECRKA4.3). These results indicate the prominent role of GCR1 in hormone biosynthesis/response.

Phosphate starvation response genes among GCR1-regulated genes
Interestingly, enriched GO analysis revealed several genes related to phosphate response/starvation as downregulated in the GCR1 mutant ( Table 2, Fig. 5c). They include SPX domain containing genes (SPX1, SPX3) and galactolipid/sulfolipid biosynthesis genes (MGD2, MGDC, SQD2), which are known to get differentially regulated on phosphate starvation. This finding would need further investigation.

Protein kinases/phosphatases and transcription factors regulated by GCR1
Functional classification of DEGs identified by enriched GO analysis revealed many genes involved in protein modification by phosphorylation/dephosphorylation (S3 Table, Fig. 5d). The list included several serine/threonine kinases, leucine-rich kinases and S-locus pectin kinase families, and phosphatases like PP2C and PPCK, some of which are also annotated under stress response. Transcription and gene regulation also figured as major categories in functional classification. Comparing them with the Plant transcription factor database (PlantTFDB 2.0) revealed more than 25 transcription factor families, most of which were up-regulated in the mutant. The most represented transcription factor families were WRKY, MYB, AVI3/VP1, bHLH and C2H2 (Fig. 6). A large number of putative DNA binding and unspecified TFs also figured in the list. Secondary metabolite biosynthesis among GCR1-regulated pathways When the DEG list was mapped to AraCyc [33] using plant MetGenMAP [34], most of the pathways belonged to biosynthesis of secondary metabolites. These included pathways like sulpholipid biosynthesis, fatty acid α oxidation, camalexin biosynthesis, leucopelargonidin and leucocyanidin biosynthesis, leucodelphinidin biosynthesis, flavonoid biosynthesis, homogalacturonan degradation and glycolipid biosynthesis as significantly changed pathways (p-value 0.05) ( Table 2). Out of the above pathways, leucopelargonindin and leucocyanidin biosynthesis, leucodelphinidin biosynthesis and flavonoid biosynthesis are part of the superpathway of flavones and derivatives biosynthesis. Both the genes encoding flavanone 3-beta hydroxylase (AT5G20550, AT2G36690) were found to be up-regulated, whereas two genes involved in glycolipid biosynthesis (AT2G11810, AT5G20410) and one each in sulpholipid biosynthesis (AT5G01220) and camalexin biosynthesis (AT3G26830) were down-regulated. The DEGs involved in all other pathways were found to be up-regulated.

Discussion
The unsettled controversy over the existence or the need for G-protein coupled receptors (GPCRs) in plant G-protein signalling has overshadowed a more fundamental quest for the overall role of GCR1 in Arabidopsis. We approached it from the gene discovery perspective, using functional genomics to investigate the genomewide role of GCR1 in Arabidopsis, including, but not limited to its more debated role in G-protein-regulated processes. For this purpose, we isolated a T-DNA knock-out mutant of GCR1 (Fig. 1, S1 Fig.) disrupted in its 2 nd intron, which turned out to be similar, though not identical to gcr1-3 reported earlier [9] and also phenotypically similar to all other known GCR1 mutants [8]. It must be noted that the slightly longer, narrower leaves of this GCR1 mutant is at variance with the rounder leaves characteristic of GPA1 and AGB1 mutants that represent G-alpha and beta subunits respectively [2,35,36]. Microarray analysis of the whole plant transcriptomes of the WT and the GCR1 mutant revealed that a total of 350 genes were differentially regulated in the mutant when compared to the wild type (Fig. 2), using a stringent cut-off value of 1.0 (geometric mean log 2 ) with a p-value of 0.05. These genes span all five chromosomes, indicating an extensive genomewide role for GCR1. A few of the upregulated and downregulated genes have been validated by RT-qPCR (Fig. 3, S2 Table) and the top 20 up/down-regulated genes are shown in Table 1. GOslim classification of differentially expressed genes (DEGs) revealed that over 50% of them are involved in metabolic processes, while 32.5% are responsive to stress, 19% to biotic and abiotic stimuli, 16.5% involved in transport and 15.4% in transcription, with several genes figuring in more than one category. A closer examination revealed that some of the GCR1-responsive genes/processes have close parallels to those attributed to G-protein signalling in Arabidopsis, whereas others have either novel roles or hither to unknown in the context of G-protein signalling, as detailed below. Novel roles of GCR1 similar to those attributed to G-protein signalling Stress emerged as one of the major categories and enriched GO as well as Mapman analyses revealed that most of these differentially regulated genes pertained to biotic stress, a novel finding hitherto unknown in relation to GCR1 or any other predicted GPCR. They include genes like PDF1.2 and MLO12, which are known to be involved in innate immune response and defense against pathogens [37,38]. We also found TIR-NBS-LRR class of proteins, which have been implicated in pathogen sensing and defense response in plants [39], as well as many receptor like kinases known to be involved in plant defense signalling. Our findings on the role of GCR1 in biotic stress in this study have some close parallels with reports that implicate heterotrimeric G-proteins in biotic stress response [40,41] and plant defense signalling [42][43][44]. One of the reports has shown that GCR1 partners with GPA1 in the regulation of root growth mediated by bacterial quorum sensing signals [45]. At least two genes identified as GCR1 responsive in our study (PDF1.2 and PAD3), have also been reported earlier in G-protein-mediated response to biotic stress and jasmonic acid [46,47]. The role of GCR1 in abiotic stress response is by no means insignificant, even though fewer genes figure in this category than in biotic stress. Our transcriptome data corroborate the role of GCR1 and/or G-proteins in several abiotic stresses reported in separate studies, such as drought [9], heat, cold, salt [48,49], oxidative stress [50].
Secondary metabolic pathways such as camalexin and flavonoid biosynthesis constitute another novel category of GCR1-responsive processes that was previously unknown in relation to plant GPCRs. In fungi, G-protein alpha subunit was recently shown to be involved in secondary metabolism/biosynthesis of secondary metabolites [51]. Therefore, it remains to be seen whether plant heterotrimeric G-proteins mediate biosynthesis of secondary metabolites and whether they partner with GCR1.
GCR1 is involved in hormonal responses, as is evident from the categorization of the DEGs. They are involved in the biosynthesis and mediation of multiple hormonal responses like cytokinin, ethylene, salicylic acid, ABA, etc. Some of them, including ABC transporter (PDR12) and transcription factors (WRKY18) mediate hormone responsive effects that are also involved in the regulation of stress response. Heterotrimeric G-proteins in general and GPA1 and GCR1 in particular have been implicated earlier in a variety of hormone responses [8,9,52,53]. Our transcriptome data not only corroborate them category-wise, but also in terms of some of the genes involved in these processes. At least 20 of them were earlier reported in G-protein-mediated ABA response [54].
G-proteins are known to be involved in the regulation of guard cell functions and root hair differentiation [9,52], as are transcription factor complexes comprised of bHLH, MYB and WD40 domains [55]. Our finding that bHLH and MYB are GCR1 responsive raises the possibility that these TF complexes mediate the regulation of guard cell and root hair differentiation via GCR1 and G-proteins, which needs experimental validation.
Further evidence regarding the upstream GCR1 ligand and downstream GEF activity in all the above processes would be needed before it can be conclusively established that GCR1 and GPA1 act together in those processes/responses in Arabidopsis.
Novel roles of GCR1 unknown or unrelated to plant G-protein signalling MYB and WRKY class of transcription factors are known to be involved in metabolic regulation, stress and cell fate [56,57]. Our data show that the above transcription factors, as well as the above processes they regulate are GCR1 responsive, is a significant and novel observation. This is because, barring stress discussed above, the role of GPCR or G-protein has not been reported in any of them. This is also true of GCR1-regulation of AP2/EREB class of TFs found in our study, which are known to be involved in storage compound and fatty acid biosynthesis [58]. If confirmed using other G-protein mutants, it could mean that GCR1-regulation of these genes/processes works independent of GPCR/G-protein signalling.
Another novel finding of this study is that GCR1 downregulates the genes involved in response to phosphate starvation. These genes include SPX domain containing proteins (SPX1 and SPX3) and IPS1 (induced by phosphate starvation1) [59]; MGD2 and MGDC, which are involved in galactolipid metabolism during phosphate starvation [60]. Till date, no study has implicated phosphate starvation to any component of G-proteins.

Conclusions
Our transcriptome analysis shows for the first time that a) GCR1 has an extensive genomewide response, revealing many hitherto unknown genes and processes including, but not limited to those regulated by the known G-proteins in Arabidopsis; b) some of its roles, such as in biotic and abiotic stress, hormone response and secondary metabolism are among those regulated by the known heterotrimeric G-proteins; and c) GCR1 also has other important roles that are either independent of, or hitherto unattributed to G-protein signalling, such as in phosphate starvation, storage compound and fatty acid biosynthesis, cell fate etc., though the possibility of some of them being regulated by other yet-to-be identified G-proteins or their components cannot be ruled out. Overall, our results point to a serious need to revisit the role of GCR1 in G-protein signalling in Arabidopsis, including its possible role as a GPCR. Further, considering that GCR1 and Arabidopsis became the main basis for contesting the existence and the role of GPCRs in plants, our results also suggest that it may be too early to write-off plant GPCRs, or their role in plant G-protein signaling.  Table. Correlation coefficients of microarray spot intensities across biological replicates of the wild type, Ws2 and the GCR1 mutant. (DOC) S2 Table. List of genes used for validation of microarray data, with their primer sequences and efficiencies. (DOC) S3 Table. Singular enrichment analysis (SEA) using AgriGO to obtain list of all GO biological processes with their significance levels.

Supporting Information
(DOC)