GABA promotes peroxisome proliferation in Triticum monococcum leaves

Abstract

Although peroxisomes are integral for both primary and secondary metabolism, how developmental changes affect activity of peroxisomes remains poorly understood. Here, we used published RNA‐seq data to analyze the expression patterns of genes encoding 21 peroxisome metabolic pathways at successive developmental stages of Zea mays and Oryza sativa. Photorespiration was the most represented pathway in adult leaf relative to the juvenile stages. Components of reactive oxygen species (ROS)/reactive nitrogen species (RNS) metabolism, NADPH regeneration, and catabolism of polyamines were also enriched at later stages of leaf differentiation. The most commonly upregulated gene in differentiated leaves across all datasets of both species was BETAINE ALANINE DEHYDROGENASE (BADH). BADH functions in catabolism of polyamines where it converts 4‐aminobutyraldehyde (ABAL) to 4‐aminobutyrate (GABA). We tested the outcome of RNA‐seq analysis by qRT‐PCR in developing Triticum monococcum ssp. monococcum (Einkorn) seedlings. Consistent with the outcomes of RNA‐seq analysis, transcription of BADH and CATALASE3 (CAT3) were upregulated in older seedlings. CAT3 is an essential peroxisome biogenesis factor and a key enzyme of ROS homeostasis. Furthermore, exogenous application of GABA resulted in higher peroxisome abundance and transcriptional upregulation of BADH and a gene encoding another peroxisome biogenesis factor responsible for peroxisome fission, PEROXIN11C (PEX11C), in leaves. We propose that GABA contributes to regulation of peroxisome fission machinery during leaf differentiation.

INTRODUCTION

Peroxisome is a single-membrane organelle involved in many essential biochemical pathways including photorespiration, reactive oxygen species (ROS) metabolism, β-oxidation of fatty acids, the acetate–malate shunt, phytohormone biosynthesis, NADPH regeneration, the mevalonic acid pathway, amino acid metabolism, and others. Recent progress in functional characterization of peroxisomal proteins and their role in plant biology was facilitated by molecular biology and bioinformatics approaches along with x-ray methods and advanced microscopy techniques (Anteghini et al., 2021; Pan et al., 2020; Reumann & Chowdhary, 2018). Enzymes such as catalase (CAT), superoxide dismutase (SOD), monodehydroascorbate reductase, GSH (glutathione) reductase, ascorbate peroxidase, and peroxiredoxin in peroxisomal lumen contribute to ROS scavenging under normal and stress conditions. Peroxisome fission is governed by PEROXINs (PEX), FISSION1 (FIS1) and DYNAMIN-RELATED PROTEIN (DRP) genes (Pan et al., 2020). Drought and heat stress increase hydrogen peroxide (H₂O₂) content in leaves, and this elevation is correlated with the peroxisome proliferation along with higher transcription activity of fission machinery (Hinojosa et al., 2019). Salt stress induces the transcription of peroxisome-associated genes under salt stress (Charlton et al., 2005). Wheat exhibits higher transcription of peroxisome fission genes under drought stress (Sanad et al., 2020). It has been proposed that peroxisome abundance could be used as a proxy of redox homeostasis during drought response (Smertenko, 2017).

Betaine alanine dehydrogenase (BADH) is a peroxisomal enzyme involved in glycine betaine (GB) synthesis, NADP⁺(H) regeneration, and polyamines catabolism (Jia et al., 2002; Li et al., 2003, 2021; Lv et al., 2016; Zhang et al., 2011). BADH activity contributes to cadmium tolerance as well as ROS scavenging (Li et al., 2021), and overexpression of BADH increases tolerance to drought stress (Zhao et al., 2021). BADH also catalyzes conversion of γ-aminobutyraldehyde to GABA (γ-aminobutyric acid). GABA modulates plant growth and development by acting in conjunction with the signaling pathways such as nitric oxide, hydrogen peroxide, calcium–calmodulin (Ca²⁺/CaM) complex, cyclic guanosine-3′,5′-monophosphate, and mitogen-activated protein kinase (Jiao et al., 2019; Suhel et al., 2022). GABA accumulates in response to salt and drought stress, and wounding in Arabidopsis thaliana (Li et al., 2021; Su et al., 2019; Xu et al., 2021), UV-B in Glycine max (Jiao et al., 2019), and advances chilling resilience in tomato and wheat seedlings (Malekzadeh et al., 2012, 2014). GABA upregulates expression of ROS scavenging genes and regulates hydrogen peroxide producing enzymes such as NADPH oxidase, peroxidase, and amine oxidase in Caragana intermedia roots under NaCl stress (Shi et al., 2010). Moreover, exogenous application of GABA significantly improves growth of seedlings by increasing photosynthesis, gas exchange, and chlorophyll biosynthesis (Li et al., 2016).

