Klf9 is a key feedforward regulator of the transcriptomic response to glucocorticoid receptor activity

The zebrafish has recently emerged as a model system for investigating the developmental roles of glucocorticoid signaling and the mechanisms underlying glucocorticoid-induced developmental programming. To assess the role of the Glucocorticoid Receptor (GR) in such programming, we used CRISPR-Cas9 to produce a new frameshift mutation, GR369-, which eliminates all potential in-frame initiation codons upstream of the DNA binding domain. Using RNA-seq to ask how this mutation affects the larval transcriptome under both normal conditions and with chronic cortisol treatment, we find that GR mediates most of the effects of the treatment, and paradoxically, that the transcriptome of cortisol-treated larvae is more like that of larvae lacking a GR than that of larvae with a GR, suggesting that the cortisol-treated larvae develop GR resistance. The one transcriptional regulator that was both underexpressed in GR369- larvae and consistently overexpressed in cortisol-treated larvae was klf9. We therefore used CRISPR-Cas9-mediated mutation of klf9 and RNA-seq to assess Klf9-dependent gene expression in both normal and cortisol-treated larvae. Our results indicate that Klf9 contributes significantly to the transcriptomic response to chronic cortisol exposure, mediating the upregulation of proinflammatory genes that we reported previously.

A principal component analysis (PCA) of the RNA-seq data indicated that 60% of the variance in gene expression among samples is captured in the first two PCs, which respectively correlate with the two treatments (chronic cortisol and absence of a GR, respectively accounting for 46% and 14% of the variance; Fig. 2B). The VBA-samples again underexpressed nr3c1, indicating that the VBA screen was effective in identifying homozygous mutants (Fig. S2B). The eight VBA+ (i.e. mixed wild type and heterozygous mutant) samples cluster according to whether they were treated with cortisol, with the two clusters (cortisol-treated and vehicle-treated controls) segregating toward opposite poles of PC1. This is not the case in the VBA-(i.e. homozygous GR 369-) fish, the four cortisol-treated replicates of which are widely dispersed along PC1. This indicates that the global effect of the cortisol treatment captured in PC1 is GR-dependent and suggests as well that gene expression is less constrained overall in larvae lacking a GR. PC2 correlates with the presence and absence of a GR (VBA+   )): VBA+ vs. VBA-(vehicle-treated controls); cortisol-treated vs. vehicle-treated VBA+; and cortisoltreated vs. vehicle treated VBA-. Upregulated genes in common respectively between the first two and the last two comparisons are listed on the right. (D) Box plots of Z-transformed expression levels of klf9 and npas4a obtained from the RNA-seq data. Sample numbers correspond to those in (B). Asterisks indicate significant differential expression (adjusted p values): ***< .0001; *< .05.
Scientific RepoRtS | (2020) 10:11415 | https://doi.org/10.1038/s41598-020-68040-z www.nature.com/scientificreports/ vs. VBA− respectively). Interestingly, along PC2 chronic cortisol treatment displaces the VBA+ transcriptome toward the VBA-pole, suggesting that the cortisol-treated fish adapt to the exposure by developing GR resistance. The regulatory roles of the GR were further assessed by analyzing differential gene expression (DGE) between pairs of treatments, using an adjusted p value (false-discovery) threshold of 0.05 as the criterion for differential expression. Comparison of VBA+ and VBA-larvae identified 405 genes affected by loss of the GR in 5-dpf larvae, 103 of which are underexpressed (Fig. 2C) and 302 of which are overexpressed in VBA-larvae compared to their VBA+ counterparts (Fig. S3A). Gene ontology (GO) analysis shows that the underexpressed genes (i.e. genes normally upregulated by the GR) are involved in sugar metabolism and response to heat, whereas the overexpressed genes (i.e. genes that are normally downregulated via the GR) are involved in basement membrane organization, epidermis development, cell adhesion, locomotion, and growth (Figs. S3 and S4; Table S1).
A DGE analysis comparing cortisol-treated VBA+ fish to their vehicle-treated VBA+ counterparts showed that in cortisol-treated larvae with a functional GR, 4,298 genes were differentially expressed (Fig. S3B), 2,177 of which were upregulated (Fig. 2C) and 2,121 of which were downregulated. GO enrichment analysis of the upregulated genes identified biological processes related to nervous system development and function as well as cell adhesion, locomotion, and growth, whereas the downregulated genes were largely involved in protein synthesis and metabolism (Figs. S5, S6; Table S2). Interestingly, the transcriptome of cortisol-treated VBA+ larvae overlapped more with that of vehicle-treated VBA-larvae than with that of vehicle treated VBA+ larvae (Fig. S3C, Table S3), and accordingly, many of the biological processes affected by the absence of GR function in VBA-larvae were similarly affected by the chronic cortisol-treatment in VBA+ larvae (Fig. S8). This again suggests that the latter larvae develop a GR resistant phenotype. This resistance is probably not associated with any effect on levels of the GR or MR in the cortisol-treated larvae, as neither nr3c1 nor nr3c2 displays significant differential expression in response to the cortisol treatment.
We reasoned that some of the transcriptomic effects of the cortisol treatment might stem from GR-induced upregulation of a GR target gene that functions as a feedforward transcriptional regulator of GR signaling. To the extent that basal expression of such a gene depends on the presence of a GR it might be expected to be underexpressed in VBA-larvae (which lack GR function) but upregulated in VBA+ larvae in response to chronic cortisol (i.e., opposite of the predominant trend noted above). Of the 2177 genes upregulated in cortisol-treated VBA+ larvae only four were basally underexpressed in VBA-larvae (Fig. 2C), two of which encode transcription factors: klf9 and per1a. Of these, only klf9 was also found to be significantly upregulated in our previous RNAseq analysis of the effects of chronic cortisol treatment in wild-type fish (Fig. S3C), being one of the most highly upregulated transcription factors 9 . A plot comparing klf9 expression in each of the conditions reveals that the GR contributes to both its normal developmental expression and its upregulation in response to the cortisol treatment ( Fig. 2D), which was confirmed by qRT-PCR in another experiment (Fig. S2C). However, the plot also suggests that cortisol affects klf9 in a GR-independent fashion, albeit more variably, as indicated by the range of expression levels in the cortisol-treated VBA-samples shown in Fig. 2D, which correlate with the spread of the cortisol-treated VBA-samples along PC1 shown in Fig. 2B.
In contrast to the situation in VBA+ larvae, only 8 genes were differentially expressed in cortisol-treated VBA-larvae compared to their vehicle-treated VBA-siblings, all of them upregulated (Fig. 2C, Fig. S3D). The genes included the immediate early genes (IEGs) npas4a, egr1, egr4, fosab, and ier2b (Fig. 2C, Fig. S9). The GR-independence of their cortisol-induced upregulation is clearly seen in a plot of the expression levels of the most highly upregulated gene of this set, npas4a (Fig. 2D), a neuronal IEG that along with the other IEGs was also found to be upregulated in our previous RNA-seq analysis of cortisol-treated larvae (Fig. S3D) 9 . This indicates that the GR mediates nearly all the transcriptomic effects of chronically elevated cortisol, except for a small subset that appears to relate to neuronal activity.
A frameshift deletion introduced into exon 1 of zebrafish klf9 eliminates the DNA binding domain and significantly reduces expression of the mature transcript. The fact that klf9 was the transcriptional regulatory gene most consistently found to be upregulated by chronic cortisol exposure in a GR-dependent way suggested that it may contribute to the transcriptomic effects of the exposure. To test this, we mutated klf9 using CRISPR-Cas9 with a gRNA that targets exon 1 (Fig. 3A,B). This resulted in a frameshift mutation upstream of the DNA binding domain (Fig. 3A,B), producing a transcript encoding a truncated protein predicted to lack function as a transcription factor. Klf9 loss-of-function mutations are viable in mice 22 and similarly, the klf9 mutant fish were viable and fertile when bred to homozygosity, although mutant embryos survive at a lower rate than wild type (data not shown).
To ask how the mutation affects klf9 expression we used qRT-PCR to compare klf9 transcript levels in wild type and klf9 homozygous mutant (hereafter referred to as klf9 -/-) larvae under both normal conditions and in response to chronic cortisol treatment. This provided further confirmation that the cortisol treatment leads to upregulation of klf9 and revealed that klf9 mRNA levels are significantly reduced in the klf9 -/larvae, probably due to nonsense-mediated decay triggered by the premature stop codon (Fig. 3C). In support of this possibility, there was no significant effect on klf9 pre-mRNA levels, measured by qRT-PCR of the intron (Fig. 3C). We conclude from these experiments that the frameshift mutation introduced into klf9 exon 1 abrogates Klf9 function.
RNA-seq shows that klf9 is required for the pro-inflammatory transcriptomic effects of chronic cortisol treatment. To identify Klf9 target genes and ask whether Klf9 contributes to the transcriptomic response to cortisol treatment we used RNA-seq to query gene expression in 5-dpf wild type and klf9 -/mutant larvae from sibling parents, developed both normally and in the presence of 1 µM cortisol. Samples were collected at the same time on day 5 post-fertilization as in the previous GR knockout experiment and processed similarly. However, PCA revealed the largest source of variance in this experiment was not due to genotype or .0005  (Table S4), suggestive of a physiological stress response (e.g. to the stress of capture). However, a further confound is that PC1 also correlates with preparation of the RNA in two batches on separate days, suggesting that it may also reflect technical variance in sample preparation (see Materials and Methods). Reassuringly, after normalizing for this variation the samples segregate along two principal components representing genotype and treatment (Fig. S10B). For subsequent DGE analysis we therefore treated the batch effect as a categorical covariate. DGE analysis using a false-discovery rate of 0.05 identified 239 genes affected by loss of klf9 function in vehicle-treated larvae, 100 of which were upregulated and 139 of which were downregulated. Gene ontology term enrichment analysis shows that the upregulated genes are largely involved in complement activation (e.g. c3b, c3c, c4, cfb ), glucose and carbohydrate metabolism (e.g. pgm1, tpi1b, and pfkmb), and nucleosome positioning (e.g., h1fx, h1f0, and smarca5; Fig. S11, Table S5). Genes downregulated in response to loss of klf9 function are largely involved in sterol metabolism (e.g. faxdc2, sqlea, tm7sf2, sc5d, cyp51, lss, and msmo1; Fig. S12, Table S5). Interestingly, several of the processes associated with carbohydrate metabolism that we identified as being positively regulated by the GR are negatively regulated by klf9 (Fig. S13).
To ask how loss of kf9 function affects the transcriptomic response to cortisol treatment, we compared the response in wild-type embryos to that in klf9 -/larvae (Fig. 4A). Strikingly, with a false discovery threshold of 0.05 the major difference between the two responses was that ~ 70% (408) of the genes upregulated in by cortisol treatment in wild type embryos were not upregulated by the treatment in klf9 -/embryos, which instead upregulated a mostly different set of genes, albeit less strongly (Fig. 4A,B). This indicates that Klf9 is a key regulator of the transcriptomic response to cortisol. Of the 408 genes upregulated by the cortisol-treatment in a klf9-dependent way (Fig. 4A,B), about a quarter (91) were also identified in our previous study 9 as being significantly upregulated by chronic cortisol exposure (Table S6). Examples of the latter include irg1l and marco ( Fig. 4C and Fig. S14), as well as irg1, irf1a, ifi35, mpeg1.1, mpeg1.2, mxc, socs1a, socs3b, stat1b, and stat4 (Table S6). Gene ontology term enrichment analysis of the 408 genes upregulated by chronic cortisol treatment in wild-type but not klf9 -/embryos revealed significant enrichment for genes involved in defense and immunity (Fig. 4D, Table S7), the same biological processes that we previously found to be the most strongly affected by the chronic cortisol treatment 9 , whereas these processes were not associated with either the 176 genes upregulated in both wild-type and klf9 -/embryos, or the 903 genes upregulated in klf9 -/embryos but not wild-type (Table S7). Notably the set of genes upregulated by cortisol in wild-type but not klf9 -/embryos included four interferon regulatory factors (irf1a, irf8, irf9, and irf10), two interleukins (il1b and il34), and four interleukin receptors (il4r.1, il6r, il13ra1, and il20ra). These results indicate that Klf9 contributes in a significant way to the proinflammatory gene expression induced by chronic cortisol exposure.
Genes found to be consistently upregulated by chronic cortisol treatment in multiple RNA-seq experiments depend on klf9 for that upregulation. As a final analysis we assessed the overlap between the transcriptomic effects of chronic cortisol treatment across all our RNA-seq experiments with wild-type and VBA+ (mixed wild-type and heterozygous GR 369-) larvae, including the experiment published previously 9 , in order to identify a set of high-confidence genes that consistently respond the cortisol treatment and then ask how loss of klf9 function affects that response. To eliminate any technical artifacts emanating from the use of different parameters in the different analyses the sequence reads from all the experiments were reanalyzed from scratch using a common pipeline (see Materials and Methods). A PCA of the variance across all experiments produced several interesting observations. First, it showed that the effect of chronic cortisol treatment across all experiments was subtle, found only in PCs 4, 5, and 6 accounting for 14.3% of the total variance. PC5 captures a cortisol-treatment effect common to all three experiments (accounting for 2.77% of the total variance) and when plotted against PC4 clearly shows segregation of the cortisol-treated and control (vehicle-treated) samples (Fig. 5A). Gratifyingly, GO analysis of a single gene list ranked by upregulation along PC5 showed the same biological response to chronic cortisol as that which we reported previously 9 , i.e. upregulation of processes related to defense, inflammation and immunity (Fig. 5B, Table S8). The fourth principal component, accounting for 8.95% of the total variance, shows a cortisol treatment effect in the two experiments reported here, but not in that which we reported previously 9 ( Fig. 5A and Fig. S15). This suggests that PC4 represents cortisol-treatment effects that are dependent on the circadian light-dark cycle under which embryos in this study developed, which was absent in our previous study in which embryos developed in the dark 23 (see below, Discussion). The sixth PC (2.57%) is the complement of PC4, suggesting cortisol-induced effects that are only apparent in the larvae that were developed in the dark (Fig. S15).
The spread of the cortisol treatment effects across three principal components likely reflects differing biological responses to the treatment among different experiments. That the effects are somewhat different under total dark compared with light-dark cycles is not surprising given what is known about the interplay of GC signaling and the circadian clock [23][24][25][26] . Indeed GO analysis of PC2, which segregates our previously published experiment 9 from the two reported here (accounting for 18.5% of the variance), clearly shows that effects on circadian and light-responsive gene expression (Table S9). Additional sources of variance are less clear. The first PC, which segregates the wild-type samples in the last (klf9 -/vs. wild-type) experiment from the wild-type samples in our previous experiment as well as the VBA + samples described here (accounting for 37.6% of the variance, Fig. S15), is heavily loaded with genes involved in synaptic signaling and neurogenesis (Table S10) suggestive of differences in neurodevelopment and/or responsiveness to stress in those samples. The third PC (14.2% of the variance) segregates replicates 1 and 2 of the Klf9 experiment from all other samples, possibly reflecting the batch effect in sample preparation noted above.
Unsurprisingly given the large amount of gene expression variance across experimental samples unrelated to the cortisol treatment (i.e. noise), only 12 genes were found to be consistently upregulated by the cortisol  (Fig. 6A). The consistently upregulated set included klf9, the only gene in that set that encodes a transcription factor. We reasoned that many more genes are affected by the treatment, albeit not consistently with statistical significance owing to the abovementioned noise (see Discussion). This was borne out by plotting the estimated fold-change of the 149 genes that were upregulated by cortisol-treatment in at least two of the three experiments (Fig. 6B, Table S11), which also revealed that the upregulation of those genes is klf9-dependent. Furthermore, this plot shows that the cortisol treatment effect is stronger in the wild-type larvae than in the VBA+ (1:2 wild-type:heterozygous GR +/-369-), suggesting haploinsufficiency of the GR for some of the effect. Query of this set for enrichment of GO biological process terms indicated chronic cortisol exposure upregulates genes associated with response to organic substance (including response to stress, defense response,  www.nature.com/scientificreports/ response to external biotic stimulus, and inflammatory response to wounding), similar to what we reported previously 9 , and gluconeogenesis (Fig. 6C, Table S12). Finally, we used HOMER motif enrichment analysis 27 to ask what transcription factor binding sites are enriched in flanking regions of the set of 149 genes upregulated by cortisol-treatment in at least two of the three experiments. The resulting set of motifs included sites for the various krüppel-like factors as well as the GR (Table S13). The most significantly enriched motif, the Klf14 binding motif RGKGGGCGKGGC, matches the Klf9 consensus motif in the JASPAR database 28,29 and would be expected to bind Klf9, which is in the same KLF subfamily as Klf14 30,31 . In contrast, HOMER analysis of the 408 genes identified as klf9-dependent in Fig. 4 recovered a somewhat different list of motifs that included binding sites for several immunoregulatory transcription factors (Table S13; note that there are 65 genes in common between the 408 identified in Fig. 4 and the 149 identified in Fig. 6), suggesting that the larger set includes more indirect targets of feedforward regulation downstream of Klf9. Consistent with this, the latter set of motifs includes sites for two immunoregulatory genes, irf1 and stat4 that are both consistently upregulated by chronic cortisol treatment (i.e. in the common list of 149 genes identified in Fig. 6) and dependent on klf9 for that upregulation (i.e. in the list of 408 genes identified in Fig. 4). Altogether these results underscore the conclusion that klf9 is a feedforward regulator of GR signaling that mediates the pro-inflammatory transcriptomic response to chronic cortisol exposure, likely involving both   www.nature.com/scientificreports/ direct engagement of Klf9 with its transcriptional regulatory targets and downstream effects that those targets in turn have on the genes that they regulate.

