Changing the strength of feedback in gibberellin production may be regulated by the circadian clock and carbon dioxide levels.

The Gibberellic Acid (GA) signaling pathway is a release of repression type hormone signaling pathway, where the receptor complex binds the ligand GA to negatively regulate DELLA repressors and activate GA-responsive transcription ([11,30], Fig 1A - black arrows). The DELLA repressors control GA homeostasis by feedback regulation of the GA pathway through at least three methods: 1) DELLAs activate transcription of the GA biosynthesis genes GA20 oxidase 2 (GA20ox2) and GA3 oxidase 1 (GA3ox1) [23], 2) DELLAs activate the transcription of the receptor GID1 [12,23], and 3) DELLAs repress their own transcription (Fig 1A - blue lines, [23,24]). We hypothesized that the strength of feedback regulation implemented by the DELLAs on GA perception and biosynthesis are important determinants of the steady state level of GA in a tissue. We further hypothesized that targeting a GAHACR to regulate the genes that mediate perception and biosynthesis (Fig 1A, blue arrows – DELLA to GA20ox, DELLA to GID1) would implement GA dependent repression and thus result in an effective decrease in the GA dependent transcriptional activation of these genes. For example, when targeted to a GA20ox promoter (Fig 1A, red dashed line), the GAHACR would repress transcription counteracting the activation through DELLA proteins in a GA-dependent manner (Fig 1B). We hypothesized that this feedback modulation would result in a decrease in the steady state levels of GA in planta as well as a change in the dynamics of GA accumulation.

To test our predictions, we simulated this perturbation by first implementing the model of GA signaling proposed by Middleton et al. in python. We then added additional differential equations to simulate GAHACR transcription, translation, and GA-dependent degradation paralleling the approach used for DELLA proteins. However, as the GAHACR is a repressor rather than an activator of the target gene, it's impact on transcription was appropriately altered (Fig 1C, Methods). We also built versions of this model without a GAHACR and without GA triggered degradation of the GAHACR, to simulate a dCAS9 repressor lacking the GA degron (No Degron CR) (Fig 1D). There is qualitative difference in both dynamic and steady state behavior of the GA biosynthesis pathway depending on if a GAHACR or a No Degron CR are targeted to a GA20ox locus. We modeled the GA4 concentration at steady state levels across a range of GAHACR expression and repression strengths (Fig 1E, 1F). The decrease in GA4 level when a GAHACR is targeted to GA20ox (compared to the no GAHACR control) is directly proportional to the expression level and repressive strength and inversely proportional to the degradation rate of the GAHACR (Fig 1E, 1F). There is a non-linear relationship between expression and repression strength's impact on GA4 levels when a GAHACR targets GA20ox as compared to the more linear relationship with a No Degron CR. This highlights the impact of having a feedback loop, and demonstrates that the GAHACR is predicted to be qualitatively different from a No Degron CR.

The use of the dCAS9 protein as the DNA-targeting domain of the HACR allows for rapid re-targeting of HACRs to genes in vivo [34]. This is useful when targeting genes in families with multiple redundant members such as the GA20ox gene family (five genes GA20ox 1–5, Fig 2A, [37]). We targeted the promoters of the four GA20ox genes that demonstrated the highest levels of transcription [11,21]) with sgRNAs assembled onto one T-DNA (GA20ox1–3 & GA20ox5, Fig 2B, S1A Fig). We introduced these T-DNAs into Arabidopsis GAHACR lines containing either a CAS9-TPLN (no degron) or a CAS9-GAdegron-TPLN (GAHACR) and observed significant reductions in primary root length (Fig 2C, 2D). We focused on the GAHACR retargeted to GA20 oxidases for simplicity, and because they allowed us to test the effect of tuning transcriptional feedback, which was our goal. We saw a significant reduction in root length for GAHACRs constructs using either GAI1 or RGA1 as the degrons (Fig 2E, 2F). We chose to further characterize two representative GAHACR lines with the RGA1 degron targeting GA20ox and observed consistent reductions in root length in both second and third generations after transformation, suggesting genetic stability of the components (sgRNA, GAHACR) of this repression scheme (Fig 2G).