Characterization of peroxisome metabolic activities has been extensively studied during stress responses. However, peroxisome metabolism could also be regulated in the developmental context. To test this hypothesis, we examined expression patterns of peroxisomal genes using RNA-seq data obtained from developing leaves in two representatives of the Poaceae family Zea mays (maize) and Oryza sativa (rice). The Poaceae family consists of diverse grass species including several most important crops: maize, rice, and Triticum aestivum (bread wheat). Although these species share extensive homologies on the genome level (Ahn et al., 1993; Matsuoka et al., 2002), metabolic analysis of maize and rice revealed significant differences (Deng et al., 2020). Previous studies have also suggested that another member of the Poaceae family Triticum monococcum ssp. monococcum (Einkorn), a relic crop with diploid genome cultivated in Turkey, Spain, Morocco, and southern Europe, contains higher antioxidant and bioactive compounds such as phenolics, flavonoids, and tocopherols (vitamin E) than durum or bread wheat (Sahin et al., 2017; Serpen et al., 2008). We reason that members of this family represent a good model system to study metabolic changes of peroxisomes during development.

We found transcriptional changes of genes encoding many peroxisomal metabolic pathways of which photorespiration was the most represented pathway in adult leaf, relative to the juvenile stages. Components of ROS/reactive nitrogen species (RNS) metabolism, NADPH regeneration, and catabolism of polyamines were also enriched at later stages of leaf differentiation. Interestingly, genes encoding fatty acid and amino acid metabolism were less represented in adult leaves. Our results demonstrate transcriptional upregulation of BADH at successive developmental stages of seedlings. In addition, our findings suggest that GABA plays a role in regulating peroxisome proliferation as the exogenous application of GABA to Einkorn seedlings resulted in an increase in peroxisome abundance. The transcriptional upregulation of PEX11C pointed out that GABA could contribute to regulation of peroxisomal abundance by promoting peroxisome fission.

MATERIALS AND METHODS

High throughput RNA-seq datasets acquisition

RNA-seq datasets were downloaded from the Gene Expression Omnibus (GEO, , accessed on October 1, 2021). Experiments were identified using keywords for development stages and drought stress, “development” AND “species name” [organism]. Four considerations were followed for selecting the RNA-Seq data, and the workflow of the study was depicted in Figure S1.

Data preprocessing

The read quality of raw sequence data was assessed using FastQC software ver. .11.9 (Andrews, 2010). After the quality check, the raw sequences that showed low quality were subjected to the read trimming tool. The adapter sequence, low-quality bases, and short reads were trimmed using Trimmomatic ver. .39 (Bolger et al., 2014). FASTQ files were screened to the level of Q30 and length >50 bases. Hisat2 ver. 2.1.0 was used as the alignment software (Pertea et al., 2016).

Z. mays reads were mapped to the B73v4 genome (); O. sativa reads were mapped to the MSU7 genome (). Counts of the transcript were obtained by the featureCounts software ver. 1.6.4 (Liao et al., 2014). The output of featureCounts, which contains actual read counts per gene (with gene ID, genomic coordinates of the gene including strand and length), was used in the differentially expressed gene analysis (DEGs) in the R programming environment.

Principal component analysis (PCA)

Global read count patterns were determined by PCA. 2D plots using PC1 and PC2 were constructed with the built-in plotPCA function provided by the R DESeq2 package. First, DESeqDataSetFromMatrix() was created with raw counts, and then all samples were normalized, and logarithm transformed with vst() function. The genes that show greater variance were selected with the rowVars() function. The top 2000 genes were used to construct the 2-D PCA. The dendrograms were generated by the hclust() function.

Determination of DEGs using Deseq2 package

Deseq2 ver. 3.10 (R Bioconductor package) was used to identify DEGs between each pair of samples (Love et al., 2014). A DESeqDataSet object was generated with the DESeqDataSetFromMatrix() function. The genes that included less than 100 counts in the rowSums() function were eliminated. The DEG analysis was run with DESeq() function. The results were retrieved by the results function (Deseq2 object, independentFiltering = TRUE, alpha = .05, lfcThreshold = 0). Genes with at least a|log2-foldchange| > 0 in expression and a Benjamini-Hochberg adjusted p-value (q-value) <0.05 were considered as DEGs.