Discussion
Several previously published studies have characterized loss-of-function mutations of the zebrafish GR, including four frameshifting indels that introduce a premature stop codon in exon 2 [13][14][15] (Fig. 1C). Here we report a new frameshift deletion (GR 369-) within exon 3 (Fig. 1A,B), which unlike the previous mutations removes all possible in-frame initiation codons upstream of the DNA binding domain (Fig. 1C), eliminating the transcriptional function of the mutant gene (Fig. 1E). We used RNA-seq to determine how this loss of GR function affects the larval transcriptome, both in larvae developed under normal conditions and in larvae treated chronically with cortisol. To overcome the low survival of embryos from homozygous mutant females and avoid nonspecific maternal effects that might be associated with poor egg quality, we took advantage of a visual background adaptation (VBA) screen to identify homozygous GR 369progeny of a heterozygous cross, as those larvae lack a VBA response (VBA-larvae). Based on the observation that VBA− larvae comprised ~ ¼ of the population, as expected for a recessive Mendelian trait, larvae that successfully mount a VBA response (VBA+ larvae) are predicted to consist of a 1:2 mixture of homozygous wild-type and heterozygous GR 369mutants, and hence to contain at least one intact nr3c1 allele encoding a functional GR. A principal component analysis of the RNAseq data revealed that absence of a functional GR has a profound effect on gene expression in both normal and cortisol-treated larvae, such that transcriptomes from larvae lacking a GR are clearly distinguished from those that have one, and that cortisol treatment produces a coherent effect only in larvae with a functional GR ( Fig. 2A). A pair-wise comparison of gene expression between VBA+ and VBA− larvae (accounting for PC2) identified with statistical significance about four hundred genes that are regulated by the GR in 5-dpf larvae at the time they were collected (midmorning). GO term enrichment analysis indicated that genes upregulated by the GR at that time are involved in metabolism and stress response, as would be expected, while genes downregulated by the GR are involved in epidermis development, cell adhesion, growth, and basement membrane formation, suggesting that the GR may function as a switch to downregulate those morphogenetic processes in late development, or to temporally segregate them to a certain time of day given the circadian dynamic of glucocorticoid signaling. Interestingly, numerous biological processes and individual genes affected by loss of GR function were similarly affected by chronic cortisol treatment in larvae with a GR. This suggests that one effect of the chronic treatment is to promote development of GR resistance, and moreover, that it does so via the GR. Such resistance might be construed as an adaptive response to the chronic exposure. Comparing gene expression in VBA+ larvae treated with cortisol versus vehicle (accounting for PC1 in that experiment) identified over four thousand genes that are differentially expressed in response to chronic cortisol treatment. This latter number is substantially larger than the 555 differentially expressed genes identified in our previous analysis 9 of the effects of chronic cortisol treatment and yielded a somewhat different result when subjected to GO term enrichment analysis (Figs. S5, S6). One major difference between the analyses reported here and that reported previously 9 is that in the latter the embryos and larvae were cultured in the dark, whereas in the present study we cultured them from fertilization in a diurnal light-dark cycle. In zebrafish larvae the circadian clock is not synchronized until the fish are exposed to a light-dark cycle 32,33 , so our previous results may have had circadian asynchrony as a confounding variable. Indeed, the impact of this difference on the transcriptome is clearly seen in the PCA of the combined analysis of all three RNA-seq experiments (Fig. S15), accounting for nearly 20% of the variance. Despite this, the RNA-seq results reported here assessing effects of chronic cortisol exposure in wild-type and klf9 -/larvae (Fig. 4) were from embryos developed with light-dark cycles, and in the wild type larvae produced an effect on proinflammatory gene expression similar to that of our earlier study, demonstrating that that effect was not an artifact of circadian asynchrony. Another difference between both studies using wild-type larvae and that depicted in Fig. 2 was that the latter measured transcriptomic effects of chronic cortisol treatment in a 1:2 mixture wild-type and heterozygous mutant larvae (VBA +); thus the results in VBA + larvae would be expected to be less sensitive to any effects of the chronic exposure for which the GR is haploinsufficient, a possibility supported by the comparison shown in Fig. 6B. Further work is required to more fully assess the effects of GR gene dosage on the transcriptome under both normal conditions and in response to chronic cortisol exposure.
The meta-analysis comparing all our RNA-seq experiments examining the transcriptomic effects of chronic cortisol treatment in wild-type or VBA + larvae (Figs. 5 and 6 and Fig. S15) provides some important insights that are broadly relevant to RNA-seq data interpretation. One is that different experiments that examine the effects of a single variable under somewhat different conditions and in a limited number of biological replicates will often produce different lists of differentially expressed genes passing an arbitrary threshold of statistical significance (e.g. adjusted p < 0.05). The reason is clear enough: biological systems are highly responsive to genetic and environmental factors that vary between experiments and affect gene expression, sometimes stochastically and/or in ways that are difficult to control and measure, especially with a limited number of biological replicates. Nevertheless, GO analyses of the lists obtained from different experiments can detect consistent biological effects following reanalysis of the data from all three experiments and their overlap. Only 8 of the 12 genes in common between all experiments were annotated with gene names (listed). (B) Violin and box plots of estimated fold change of the 149 genes that are upregulated in at least 2 of the three experiments shown in (A); the differences between each experiment are statistically significant (Table S14). (C) Treemap generated by REVIGO 57 (https ://revig o.irb.hr/) of GO biological process terms found by GOrilla 34 to be associated with the 149 genes upregulated by chronic cortisol treatment in at least two out of the three experiments.
Scientific RepoRtS | (2020) 10:11415 | https://doi.org/10.1038/s41598-020-68040-z www.nature.com/scientificreports/ even if the gene lists differ in the individual genes that they include, particularly if methods such as GOrilla 34 are used to test for statistically significant enrichment of GO terms toward one end or the other of a single list of genes ranked by some measurable criterion (e.g. a principal component that accounts for a specific condition as in Fig. 5). This underscores the important but often unappreciated point that statistical significance does not equate to biological significance, and generally is not a good sole criterion to assess the effects of a given condition on the expression of a given gene using high throughput methods such as RNA-seq. By comparing across the three experiments, we were able to identify both a small but "very high confidence" set of 12 genes as well as a larger "high confidence" set of 149 genes consistently affected by the cortisol treatment that would not have been discernible without our integrated meta-analysis (Fig. 6). Moreover, use of unbiased approaches such as PCA to parse the variance in the data can help identify robust condition-specific effects and provide insight into the biology underlying those effects when combined with GO term enrichment analysis. For the experiments reported here this approach validated our earlier finding that chronic cortisol exposure leads to upregulation of pro-inflammatory gene expression and extended that result by showing that the upregulation depends on the GR target gene klf9. The set of "very high confidence" genes identified in the meta-analysis included klf9, one of only four genes showing GR-dependent expression in normal 5-dpf larvae that were also upregulated in VBA+ larvae in response to chronic cortisol, and only one of two that encode transcription factors, the other being per1a (Fig. 2B). Both klf9 and per1a are involved in circadian regulation, and both have been shown to be GR targets in other vertebrate models 16,[35][36][37] . Interestingly, klf9 was also found to be the most commonly upregulated transcription factor gene in a recent meta-analysis of glucocorticoid-induced gene expression in the mammalian brain 38 . We have found by ATAC-seq that the promoter region of klf9 is one of the most differentially open regions of chromatin in blood cells of adults derived from cortisol-treated embryos (Hartig et al., submitted). In mice klf9 was recently shown to mediate glucocorticoid-induced metabolic dysregulation in liver 39 . Among other things Klf9 functions as a transcriptional repressor 40,41 , and in mouse macrophages as an incoherent feedforward regulator of the GR target klf2 42 , which functions to control inflammation 43 . Further work is needed to determine how klf9 contributes to pro-inflammatory gene expression in response to chronic cortisol exposure, which could either be directly as a feedforward activator (possibly via effects on metabolism), indirectly as a feedforward repressor of an antiinflammatory regulator like Klf2, or both. GO analysis also showed that genes involved in sterol biosynthesis are downregulated by loss of klf9, and more experiments are required to determine if klf9 regulates the metabolism of cortisol and other steroid hormones. Our motif enrichment analysis of flanking sequences from the set of 149 genes upregulated by chronic cortisol in at least 2 of our 3 RNA-seq experiments indicated enrichment for KLF binding sites (Table S13); further work involving chromatin immunoprecipitation is needed to determine whether any of those sites are bound by Klf9 or Klf2. Additional studies are also required to determine whether loss of Klf9 alters the function of immune cells or the inflammatory response to injury or infection.
Perhaps unsurprisingly, our results ( Fig. 2 and Fig. S3) indicate that the GR is required for nearly all the transcriptomic effects of chronically elevated cortisol. However, eight genes upregulated by the treatment were found in the RNA-seq analysis to be upregulated in both VBA+ and VBA− larvae (Fig. 2B), indicating that the GR is not required for their upregulation. Interestingly most of these genes are known IEGs, and include the neuronal activity-dependent gene npas4a, the mammalian homologue of which is directly repressed by the GR 44 , as well as egr1 which has been shown to differentially regulate GR in rat hippocampus depending on level of maternal care during development 45 . One possible explanation is that the genes are upregulated by increased MR activity, which was recently shown to contribute to stress axis regulation in zebrafish larvae 14 . Further work in MR mutant fish 14 will be needed to test this.
Finally, gene ontology analysis of genes upregulated by chronic cortisol treatment in VBA+ progeny of the GR +/369cross indicated a strong effect on biological processes associated with nervous system development and function. This is consistent with a recent report that injection of cortisol into eggs leads to increased neurogenesis in the larval brain 10 . In this regard it is interesting that klf9 is a stress-responsive gene that regulates neural differentiation and plasticity 35,46 . The long-term dysregulation of the HPA/I axis caused by early life exposure to chronic stress and/or chronically elevated cortisol suggests that the exposure perturbs brain development and activity. Given its role in regulating plasticity in brain regions relevant to neuroendocrine function, it will be interesting to determine whether klf9 contributes to those effects.

Materials and methods
Zebrafish strains, husbandry, and embryo treatments. The AB wild-type strain was used for all genetic modifications. Husbandry and procedures were as described previously 9 . All animal procedures were approved by the Institutional Animal Care and Use Committee (IACUC) of the MDI Biological Laboratory, and all methods were performed in accordance with the relevant guidelines and regulations. Embryo culture and cortisol treatments were performed as previously described 9 , with one difference: embryos were cultured in a diurnal light-dark cycle (14 h light-10 h dark). Briefly, fertilized eggs were collected in the morning, disinfected and at ~ 4 h post fertilization placed in dishes with either 1uM cortisol or vehicle (DMSO) added to embryo media. Embryos developed in a 28.5° C incubator with a 14/10 light/dark cycle synchronized with the core fish room. Media was changed daily.
Construction of nr3c1 and klf9 mutant lines. To mutate nr3c1 and klf9 we used CRISPR-Cas9 47 , injecting zygotes with multiple guide RNAs for each gene and mRNA encoding Cas9. Guide RNAs were designed using the CHOP-CHOP algorithm 48,49 .
To generate the GR 369mutant line fertilized wild-type AB embryos were injected at the 1 cell stage with 1-2 nL of a gRNA cocktail targeting nr3c1 exons 2 and 3 (final concentration 40 ng/µL for each gRNA, 230 ng/µL Scientific RepoRtS | (2020) 10:11415 | https://doi.org/10.1038/s41598-020-68040-z www.nature.com/scientificreports/ Cas9 mRNA, 0.1 M KCl, and 0.01% phenol red indicator dye). Individual whole injected larvae were screened for mutations of the targeted regions by high resolution melt analysis (HRMA) of PCR amplicons containing those regions 50 . Detected mutations were then verified by Sanger sequencing of a PCR amplicon containing the targeted region. F0 adults bearing mutations were identified by HRMA of DNA extracted from tailfin clips, and germline mutations were then identified by PCR and HRMA/sequencing of sperm. F0 males with germline mutations were outcrossed to AB females, and heterozygous progeny were screened via sequencing from tailfin clips. The GR 369mutation was identified and selected for by breeding over 2 additional generations to yield F3 homozygous progeny. Subsequent generations were maintained as heterozygotes for health and breeding purposes.
To generate the klf9 -/mutant line fertilized wild-type AB embryos were injected at the 1-cell stage with < 20% cell volume of injection mix consisting of 200 ng/ul Cas9 mRNA, 100 ng/ul guide RNA, 0.05% phenol red dye, and 0.2 M KCl. Individual F0 larvae were screened with HRMA to confirm CRISPR efficacy. Larvae were placed into system and raised. Young adult fish were genotyped via HRMA using DNA extracted from fin clips, and mutations were confirmed by Sanger sequencing. F0 fish positive for mutation were outcrossed to wild-type (AB) fish. F1 offspring of this cross were screened by HRMA to confirm germline transmission, and Sanger sequencing identified the 2 bp frame-shift mutation in one female founder. This female was out crossed with WT AB males. The resulting F2 fish were screened as young adults via fin clip and HRMA, and males positive for mutation were back crossed to the F1 founder female. Resulting F3 generation fish were screened via fin clip HRMA and sequenced to identify homozygous mutants as well as homozygous wild-type siblings. F3 generation was Mendelian 1:2:1 wild-type:heterozygote:homozygous-mutant ratio.
In vitro mRNA transcription and injection. Total RNA was extracted from 5dpf wild-type larvae using Trizol. RNA was treated with DNase I (NEB), and full-length nr3c1 cDNA was reverse transcribed using a specific primer (ggtcaaggttagtttaatgaattagtctgac) and Primescript RT kit (TaKara). Template DNA for transcripts was amplified from this cDNA using Q5 high fidelity polymerase (NEB), full-length (agtaatgcaaaatggatcaaggagg), truncated 310-(ctctttgggaacagctcgcc) or 369-(gggccagtttatgcttttcca) forward primers with upstream T7 promoters, and a common reverse primer (catcgtgtcctgctgttggg) downstream of the stop codon. Template DNA was run through an agarose gel to verify size, extracted and purified using E.Z.N.A. kit from Omega Biotek. Transcription reactions were run using mMessage mMachine Ultra T7 kit from Invitrogen, and yield quantified by spectrophotometry. Xenopus elongation factor 1a transcript from pTRI-Xef included in the mMessage mMachine kit was used as a control. Homozygous GR 369mutant embryos were injected at the one-cell stage with a mix containing 200 ng/ul mRNA and 0.05% phenol red dye, with an injection volume of ~ 20% cell volume. Injected embryos were snap-frozen at 6 h post fertilization for RNA extraction and qRT-PCR.
Visual background adaptation (VBA) screen. To identify larvae lacking a functional GR a VBA screen was performed on 4-day old larvae as described 12 . Briefly, larvae were incubated for 20 min in a dark incubator then transferred to a white background and immediately examined under a stereomicroscope with brightfield optics. Larvae that failed to mount a VBA response were identified by the failure of melanophores to disperse, remaining clustered in a dark patch on the dorsal surface 12 . Larvae were segregated as VBA+ and VBA-cohorts and returned to culture for an additional day before they were collected for RNA-seq.
RNA-seq and data analysis. At 3 h zeitgeber time (post lights-on) on day 5 post-fertilization four biological replicates of cortisol-treated and control embryos were collected as follows. For the first RNA-seq experiment, one by one, a dish of larvae corresponding to a single condition was removed from the incubator, and 8 larvae per replicate (4 replicates) were collected in a 1.5 mL tube with minimal water and immediately snap frozen in liquid nitrogen. All replicates from one condition were collected sequentially before moving on to the next condition. Total RNA was extracted using the Qiagen RNA-Easy Plus mini kit (Qiagen). For the second RNA-seq experiment, four replicates of n = 10 larvae were collected from a single dish for each condition and immediately snap frozen in liquid nitrogen. Collection of 16 samples occurred over 22 min. RNA was prepared as described above on two different days. On the first day (experimental replicates 1 and 2) the lysis buffer was added to all 8 frozen samples before homogenization, while on the second day the lysis buffer was added to each sample which was then immediately homogenized. This difference in sample preparation likely accounts for the large variance between samples prepared on day 1 and day 2. RNA was sent to the Oklahoma State Genomics Facility for Illumina library preparation and single-end sequencing.
RNA-seq libraries were generated with Illumina-compatible KAPA libraries and sequenced on an Illumina NextSeq 500 High Output sequencer. klf9 -/and matched control samples were sequenced as single end 75-bp reads. VBA+ and VBA-samples were sequences as paired-end 75-bp samples.
Fastq formatted read files were preprocessed with Trimmomatic 51 version 0.38 with default options, and then aligned to the Zebrafish genome version 11 as presented in Ensembl 52 version 93, using the STAR aligner 53 version 2.6.1b. The Ensembl transcriptome was preprocessed with a splice junction overhang of 100 nt. Following alignment, the resulting BAM files were processed with RSEM 54 version 1.3.0 for isoform and gene-level expression estimates. The resulting gene-level expression values were merged into a single expression matrix with an in-house python script. EDASeq 55 carried out in R version 3.6.1 was used to further normalize data for systematic effects, using gene-level length and GC-content as downloaded from Ensembl version 98 using EDAseq's included scripts. "WithinLane" normalization with GC content was judged as superior to that based on length. Final normalized gene-level counts (which = "full") were generated using GC-based WithinLaneNormalization followed by BetweenLaneNormalization. Subsequent differential expression analysis was carried out in R version 3.6.1 with the DESeq2 56 version 1.24.0, using either treatment (DMSO/cortisol) or genetics (VBA+/VBA− or klf9 -/-/WT) as the comparison. The klf9 -/experiment analysis also included a Scientific RepoRtS | (2020) 10:11415 | https://doi.org/10.1038/s41598-020-68040-z www.nature.com/scientificreports/ two-level categorical batch covariate. DESeq2 was also used to generate a rlog-matrix which was subsequently Z-transformed to normalize each gene across all samples. The Z-transformed rlog matrix was then loaded into JMP version 15 and used for PCA (using the "Wide Method") after thresholding to an average rlog expression value of 7.5. Both RNA-seq dataset have been deposited in the NCBI Gene Expression Omnibus database, under accession numbers GSE144884 (GR+/− experiment) and GSE144885 (Klf9+/− experiment).
The previously published RNA-seq data set was reprocessed as just described, and an overall matrix of only wild-type or VBA+ samples was generated and jointly normalized with EDASeq and then subjected to PCA as described in the previous paragraph.
The GOrilla algorithm 34 (https ://cbl-goril la.cs.techn ion.ac.il/) was used for Gene Ontology term enrichment analysis, and the data were visualized using REVIGO 57 (https ://revig o.irb.hr/). Venn diagrams were generated using Venny 2.1 (https ://bioin fogp.cnb.csic.es/tools /venny /). HOMER motif enrichment analysis 27 was used to compare incidence of known vertebrate motifs in a list of promoters of interest with incidence in a background list of all zebrafish promoters by running the findMotifs program (https ://homer .ucsd.edu/homer /micro array / index .html) using default settings except that sequence from − 1500 to + 500 bp relative to the transcription start site was searched for motifs from 10 to 18 bp in length.
Quantitative reverse transcription and polymerase chain reaction (qRT-PCR). Total RNA purified from snap-frozen larvae using the Trizol method and used as template to synthesize random-primed cDNA using the Primescript cDNA synthesis kit (TaKaRa). Relative gene expression levels were measured by qRT-PCR, using the delta-delta Ct method as described previously 9 , and eif5a as a reference gene. Examination of the results of multiple RNA-seq data sets indicated that eif5a activity was highly stable across treatments and genotypes. In many experiments beta-actin was also used as a reference gene, and this did not substantially change the results.
Graphics. RNA-seq results (PCA plots, box plots, scatter plots, and violin plot) were graphed using JMP version 15 from SAS. qRT-PCR results were graphed using Microsoft Excel. All figures were drafted using Adobe Illustrator CS4. Some of the graphs (Figs. 2B,D, 4C, and 5A) were redrawn in Illustrator exactly without modifying the depiction of the data.