We observe the predicted reduction in primary root length in lines that have a GAHACR targeted to GA20ox compared to parental controls without the gRNAs, consistent with lowered levels of GA (Fig 2D–2F). We tested the ability of exogenously applied GA to rescue the effect of the GAHACR repression of GA20ox and observed a rescue of root elongation (Fig 2H). This rescue by exogenous GA supports the effectiveness of the GAHACR as a conditional repressor. We sought additional support that GA20ox repression in adult plants resulted in GA-deficient phenotypes, and so we measured the length of time required for flowering as GA plays a positive role in this developmental transition. We observed a delay in flowering time in lines that have a GAHACR targeted to GA20ox compared to parental controls without the gRNAs in long day conditions (Fig 2I). As the contribution of GA to flowering time in long day conditions is relatively minimal compared to other inputs [7], we evaluated the effect of our perturbation in plants grown in short day conditions where the role of GA is much more pronounced. We observed a more significant delay in flowering times under short day conditions (Fig 2J), corroborating the previous observations under long day conditions.

The mathematical model predicted that repression of GA20ox genes would reduce the cellular concentration of bioactive GA (Fig 1). Phenotypic analysis of the lines with GAHACRs targeting GA20ox promoters demonstrated three lines of evidence that GA20ox activity had been reduced: (1) reduced root elongation (Fig 2C–2F) that is (2) rescued by exogenous GA application (Fig 2H), and (3) an extension in time before flowering (Fig 2I, 2J). In order to map the transcriptional landscape that is modulated by tuning down feedback between DELLAs and GA20 oxidases, we performed RNA-seq to quantify changes in the transcriptome. Seedlings were grown under long day conditions and RNA-seq was performed on 10-day old seedlings at ZT8 (Fig 3A). We observed a small number of differentially regulated genes (DEGs, 172 genes) that demonstrated clear overlap in behavior between the two biological replicates (Fig 3B, S2 Fig, S1 Table). When we examined this list of DEGs and intersected them against genes with altered expression in GA mutants, GA treatments, and Paclobutrazol treatments [38 –40]. We observed a partial overlap with known GA-related genes (82/172 genes, Fig 3C). When we examined the upregulated DEGs by Gene Ontology (GO) analysis we observed enriched GO terms related to photosynthesis, light perception or response, and circadian rhythms (Fig 3D, S3 Fig).

To better understand the network-level changes induced by the GAHACR on GA20ox activity we examined the relationships between DEGs to identify transcriptional regulators that could be responsible for the altered GO terms identified using Cytoscape [41] (S4 Fig). We identified a specific cluster of DEGs with related GO-terms 'Response to light' and 'Rhythmic process' that include well-studied components of the circadian clock PSEUDO-RESPONSE REGULATOR 5 (PRR5), PSEUDO-RESPONSE REGULATOR 7 (PRR7), GIGANTEA (GI), and REVEILLE 7 (RVE7). We also observed that two transcriptional regulators PHYTOCHROME INTERACTING FACTOR 4 (PIF4), and PHYTOCHROME RAPIDLY REGULATED1 (PAR1) are upregulated by the GAHACR intervention and could be causal to the changes in transcription related to photosynthesis and light perception. In order to verify these observations, we performed an independent experiment to test RNA levels of selected genes in 10-day old seedlings to examine whether we could replicate these changes in expression by qRT-PCR (Fig 3F). First, we examined the extent to which our GAHACR intervention reduced the expression of the GA20oxidases that they are directly targeting. We observed reductions in RNA abundance by ~50% for GA20ox1, GA20ox3 and GA20ox5, with a much weaker reduction in GA20ox2 abundance (Fig 3F). Next, we examined expression levels of GI, PIF4 and PRR7 and again saw small but significant increases in transcript abundance (Fig 3F). Because our data pointed to genes in the circadian clock having altered expression levels, we also tested the abundance of core clock genes CCA1 and TOC1 and observed that their abundance was slightly lowered (Fig 3F), suggesting that the amplitude or phase of clock components may be slightly altered by this GAHACR intervention.

To understand the connection between elevated atmospheric carbon dioxide and GA-responsive gene expression, we subjected our lines with GAHACR targeted to GA20oxidase to growth at elevated carbon dioxide levels. The original experiments demonstrated that elevated carbon dioxide levels rescued seedlings treated with the GA biosynthesis inhibitor paclobutrazol (PAC) in Arabidopsis [31,32], which inhibits ent-kaurene oxidase, an enzyme upstream of the GA20oxidase enzyme [11,42]. Therefore, we wanted to test whether elevated carbon dioxide levels would have the same effect on our genetic intervention on GA20oxidase activity. We grew seedlings on plates at elevated carbon dioxide levels (750 ppm) for 10 days and evaluated the effect of primary root growth (Fig 3G). We observed no rescue of primary root growth by elevated carbon dioxide levels at this time point, however we did note that primary root length was increased for both parental and experimental seedlings at elevated carbon dioxide levels (Fig 3G), consistent with previous studies that demonstrate increased biomass and root growth under elevated carbon dioxide levels due to higher photosynthetic rates [43]. However, we did observe a marked rescue of the delay in flowering time under long days at carbon dioxide levels (Fig 3H).