GO term enrichment analysis

GO term annotations of the genes were downloaded from , accessed on October 1, 2021. The topGO R Bioconductor package was used for the enrichment analysis of the DEGs via Fisher method (Alexa & Rahnenfuhrer, 2021). First, a factor vector was formed with the factor () function to create a new (“topGOdata,” ontology = “BP,” allGenes = factor_vector, nodeSize = 1, annot = annFUN.gene2GO, gene2GO = your_GO_table) object built-in topGO package. Fisher test was applied to data by the runTest function (object, algorithm = “classic,” statistic = “fisher”). False discovery rate of p-values was adjusted by the Benjamini-Hochberg false discovery rate method.

Determination of peroxisome gene homologs in species

Peroxisomal genes were determined according to Kaur and Hu (2011) in O. sativa and A. thaliana. Orthologous homologs of peroxisomal genes encoding peroxisomal proteins were determined in other species by using OrthoFinder software ver. 2.5.1 (Emms & Kelly, 2019). The list of genes is shown in Table S1. The homologs were categorized according to 21 metabolic pathways and processes in peroxisomes (Pan et al., 2020). Genes with a role in two or more processes were categorized into a new group. Heatmaps were generated with the pheatmap() function (NMF ver. .17.6) using vst values with a z-score transformation.

Cultivation of plants

T. monococcum ssp. monococcum seeds were kindly provided by the US Department of Agriculture, Agricultural Research Service. T. monococcum ssp. monococcum seeds were grown in a growth chamber under 60% humidity, a 16/8 h light/dark cycle, 22°C during the day and 18°C at night, and light intensity at the plant level of 400 𝜇mol quanta m⁻² s⁻¹. Seeds were planted into 5 × 5 × 5.25″ square pots with Sungro6 peat moss potting soil. All pots were watered daily.

Exogenous GABA treatment

Two centimeters long fragments of leaf basal part collected from 3-week old seedling of T. monococcum were submerged in GABA (Sigma, St Louis, MO, USA) solution at concentration 50 mgL⁻¹ for 0, 10, or 30 min. The fragments were flash-frozen in liquid nitrogen and used for measuring peroxisome abundance or isolation of RNA. Three replicates were performed per treatment with each replicate comprising leaf fragments from eight seedlings.

Measurement of peroxisome abundance

Peroxisome abundance was measured using small fluorescent probe Nitro-BODIPY using previously published procedure (Fahy et al., 2017). A 2-cm fragment of the leaf basal part was transferred into deep 96-well plates immersed in a liquid nitrogen bath and ground to fine powder using a tissue grinder (TissueLyser II, Qiagen, Venlo, Netherlands). Total leaf protein was extracted using 0.8 mL of the extraction buffer A (EBA; 20 mM Tris HCl, pH 7.4, 500 mM NaCl, 7 M Urea) by rotating the plates for 1 h. The debris was cleared by centrifugation at 4000g for 30 min. Then 20 μL of the extract was added to 100 μL of freshly prepared 2 mM solution of N-BODIPY and 80 μL of water in 96-well plates and incubated for 10 min. The fluorescence intensity was measured at 490 nm excitation wavelength and 530 nm emission wavelength using Synergy Neo B spectrofluorometer (Biotek Instrument, Inc.). Five biological replicates (individual plants) with three technical replicates were performed per genotype and treatment. The background was measured as (i) 20 μL of each protein extract in 180 μL of water and (ii) 100 μL of N-BODIPY solution supplemented with 80 μL of water and 20 μL of extraction buffer per each 96-well plate. Both background values were subtracted from the N-BODIPY fluorescence signal value. The fluorescence intensity was normalized by the protein concentration measured with the Bradford Reagent (Bio-Rad Laboratories) using a calibration curve constructed with solutions of known concentration of bovine serum albumin. Fluorescence intensity was calculated in arbitrary units per 1 mg of protein. Statistical analysis was performed using ANOVA test. Results were presented as mean ± SD. A p-value of less than .05 was considered to be statistically significant.

Quantitative reverse transcriptase PCR analysis

Total RNA was extracted from three biological replicates using RNeasy plant kit (Qiagen). cDNA was synthesized using Maxima H Minus First Strand cDNA Synthesis Kit (Thermo Fisher Scientific). Each qPCR primer was designed to target all three homeologs. The primers are listed in Table S2.