We performed an RNA-seq experiment on 10-day old seedlings grown under elevated carbon dioxide levels and observed differential gene expression specific to the GAHACR intervention (Fig 3I, S2 Table). These DEGs appear to be carbon dioxide level-specific as they are not fully shared with the ambient DEGs (S5, S6A, S6B Figs), suggesting that the molecular phenotype of the GAHACR intervention has been rescued at this stage even through differences in primary root growth have not (Fig 3G). Indeed, in an independent experiment we observed a rescue of GA20ox1, GA20ox2, and GA20ox3 mRNA levels, but not GA20ox5 (Fig 3J). Consistent with this observation we saw a five-fold increase in LUC gene expression in lines with GAHACR retargeted to GA20ox genes (Fig 3J). Furthermore, this rescue was observed for ambient DEGs GI, PIF4 and PRR7 as well as the core clock genes CCA1 and TOC1 (Fig 3J, S7 Fig). Additionally, an auxin-related up regulated DEG SAUR15 was also rescued by elevated carbon dioxide levels, suggesting the rescue of GAHACR intervention does not only rely on GA signaling (Fig 3J). When we compared our ambient and elevated carbon dioxide level RNA-seq data sets we observed a large number of differentially expressed genes (>5,000 DEGS, S6C Fig) with many that are circadian regulated using DIURNAL [44], which led us to question whether there was a misalignment in timing between these datasets, or whether carbon dioxide level could influence either the phase or amplitude of the circadian clock. To test this hypothesis we employed a transcriptional reporter of the circadian clock the pCCA1∷LUC under ambient and elevated carbon dioxide levels [45]. Surprisingly, under long day entrainment we observed a dampening of the amplitude of the pCCA1∷LUC reporter under elevated carbon dioxide levels (Fig 3K, light purple), and a slightly earlier max peak in phase (~30 minutes, difference in dashed lines) compared to ambient levels (Fig 3K, light blue). Taken together with the differences in RNA-seq data sets as well as independent RT-qPCR analyses, we suggest that atmospheric carbon dioxide levels may slightly influence the circadian clock, and that the GA feedback pathway appears to be downstream of this modulation.

To build our models we first implemented the differential equations described in the Middleton et al. in a python environment for ease of use. To simulate the GAHACR, we incorporated additional equations that were identical to those used to simulate the transcription, translation, and GA induced degradation of DELLA. A version of the model that set the GAHACR protein and mRNA concentration to zero was used to simulate wildtype regulation. Additionally, a version of the model where the GA dependent degradation terms associated with the GAHACR were set to zero was used to simulate the No Degron CR. To simulate the GAHACR's or No Degron CR's repression of GA20ox expression, we incorporated a term that captures GAHACR protein concentration scaled by a user defined repression strength constant into the denominator of the hill function used to simulate the activation of GA20ox by the DELLA protein [59]. The same constants used in the original Middleton et. al model was preserved in our model. Those associated with the GAHACR either mimicked those of the DELLA proteins or were programmatically varied across a range as reported in Fig 1C. Additional terms are described in the supplemental materials under Supplemental Model Information.

Expression cassettes for the gRNAs were cloned by golden gate methodology as described in [60]. The construction of the GAHACR was described previously [34]. The sgRNA expression cassettes contain sgRNAs driven by the U6 promoter and have a U6 terminator and were cloned as described in [60], with the modification that the dCAS9 cassette was eliminated from the final T-DNA by mutagenesis PCR. The guide RNA sequences in the GA20ox sgRNA plasmid are as follows: guide1: GA20ox5 – TATGATTTTGAAAACAAACG, guide2: GA20ox3 – AATTTATTTAGTGGCTGAAC, guide3: GA20ox1 – ATGGTCCTTTTAGTCTTTAT. The guide RNA sequences in the GID1 sgRNA plasmid are as follows: guide1: GID1A – GCGTTGAGGGATGAGTAGGG, guide2: GID1B – AAGAAAGACCAATCGGACGG, guide3: GID1C – GTGGATAAGAATATCGGCGT.

All GAHACR reporter lines were built by transforming the yeast artificial chromosome plasmids described above into Agrobacterium tumefaciens (GV3101) and using the resulting strains to transform a Columbia-0 background by floral dip [61] as described in [34]. Transformants carrying guide RNA constructs were then selected using a light pulse selection [62]. Briefly, this involves exposing the seeds to light for 6 hours after stratification (4°C for 2 days in the dark) followed by a three-day dark treatment. Resistant seedlings demonstrate hypocotyl elongation on Hygromycin. After selection seedlings were transplanted to soil and grown in long day conditions at 22°C. Two different auxin HACR backgrounds (independent transformant lines) were transformed with a sgRNA targeting GA20ox or GID1.