qRT-PCR reactions were performed using Fast SYBR™ GreenMaster Mix (Thermo Fisher Scientific, Waltham, MA) in 96-well plates on ViiA 7 Real-Time PCR System with default ViiA™ 7 SYBR conditions. Reactions were replicated three times and analyzed in QuantStudio™ Real-Time PCR Software v1.3; transcription levels were normalized to the housekeeping gene RNase L inhibitor-like protein (Giménez et al., 2011).

RESULTS AND DISCUSSION

Analysis of genes transcription during seedling development

The RNA-seq datasets of developing Z. mays and O. sativa leaves were selected using the following four criteria: (1) The experimental design included at least two biological replicates; (2) RNA was extracted from the aboveground organs; (3) the plants were wild type; and (4) the material for RNA-extraction was collected at different stages of leaf development. Of 33 datasets in the GEO dataset collection that broadly satisfy these criteria (Table S3), two datasets were selected (Table S4). Clean paired- and single-end reads of Z. mays and O. sativa datasets were mapped to the corresponding reference genomes with the alignment efficiency of 70% and 89%, respectively (Table S5). This quality is sufficient for analysis of peroxisomal gene transcription during successive developmental stages.

First, to generate an overview of transcription patterns, we clustered the samples using the top 2000 genes that show the greatest expression variance. PCA analysis and scree plot of the samples are shown in Figure S2A-C; Figure S3–4. Z. mays samples clustered according to the leaf number and the leaf developmental stage (Figure 1a). We classified the 8th, 11th, and 13th leaves as mature leaves and all previous leaves as juvenile. The transcriptome of the 11th and 13th leaves sampled on the developmental stage of V9 were similar. Interestingly, the transcriptome of the 8th leaf sampled at stage V9 grouped with the 13th leaf of adult stages VT and R2. Further, the transcriptome of the immature leaf samples was similar to the 11th and 13th leaves sampled on the developmental stage of V9.

[IMAGE OMITTED. SEE PDF]

Differential expression analyses were carried out using 16 comparisons each based on generalized linear model implemented in DESeq2 package (version 3.6; see File S1 for the DEG locus IDs). The number of DEGs was higher at later developmental stages (Table S6). Next, we compared the GO term enrichment data of each cluster with the immature leaf cluster (Figure 1b; see File S2 for all comparisons). The immature leaf samples were used as a reference point for analysis of differentially expressed genes. The most common DEGs groups were RNA regulation, cell division, and metabolic regulation (Figure 1c). Clustering analysis of O. sativa samples showed high similarity on the fourth- and fifth-day transcriptome profile of the seedlings (Figure 2a,b). Ten GO terms were enriched between the fifth- and seventh-day transcriptomes (Figure 2c). GO terms related to biosynthetic process and response to the signals were more common between the seventh- and ninth-day (Figure 2d). GO terms related to cell cycle were more common between the sixth- and ninth-day (Figure 2e). The GO term enrichment analysis demonstrates that gene expression patterns in leaves at the juvenile stages share more similarities than at the adult stages.

[IMAGE OMITTED. SEE PDF]

Analysis of peroxisome metabolism and fission genes

Peroxisomal genes were categorized in 21 different processes (Table S8). Multifunctional genes were assigned a new category. For example, CAT3 that contributes to both photorespiration and ROS/RNS metabolism was assigned to a new category “photorespiration and ROS/RNS metabolism”. The numbers of differentially expressed genes in each process and how many DEGs are upregulated in each process are shown in Table S7 and Table 1.

TABLE 1 Ratio of differentially expressed genes (DEGs) and upregulated genes in the peroxisomal metabolism categories. Each row refers to each peroxisomal metabolism category. The comparisons were selected according to five developmental comparison column groups. Each column presents the comparison groups. Abbreviations are used for the column names. The numbers are the count of DEGs and that of upregulated genes, respectively. If only single gene was upregulated, the gene symbol was given in the Table 1 amino acid metabolism, 2-biosynthesis of phylloquinone, biotin, CoA and ubiquinone, 3-catabolism of polyamines, urate, pseudouridine, sulfite and methylglyoxal, 4-docking complex, 5-fatty acid breakdown, 6-fatty acid breakdown and phytohormone biosynthesis, 7-glyoxylate cycle and acetate–malate shunt, 8-NADPH regeneration, 9-NADPH regeneration and catabolism of polyamines, urate, pseudouridine, sulfite and methylglyoxyal, 10-peroxisomal import protein, 11-peroxisomal solute transporters, 12-peroxisomal solute transporters and biosynthesis of phylloquinone, biotin, CoA and ubiquinone and phytohormone biosynthesis, 13-peroxisome fission gene, 14-photorespiration, 15-photorespiration and reactive oxygen species (ROS)/RNS metabolism and autophagy, 16-phytohormone biosynthesis, 17-quality control and proteome remodeling, 18-recycling of PEX5, 19-recycling of PEX5 and pexophagy and phytohormone biosynthesis, 20-ROS/RNS metabolism, 21-signaling Ca²⁺, and protein phosphorylation.

	Total	11_im_v9	8_im_v9	13_im_v9	vT_v9_13	v9_R2_13	vT_R2_13	11_8_v9	13_11_v9	13_8_v9
1	14	64%| 43%	64%|43%	58%| 36%	71%|36%	57%| 21%	14%| 14%	43%| 31%	21%| 7%	57%| 36%
2	22	50%| 36%	46%|36%	50%| 36%	41%| 22%	50%| 32%	14%| 14%	36%| 23%	14%–0%	50%| 18%
3	13	46%| 39%	62%| 54%	46%| 31%	54%| 23%	46%|23%	8%| 8%	31%| 15%	0| 0	39%|15%
4	2	PEX14	PEX14	0| 0	0|0	PEX14	0| 0	0|0	0| 0	50%| 0%
5	19	79%| 68%	63%| 53%	47%| 42%	53%| 32%	58%|42%	5%| 5%	21%| 21%	0| 0	32%|21%
6	5	60% - 60%	80%| 40%	40%|40%	60%| %20	60%| 20%	0| 0	60%| 40%	0|0	60%| 40%
7	5	CSY3	80%| 40%	60%| 60%	80%|40%	80%| 40%	0| 0	40%|40%	20%| 20%	60%| 40%
8	13	39%| 31%	62%| 39%	39%| 23%	62%|39%	77%| 39%	0| 0	31%| 15%	20%| 0%	46%| 31%
9	2	BADH	BADH	BADH	BADH	BADH	50%| 0	BADH	0|0	BADH
10	6	33%| 33%	50%| 50%	0| 0	PEX19a	50%|50%	0|0	0| 0	0|0	17%| 0%
11	9	33%| 22%	78%| 56%	33%| 33%	44%|33%	22%| 22%	0| 0	11%| 0%	22%| 11%	56%|11%
12	1	0| 0	PXA1	0|0	PXA1	100%| 0%	0| 0	PXA1	0|0	PXA1
13	8	50%| 50%	62.5|50%	50%| 50%	62.5%| 37.5%	50%|37.5%	0| 0	DRP5B	12.5%| 0%	62.5%| 25%
14	10	80%| 70%	70%| 60%	70%|60%	90%| 60%	90%|60%	20%| 20%	70%|20%	40%| 0%	70%| 20%
15	3	67%| 67%	100%|67%	33%| 33%	100%|33%	100%| 33%	33%|33%	67%| 33%	33%| 0%	67%|33%
16	10	30%| 30%	40%| 40%	40%|30%	30%| 20%	40%| 30%	0| 0	0|0	0| 0	10%| 0%
17	3	67%| 67%	67%|67%	67%| 67%	67%|67%	LON2	0| 0	0|0	0| 0	33%| 0%
18	6	83%| 67%	83%| 83%	50%| 50%	50%|50%	33%| 33%	0| 0	17%| 0%	33%|0%	50%| 0%
19	2	0| 0	0|0	0| 0	PEX6	0|0	0| 0	0|0	0| 0	0|0
20	11	55%| 46%	82%| 73%	64%| 55%	55%|55%	46%| 46%	0| 0	36%–0%	9%| 0%	55%|0%
21	3	67%| 67%	67%|67%	67%| 67%	67%| 33%	67%|33%	0| 0	33%| 0%	0|0	33%|0%