To characterize root growth phenotypes in plant lines with and without a GAHACR regulating target genes, we selected T2 transformants for lines that had a gRNA targeting GA pathway genes and the parental HACR background that had no gRNA. These seedlings were grown on selection plates for 4 days before being transplanted to media lacking selection. In all cases the parental controls that lack a gRNA and the lines derived from them were all grown in parallel and phenotyped on the same day to ensure the data collected was comparable. For flowering time experiments, T2 plants were grown on selection media for 7 days before being transplanted to soil and grown until flowering in a Conviron growth chamber with either long or short daylight settings. All plants that were phenotyped in soil were grown on Sunshine #4 mix soil in rose pots and watered every other day and checked every day for the presence of the inflorescence bolt at greater than 1 cm above the rosette, at which point they were scored as flowering. For RNA-seq experiments on seedlings grown at elevated carbon dioxide, T4 plants were grown on media without selection for 4 days to ensure proper germination before being transplanted to new plates lacking selection and grown in a Conviron growth chamber with either ambient or elevated carbon dioxide settings and long day light settings.

All qPCR assays were performed on seedlings grown in their stated conditions and then flash frozen in liquid nitrogen. Seedlings were homogenized in 2 mL tubes using a mixer mill bead beater with one steel ball bearing, at max speed for two 3-minute beating sessions, separated by a bath in liquid nitrogen to maintain freezing conditions. RNA was extracted from these seedlings using the Illustra RNAspin Mini Kit from GE (Via Millipore/Sigma Aldrich, Cat. No. GE25-0500-71). cDNA was then prepared from 1 μg of RNA using the iScript cDNA synthesis kit (Biorad, USA) and then used to run a qPCR with the iQ SYBR Green Supermix (Biorad, USA) on a qPCR machine (Biorad, USA). Each sample was analyzed for expression of target genes and AtPP2AA3 subunit which was used to normalize mRNA levels. A standard curve was generated using the pooled cDNA samples for each primer set to determine amplification efficiency and only primers with >90% efficiency were used. The primers used are listed in the supplementary primers table.

For ambient carbon dioxide concentrations, T4 seedlings were cultured on plates for 4 days without drug selection, and then transplanted to LS0 plates without drug for 10 days in either ambient or carbon dioxide supplemented Conviron (Pembina, North Dakota) growth chambers. RNA was extracted using the Illustra RNAspin Mini Kit from GE, and RNA-Seq was performed by Amaryllis Nucleics (Oakland CA). In brief, RNA was checked for quality and quantity on a Bioanalyzer, and poly(A) mRNA was purified from total RNA. The library was constructed using the Amaryllis proprietary construction kit, followed by QC by E-gel & Bioanalyzer. Libraries were pooled and run on an Illumina NextSeq 500 SR75, with a target of approximately 20 million reads per sample. Read preprocessing, mapping and differential gene expression was performed by Amaryllis, in addition to independent differential gene expression analysis in house using EdgeR and DESeq2. Network analysis was performed in cytoscape. GO terminology enrichment was calculated using g:profiler (biit.cs.ut.ee/gprofiler/gost) and plotted in R using ggplot2 (see code deposited at https://github.com/achillobator/GAHACR↗).

Luciferase based time course assays were used to characterize the dynamics of the CCA1 promoter driving luciferase (CCA1::LUC) under ambient and elevated carbon dioxide concentrations. CCA1∷LUC seedlings were grown in sterile 24 well plates in 0.5x LS0 media in either ambient or 750 ppm carbon dioxide levels for 10 days in long day conditions (6h-22h in light). These were then sprayed with luciferin (5 µM in water, Biosynth, United Kingdom, Cat. No. L-8220) in the evening of the final day (approximately 7 hours before midnight) and loaded into a Tecan (Spark model, Tecan Life Sciences, Switzerland), where they were sampled over time for 24 hours starting at time zero (midnight). The acquisition settings were as follows: Kinetic mode, Luminescence (Luciferin1), Interval time – 15 minutes, settle time – 0, Integration time – 1000 milliseconds, Output – counts/s. Data was collected for three separate runs of the experiment per genotype with 24 biological replicates per experiment. Data was analyzed in R, using the ggplot2 package (see code deposited at https://github.com/achillobator/GAHACR↗).