Analysis of peroxisomal genes demonstrated the greatest difference between juvenile and mature leaves in both Z. mays and O. sativa (Figure 3a; Figure S6A–P). The early leaf developmental stages S1 (4 days after germination), S2 (5 days after germination), and S3 (6 days after germination) in O. sativa were more similar than later stages S4 (7 days after germination) and S5 (9 days after germination) (Figure 3b). The most represented pathways in 7–9 days old seedlings were biosynthesis of phylloquinone, biotin, CoA, and ubiquinone and photorespiration (S4 and S5, respectively) relative to 4–6 days old ones (S1, S2, and S3, respectively; Table 2). Similar correlation of peroxisomal gene expression with the developmental stage was observed in Z. mays (Figure 3c,d). Photorespiration was the most represented pathway in adult leaves relative to the juvenile stages (Table 2). Components of ROS/RNS metabolisms, NADPH regeneration, and catabolism of polyamines were also enriched at later stages of leaf differentiation. However, genes involved in glyoxylate cycle and acetate–malate shunt, and peroxisomal solute transporters were not enriched in the adult phase. During transition and reproductive stages, the most upregulated pathways were fatty acid breakdown and phytohormone biosynthesis. Interestingly, pathways responsible for catabolism of polyamine, urate, pseudouridine, sulfite, and methylglyoxal were enriched at all stages. This comprehensive analysis highlights developmental control of peroxisomal metabolic activities during transition from juvenile to adult stages.

[IMAGE OMITTED. SEE PDF]

TABLE 2 Constantly upregulated genes in the categories. All comparisons were assigned into three main developmental stages.

Category	Adult versus immature leaves	Leaves comparison in adult and reproductive phase	Early seedling comparison
Amino acid metabolism pathway	GLUTAMATE:GLYOXYLATE AMINOTRANSFERASE (GGT1)		GLUTAMATE:GLYOXYLATE AMINOTRANSFERASE (GGT1)
POWDERY MILDEW RESISTANCE GENE (PM16)	POWDERY MILDEW RESISTANCE GENE (PM16)
	SRY-TYPE HMG BOX (SOX)	SRY-TYPE HMG BOX (SOX)
	COBALAMIN-INDEPENDENT METHIONINE SYNTHASE (ATMS1)
	ANGIOTENSINOGEN 2 (AGT2)
The biosynthesis of Phylloquinone, biotin, Coa and ubiquinone	ACYL-ACTIVATING ENZYME (AAE) -AAE1, AAE12, AAE14, AAE17-	ACYL-ACTIVATING ENZYME (AAE)- AAE7, AAE12	ACYL-ACTIVATING ENZYME (AAE) AAE1, AAE7, AAE12
4-COUMARATE-COA LIGASE-LIKE (4CL1)	ATP CITRATE LYASE 1 (ACL1)
	CYTOSOLIC L-ASPARAGINASE I (ASP3)
	NICOTINAMIDASE 1 (PNC1)
	DDB1-CRBN E3 UBIQUITIN LIGASE (4CI3)	DDB1-CRBN E3 UBIQUITIN LIGASE (4CI3)
Catabolism of polyamines, urate, Pseudouridine, sulfite and Methylglyoxyal	SULFITE OXIDASE (SO)	POLYAMINE OXIDASES 4 (PAO4)	POLYAMINE OXIDASES 4 (PAO4)
6-PHOSPHOFRUCTOKINASE ISOZYME 2-TYPE CARBOHYDRATE KINASE (PFK2)		COPPER AMINE OXIDASe (CuAO)
GLYOXALASE 1 (GLX1)
Component of the protein docking complex	PEROXIN14 (PEX14)
Fatty acid breakdown and Phytohormone biosynthesis	ACYL-COA OXIDASE 5 (ACX5)	ACYL-COA OXIDASE (ACX) ACX3, ACX4, ACX5, ACX6	ACYL-COA OXIDASE (ACX) ACX4, ACX6
CITRATE SYNTHASE (CSY3)	CITRATE SYNTHASE (CSY3)
3-KETOACYL-COA THIOLASE 1 (KAT1)	3-KETOACYL-COA THIOLASE 1 (KAT1)
	ENOYL-COA DELTA ISOMERASE (ECI)
	ENOYL-COA HYDRATASE/ISOMERASE A (ECHIA)
	ENOYL-COA HYDRATASE (ECH2)
	DIENOYL-COENZYME A ISOMERASE (atDCI)
	SMALL CDPK-INTERACTING PROTEIN 2 (SCP2)
	ACETYL-COA ACETYLTRANSFERASE (ACAT1.3)	ACETYL-COA ACETYLTRANSFERASE (ACAT1.3)
	D-SPECIFIC MULTIFUNCTIONAL PROTEIN 2 (MFP2)
	GLYOXYSOMAL PROTEIN KINASE 1 (GPK1)
Glyoxylate cycle and acetate–malate shunt		MALATE SYNTHASE (MLS)	MALATE SYNTHASE (MLS)
	ISOCITRATE LYASE (ICL)	ISOCITRATE LYASE (ICL)
6-PHOSPHOGLUCONATE (6PGL)	6-PHOSPHOGLUCONATE (6PGL)
Peroxisomal protein import	PEROXIN19a (PEROXIN19a)	PEROXIN (PEX26)
NADPH regeneration and catabolism of polyamines, urate, Pseudouridine, sulfite, and Methylglyoxyal	BETAINE ALDEHYDE DEHYDROGENASE (BADH)	BETAINE ALDEHYDE DEHYDROGENASE (BADH)
	NAD(P)H DEHYDROGENASE B1 (NDB1)
	ZINC-BINDING DEHYDROGENASE (ZnDH)
	6-PHOSPHOGLUCONATE DEHYDROGENASE (6PGDH)	6-PHOSPHOGLUCONATE DEHYDROGENASE (6PGDH)
	SERINE-ASPARTATE REPEAT-CONTAINING PROTEIN C (sdrC)
		NAD(H) KINASE 3 (NADK3)
Peroxisomal solute transporters		PERIPHERAL MYELIN PROTEIN 22 (PMP22)
		LONG-CHAIN ACYL-COA SYNTHETASE (LACS) 6, LACS7	LONG-CHAIN ACYL-COA SYNTHETASE (LACS) 6, LACS7
Peroxisome fission genes	PEROXIN11a (PEX11a)	PEROXIN11a (PEX11a)
PEROXIN11b (PEX11b)		PEROXIN11b (PEX11b)
PEROXIN11c (PEX11c)	PEROXIN11c (PEX11c)
		PEROXIN11d (PEX11d)
Quality control and proteome remodeling	ENDOPEPTIDASE 15 (DEG15)
LON PEPTIDASE 2 (LON2)
PEX5 recycling	PEROXIN1 (PEX1)
PEROXIN12 (PEX12)
PEROXIN22 (PEX22)
Photorespiration	MALATE DEHYDROGENASE 1 (MDH1)	MALATE DEHYDROGENASE 1 (MDH1)	MALATE DEHYDROGENASE 1 (MDH1)
SERINE:GLYOXYLATE AMINOTRANSFERASE (SGAT1)	SERINE:GLYOXYLATE AMINOTRANSFERASE (SGAT1)	SERINE:GLYOXYLATE AMINOTRANSFERASE (SGAT1)
HYDROXYACID OXIDASE (HOAX1)	HYDROXYACID OXIDASE (HOAX1)	HYDROXYACID OXIDASE (HOAX1)
HYDROXYPYRUVATE REDUCTASE 1 (HPR)		HYDROXYPYRUVATE REDUCTASE 1 (HPR)
Photorespiration and ROS/RNS metabolism and autophagy	CATALASE (CAT3)	CATALASE (CAT3)
Phytohormone biosynthesis	INDOLE-3-BUTYRIC ACID (IBR) IB1	INDOLE-3-BUTYRIC ACID (IBR) IB1	INDOLE-3-BUTYRIC ACID (IBR) IB3
SULFOTRANSFERASE 1 (ST1)	ABSENT IN MELANOMA 1 (AIM1)
SULFOTRANSFERASE (ST5)	SULFOTRANSFERASE (ST5)	SULFOTRANSFERASE (ST5)
ROS/RNS metabolism	CU–ZN SUPEROXIDE DISMUTASE 3 (CSD3)		CU–ZN SUPEROXIDE DISMUTASE 3 (CSD3)
ASCORBATE PEROXIDASE 3 (APX3)
MULTIDRUG RESISTANCE (MDAR4)	MULTIDRUG RESISTANCE (MDAR1)
PROTEIN GAMMA RESPONSE 1 (GR1)	PROTEIN GAMMA RESPONSE 1 (GR1)
DEHYDROASCORBATE REDUCTASE (DHAR)
MULTIDRUG RESISTANCE 1 (MDAR1)
Signaling Ca2+ and protein phosphorylation	OPDA REDUCTASE 3 (OPR3)	OPDA REDUCTASE 3 (OPR3)

GABA promotes peroxisome proliferation

We reasoned that developmental regulation of each peroxisomal process can be assessed arbitrarily by the fraction of genes that exhibit differential expression during successive leaf development stages. The processes were ranked according to the percentage of DEGs at each developmental stage and the fraction of upregulated genes (Table 1; see Table S8 for the rank). Photorespiration was on the top of the list. This outcome is consistent with higher photosynthetic activity in mature leaves. Seven metabolic processes in mature leaves harbored over 50% of upregulated genes in at least one comparison group. Among the most consistently upregulated genes in all comparison groups were BADH and genes encoding components of ROS metabolism (Figure S5). Transcriptional upregulation of genes encoding peroxisomal pathways in mature leaves is consistent with the role of peroxisomes in photosynthesis as maintaining a higher number of peroxisomes per cell contribute to greater photosynthetic activity.

Upregulation of BADH transcription during leaf differentiation was verified by qRT-PCR in T. monococcum. CAT3 was used as a positive control. Both BADH and CAT3 were upregulated in adult leaves relative to the juvenile stages (Figure 4a). Upregulation of CAT3 is consistent with the requirement for scavenging of ROS produced during photorespiration. Upregulation of BADH could contribute to increasing GABA synthesis. The role of GABA in peroxisome biology remains unknown. Perhaps GABA functions in regulating peroxisome abundance during leaf development. To test this hypothesis, we treated juvenile T. monococcum leaves with 50 mgL⁻¹ GABA solution for 0, 10, and 30 min. Peroxisome abundance increased after 10 min of GABA incubation (Figure 4b,d; see Tables S9 for raw values and ANOVA test). GABA-induced peroxisome proliferation was accompanied by higher transcription of both peroxisome fission gene PEX11C and BADH (Figure 4c). On the other hand, peroxisome abundance was indistinguishable from the control after 30 min of exposure to GABA. GABA could be an initial trigger of peroxisome proliferation, and sustaining peroxisome abundance may require a different signaling mechanism.

[IMAGE OMITTED. SEE PDF]

Peroxisomal BADH mutant in A. thaliana lacked a discernible morphological phenotype and showed reduced root growth, and necrotic and purpling leaves under salinity stress (Ludewig et al., 2008; Zarei et al., 2016). We found that upregulation of ROS scavenging genes and BADH is accompanied by upregulation of genes involved in photorespiration. Therefore, BADH could play a role in ROS homeostasis by increasing synthesis of GABA. Further, both GABA and peroxisomes accumulate under stress conditions (Hickey et al., 2022; Su et al., 2019), suggesting that GABA contributes to ROS homeostasis by serving as a signaling molecule for activation of peroxisome proliferation. Concomitant increase of peroxisome abundance and accumulation of PEX11C transcript after GABA treatment demonstrate activation of peroxisome fission machinery. The next steps will focus on untangling the signaling pathways between GABA and peroxisome biogenesis.

AUTHOR CONTRIBUTIONS

Conceptualization, methodology, and formal analysis: Yunus Şahin, Andrei Smertenko, and Ercan Selçuk Ünlü. Investigation: Yunus Şahin and Taras Nazarov. Writing—original draft preparation: Yunus Şahin. Writing—review and editing: Yunus Şahin, Ercan Selçuk Ünlü, Andrei Smertenko, Taras Nazarov, and Nusret Zencrici. Supervision: Andrei Smertenko, Ercan Selçuk Ünlü, and Nusret Zencrici. Funding acquisition: Nusret Zencrici, Ercan Selçuk Ünlü, Yunus Şahin, and Andrei Smertenko. All authors have read and agreed to the published version of the manuscript.

ACKNOWLEDGMENTS

The authors are grateful to Kathleen Hickey for with construction of charts. This work was supported by 2214-A International Research Fellowship Programme for PhD Students granted by The Scientific and Technological Research Council of Turkey (TÜBİTAK) to YS; Scientific Research Projects Council of Bolu Abant Izzet Baysal University under Grant 2020.03.01.1435 to NZ; and USDA-NIFA hatch project #1015621, The Orville A. Vogel Wheat Research Fund, and Washington Grain Commission to AS.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.

PEER REVIEW

The peer review history for this article is available in the Supporting Information for this article.

DATA AVAILABILITY STATEMENT

The scripts were used in the study can be reached by the link; .

References

Ahn, S., Anderson, J. A., Sorrells, M. E., & Tanksley, S. D. (1993). Homoeologous relationships of rice, wheat and maize chromosomes. Molecular and General Genetics MGG, 2415(241), 483–490. [DOI: https://dx.doi.org/10.1007/BF00279889]