KLF6 and STAT3 co-occupy regulatory DNA and functionally synergize to promote axon growth in CNS neurons

The failure of axon regeneration in the CNS limits recovery from damage and disease. Members of the KLF family of transcription factors can exert both positive and negative effects on axon regeneration, but the underlying mechanisms are unclear. Here we show that forced expression of KLF6 promotes axon regeneration by corticospinal tract neurons in the injured spinal cord. RNA sequencing identified 454 genes whose expression changed upon forced KLF6 expression in vitro, including sub-networks that were highly enriched for functions relevant to axon extension including cytoskeleton remodeling, lipid synthesis, and bioenergetics. In addition, promoter analysis predicted a functional interaction between KLF6 and a second transcription factor, STAT3, and genome-wide footprinting using ATAC-Seq data confirmed frequent co-occupancy. Co-expression of the two factors yielded a synergistic elevation of neurite growth in vitro. These data clarify the transcriptional control of axon growth and point the way toward novel interventions to promote CNS regeneration.


Spinal injury and retrograde identification of corticospinal neurons. All animal procedures
were approved by the Marquette University Institutional Animal Care and Use Committee and complied with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. Cervical dorsal hemisections were performed as in 15,18 . Briefly, adult female C57/Bl6 mice (>8wks age, 20-22 g) were anesthetized by Ketamine/Xylazine and the cervical spinal column exposed by incision of the skin and blunt dissection of muscles. Mice were mounted in a custom spine stabilizer. Using a Vibraknife device 19 , in which a rapidly vibrating blade is controlled via a micromanipulator, a transection was made between the 4 th and 5 th cervical vertebrae, extending from the midline beyond the right lateral edge of the spinal cord, to a depth of 0.85 mm. To retrogradely identify CST neurons, immediately after injury a 2% solution of CTB-488 in sterile 0.9% NaCl (C22841-Thermofisher, Waltham, MA) was injected through a pulled glass micropipette fitted to a 10 μl Hamilton syringe driven by a Stoelting QSI pump (catalog # 53311) and guided by a micromanipulator (pumping rate:0.04ul/min). 0.5 μl of CTB-488 was injected at 2 sites located 0.5 mm rostral to the injury, 0.2 mm lateral to the midline, and to a depth of 0.8 mm. Unilateral pyramidotomy was performed as described in , Wang et al. 2015. Briefly, a ventral midline incision was made to expose the occipital bone, the ventrocaudal part of which was removed using fine rongeurs. The dura was punctured and the right pyramid cut completely using a micro feather scalpel. Viral delivery to cortical neurons. Cortical neurons were transduced using intracerebral microinjection as described in Wang et al., 2015). Mice were anesthetized with ketamine/xylazine (100/10 mg/kg, IP), mounted in a stereotactic frame, and skull exposed and scraped away with a scalpel blade. 0.5 μl of virus particles were delivered at two sites located 0 and 0.5 mm anterior, −1.3 mm lateral from Bregma and at a depth of 0.55 mm, at a rate of 0.05 ul/min using a pulled glass micropipette connected to a 10 µl Hamilton syringe driven by a programmable pump (Stoelting QSI). After each injection, the needle was left in place for 1 minute to minimize viral solution flow upwards.
Horizontal ladder. Animals were pre-trained in two sessions on a horizontal ladder (30 cm long, 1 cm rung spacing) which included four crossing events across the ladder, until their performance plateued as previously described 18 . Following gene treatment, the animals were tested 4, 8 and 12 weeks post treatment in four-crossing sessions with percent error quantified on the last three runs. All trials were video-recorded and errors were analyzed by blinded observers. Errors were defined as slips by the right forelimb (wrist breaks plane of ladder) and quantified as percentage of total steps by the right forelimb.
Immunohistochemistry. For KLF6 expression profiling across development, brain tissue was snap-frozen in liquid nitrogen and cryosectioned (25 µm of EGFP+ axons quantified in transverse sections of medullary pyramid, as in 15,18 . The number of EGFP+ profiles that intersected virtual lines at set distances from the injury site or midline of the spinal cord, normalized to total EGFP+ CST axons in the medullary pyramid, was quantified by a blinded observer on an Olympus IX81 microscope 15,20,21 . Exclusion criteria for pyramidotomy experiments was animals with less than 80% decrease in PKCγ in the affected CST, and for spinal injuries were 1) lesion depth less than 800 μm 2) axons with straight white matter trajectory distal to the lesion, and 3) EGFP+ axons in thoracic cord (too far for regenerative growth). Western Blotting. Frontal cortex was dissected from early postnatal and adult mice and sonicated ~5 s in 1 ml RIPA lysis buffer (N653-Amresco,Solon,OH) with serine and cysteine protease inhibitors (4693159001-Roche, Boston, MA). After 30 minutes of incubation on ice, cell debris was removed by centrifugation at 21,000 G at 4° and protein concentration estimated using the BCA method. 25 μg of protein was loaded and separated in a 4-20% SDS -PAGE gel (4561094-Bio-Rad, Hercules, CA) gel and transferred onto PVDF membrane (1620177-Bio-Rad, Hercules, CA) and blocked with 5% skim milk in 0.1% Tween 20 in PBS (TBST). The membranes were incubated with primary antibodies [mouse anti-KLF6 (SC-365633, Santa Cruz, Dallas, TX, 1:250)] in 5% skim milk in TBST overnight at 4 °C. The membranes were incubated with HRP-conjugated secondary antibody (115035146-Jackson ImmunoResearch, West Grove, PA) at a 1:1000 dilution for 2 hours at room temperature. The labeled proteins were detected using the ECL agent (34077-Pierce, Rockford, IL), following the supplier's manual. Briefly, a 1:1 solution of Supersignal west pico stable peroxide solution and enhancer solution was mixed and membranes were incubated in this solution for 7 minutes prior to detection. Detection and band quantification, was performed using an Odyssey Imaging system (LI-COR-Lincoln, NB) and band intensities were quantified using GAPDH (ab9484-Abcam, Boston, MA) as a reference with the built-in program. Quantification was the mean obtained from at least 3 biological replicates.
For the aggregated axon elongation assay, postnatal cortical neurons were first dissociated and transfected as described above. After transfection with test plasmids and an EGFP tracer, neurons were resuspended in Matrigel at a density of 133,000 cells/ul, and two microliter aggregates were placed on in a 35 mm petri dish, adjacent to a 2 ul line of matrigel that contained no cells. Explants were covered in 3 ml ENB media and maintained for 1 week, after which confocal microscopy produced an image that spanned the entire explant and extended 100 um above the petri dish surface (11 optical sections, spaced 10 um). The number of EGFP+ axons that crossed virtual lines at 1 mm intervals from the edge of the cell aggregate was quantified.
RNA-seq. P5 cortical neurons were prepared and maintained as described above. AAV8-EBFP or AAV8-KLF6 (titer-matched, 5*10 6 total units) were added to each well at the time of plating (300,000 cells/ well). To enrich for neurons, 100 nM 5-Fluoro-2′-deoxyuridine (FuDR) (F0503-Sigma-Aldrich, St.Louis,MO) was added 1 day post-transduction. NeuN immunohistochemistry confirmed >95% purity of nuclei (245 of 257 manually counted nuclei were NeuN+). At 3DIV, total RNA was extracted from neurons in culture using TRIzol reagent (10296028-Thermofisher, Waltham,MA) according to manufacturer's instructions followed by DNAase-I treatment (EN052-Thermofisher, Waltham,MA). Total RNA quantification (Q32852-RNA HS assay, Qubit Fluorometer-ThermoFisher, Waltham,MA) and quality assessment (50671512-RNA nano chips, Bioanalyzer, Agilent, Santa Clara, CA) were performed according to guidelines recommended by ENCODE consortia (goo. gl/euW5t4). All RNA samples used for library construction had RIN scores >=9. 100 ng of total RNA from three separate cultures were used to construct replicate libraries for each treatment. The polyadenylated fraction of RNA was enriched by bead-based depletion of rRNA using TruSeq Stranded Total RNA Library Prep Kit with Ribo-Zero according to manufacturer's instructions (RS-122-2201, Illumina technologies, San Diego,CA), and sequenced by Illumina HiSeq 2000 (University of Miami Genomics Core) to an average of 40 million paired-end reads. Preparation of cells for transduction and isolation of RNA were performed consecutively and by the same individuals to alleviate batch effects in sample preparation and sequencing.
Read quality post sequencing was confirmed using FASTQC package (per base sequence quality, per tile sequence quality, per sequence quality scores, K-mer content) 22 . Trimmed reads were mapped to the rat reference genome [UCSC, Rat genome assembly: Rnor_6.0] using HISAT2 aligner software (unspliced mode along with-qc filter argument to remove low quality reads prior to transcript assembly) 23 . An average of 80-85% reads were present in map quality interval of >=30, and between 8-10% reads were excluded due to poor map quality (<10). Transcript assembly was performed using Stringtie and assessed through visualization on UCSC Genome Browser. Expression level estimation was reported as fragments per kilobase of transcript sequence per million mapped fragments (FPKM).
Differential gene expression analysis was performed using Ballgown software 23 . Isoforms were considered significant if they met the statistical requirement of having a corrected p-value of <0.05, FDR <0.05. Gene ontology analysis on the top 100 significant differentially expressed genes (up and down-regulated) was performed using default parameters in EnrichR software 24 . Transcription factor binding site/motif analysis on KLF6 target gene promoters was performed using opossum v3.0 software 25 . Search parameters used were -JASPAR CORE profiles that scan all vertebrate profiles with a minimum specificity of 8 bits, conservation cut-off of 0.40, matrix score threshold of 90%, upstream/downstream of tss -1000/300 bps, results sorted by Z-score >=10. Network analysis on significant differentially expressed genes (upregulated) was performed on the cytoscape platform v3.5.1 26 , integrating multiple apps/plug-ins such as BiNGO, CyTransfinder, iRegulon, ClueGO, CluePedia and GeneMania 26-32 , using default parameters in a custom script.
TF footprinting. Published ENCODE consortia ATAC-Seq datasets generated using postnatal day 0 mouse cortices were used for TF footprinting analyses (ENCSR310MLB) 35 . Processed bam files were used as input for calling peaks using MACS2 software using the following parameters -macs2 callpeak -t <p0_atac>.bam -f BAM -g mm10 -n ATAC -q 0.01 36 . Position Weight Matrices for KLF6 and STAT3 downloaded from Cis-BP database 37 , narrow peak files generated using MACS2 and UCSC genome annotation mm10 were included as input to the CENTIPEDE footprinting algorithm, which was executed using default parameters and statistical thresholds 38 . CruzdB package was used to annotate genomic loci containing footprints into either promoter (2000 bp upstream of TSS/300 bp downstream of TSS) regions or enhancer regions (intersect chromosomal locations with mouse ENCODE cis-element database) 35,39 . Genes were separated into two categories -(1) Genes with only KLF6 footprints (2) Genes bound by both KLF6 and STAT3 in their promoter/enhancer regions, irrespective of binding distance. These genes were fed into the network analyses pipelines described above.

Experimental Design and Statistical Analysis.
For all cell culture experiments, neurite length from a minimum of 200 cells per treatment were averaged, and each experiment was repeated a minimum of three times on separate days. These averaged values were the basis for ANOVA with Sidak's multiple comparisons.
In vivo experiments were performed in a double-blind fashion, with non-involved lab personnel maintaining blinding keys. For in vivo experiments, animal numbers were determined by power analysis with G*Power software (α = 0.05, Power = 0.8) using historical estimates of variability and adjusted for prior rates of animal attrition 15,18 . Animals were randomized prior to viral treatment, with each surgical day including equal numbers from each group. The in vivo pyramidotomy experiment was initiated with ten female animals per group. One animal (control) was lost to mortality, and three animals (one control, two KLF6) animals were excluded on the basis of failed pyramidotomy (less than 80% reduction of PKCgamma signal in the affected tract). Spinal injury experiments were initiated with 18 control and 18 KLF6 animals. Three animals (two control, one KLF6) were lost to mortality, and three (two control, one KLF6) were excluded on the basis of insufficient injury depth and visibly spared axons. Axon counts (four sections per animal) at varying distances were compared between groups by RM ANOVA with Sidak's multiple comparison's test. For KLF6 expression across time, tissue was harvested from littermates at different time points. For the RNA-seq data validation, comparisons in fold change were analyzed with One-way ANOVA with Sidak's multiple comparison test. For all software used for bioinformatic analyses, default inbuilt statistical models were employed as described above. All statistical tests were performed with Graphpad Prism software v7.0(San Diego, USA) and all error bars are ±SEM and statistical significance was accepted with p-value <0.05.

Cortical KLF6 expression declines developmentally and is not increased by spinal injury.
Consistent with a positive role in axon growth, KLF6 transcript is highly expressed in the embryonic brain during periods of axon growth and then developmentally downregulated throughout the CNS during postnatal development, including in CST neurons 40,41 . A more recent study, however, reported detection of KLF6 protein by immunohistochemistry in the adult brain 42 . We therefore re-examined KLF6 expression at the protein level across early postnatal development using immunohistochemistry western blotting, and quantitative PCR. As expected, KLF6 was readily detectable by immunohistochemistry (red) in early postnatal forebrain neurons (Fig. 1A), and western blot analysis of whole cortex detected protein of ~40 kd, consistent with KLF6's expected size 43 (Fig. 1A,C). Consistent with a prior report 42 KLF6 protein co-localized with NeuN, indicating neuronal expression (Fig. 1A, insets). In the adult mouse brain, KLF6 was also detectable by immunohistochemistry (red), but with intensity much dimmer than postnatal tissue processed identically (Fig. 1B). KLF6 was also detected by western blot in the adult cortex, with expression about 10% of that in postnatal cortex (Figs 1C and S5) (p < 0.001, ANOVA with post-hoc Tukey's). Consistent with this, qPCR analysis of P1, P8, P15, P21, and adult cortex showed a 70% reduction in KLF6 transcript during early postnatal development (p < 0.01, ANOVA with post-hoc Tukey's, Fig. 1E). Combined, these data confirm a substantial downregulation of KLF6 protein and transcript in the developing cortex.
To determine whether KLF6 expression in CST neurons is affected by spinal injury, fifteen C57/BL6 wildtype animals received sham or cervical C4/5 hemisection and injections of the retrograde tracer CTB-488 to cervical spinal cord. Animals were sacrificed between three days to four weeks post-injury (3 animals per time point), and retrogradely identified CST neurons in the cortex were examined with immunohistochemistry for KLF6. Immunohistochemistry was performed simultaneously on all tissue sections, which were visualized with identical acquisition parameters. To account for residual variability, KLF6 intensity in the nuclei of identified CST neurons was normalized to the within-slice average intensity of non-injured cells in layer 4 (Fig. S1A). KLF6 signal remained dim in injured CST neurons at all time-points post-injury ( Fig. S1C-H). Thus, similar to KLF7 and other regeneration-associated genes 14,15,44 , KLF6 is developmentally downregulated and expression levels remain depressed in adult CST neurons after spinal injury. Forced KLF6 expression promotes CST sprouting and regeneration after spinal injury. We next tested the effects of overexpressing KLF6 on neurite outgrowth in postnatal rat cortical neurons, a well-characterized model system 45,46 . Rat neurons were used for cell culture experiments to improve per-animal cell yield and viability compared to mice; gene function in this rat system has repeatedly proven to predict pro-regenerative activity in mice 6,14,15 . First, to establish baseline expression of endogenous KLF6 expression, cortical neurons were cultured for two days and immunohistochemistry was performed with antibodies to KLF6 and neuron-specific β-III tubulin. High content-screening microscopy was used to quantify endogenous KLF6 signal (Fig. S2E) and neurite growth in thousands of individual β-III+ neurons (Fig. S2D). Neurons displayed a range of endogenous KLF6 signal intensity. Consistent with a pro-growth role for KLF6, neurons in the highest quartile of KLF6 intensity produced neurites that were significantly longer than the lowest quartile (top quartile 116.2% of bottom, p < 0.001, ANOVA with post-hoc Sidak's) (Fig. S2A,B).
To test more directly the effect of KLF6 overexpression, postnatal cortical neurons were transfected with plasmid DNA encoding KLF6 or EBFP (enhanced blue fluorescent protein) control, and automated microscopy was used to quantify neurite outgrowth (Fig. S2F). To distinguish transfected cells a fluorescent EGFP (enhanced green fluorescent protein) reporter was co-transfected at a ratio of 1 part EGFP to 4 parts test plasmid; this ratio was previously established to result in greater than 90% co-expression 45 . To determine how the response to KLF6 might vary during early postnatal development, neurons were collected at either postnatal day 3 (P3) or at P7, the oldest age at which neurons could be reliably harvested and transfected in our hands. Notably, in baseline control conditions P3 neurons extended neurites that were significantly longer than P7 (Fig. S2H). At P3, KLF6 overexpression produced a non-significant trend towards increased neurite length, but significantly increased neurite lengths in P7 neurons (Fig. S2H) (EBFP 201.7 ± 17.5 µm, KLF6 283.4 ± 28.2 µm p < 0.05, ANOVA with Sidak's). We attempted to simplify labeling of KLF6-expressing cells by creating a construct in which a T2A "self cleaving" oligopeptide linker was inserted between KLF6 and an mCherry red fluorescent protein 47 . This approach creates two distinct proteins from a single vector and has been used by our lab and others to study regeneration-associated genes 15,17,18 . We found, however, that that growth promotion by KLF6-2A-mCherry was substantially lower than untagged KLF6 (KLF6 137% ± 2.7SEM, p < 0.01 versus mCherry control; KLF6-2A-mCherry 111% ± 7.9SEM, p > 0.05 versus mCherry control, ANOVA with post-hoc Sidak's). The mechanism for this difference is unclear, but based on these data we elected to use untagged KLF6 constructs for all subsequent experiments, relying on EGFP from a separate plasmid to mark transfected neurons. Cell-based HCS (high-content screening) analysis of immunohistochemistry confirmed increased expression of KLF6 in transfected neurons (Fig. S2G). Collectively these data establish a correlation between endogenous KLF6 expression and total neurite length in cortical neurons, and confirm that elevated KLF6 expression is sufficient to enhance neurite length in postnatal cortical neurons.
We next tested the ability of forced KLF6 expression to enhance sprouting by adult CST neurons in vivo. Twenty adult female mice received unilateral pyramidotomy to deprive the right side of the spinal cord of CST input, followed by cortical injections of AAV8-EBFP control or AAV8-KLF6 to the right cortex in order to target uninjured CST neurons that project to the left spinal cord. AAV8-EGFP was co-injected as an axon tracer, using a 1:3 ratio shown previously to result in >90% co-expression 15,18 . Twelve weeks later, transverse sections of cervical spinal cord were prepared, and PKCgamma immunohistochemistry was performed to confirm successful ablation of the left CST ( Fig. 2A,B insets) 15,20,21 . One animal (control) was lost to mortality, and three animals (one control, two KLF6) were excluded for incomplete pyramidotomy. The total number of EGFP+ CST axons was similar between groups (EBFP: 2031.3 ± 158.6 SEM, KLF6: 2168.0 ± 161.3 SEM, p > 0.05, t-test). Compared to EBFP control, KLF6 expression caused significant elevation of cross-midline sprouting by CST axons at distances up to 600 µm from the midline (N = 8 per group, p < 0.05 RM ANOVA with Sidak's multiple comparison Fig. 2C). To test whether this enhanced sprouting was accompanied by improvements in forelimb function, animals were tested on a horizontal ladder task, which quantifies the percent of steps by the affected right forelimb that are mistargeted. As expected, pyramidotomy injury significantly increased the slip rate, but animals treated with KLF6 did not differ significantly from control animals at any timepoint post-injury (Fig. 2D). Thus forced KLF6 expression enhances compensatory sprouting by intact CST neurons, but does not yield improved forelimb targeting in the horizontal ladder task.
Next, to test the effect of forced KLF6 expression on regeneration by directly injured neurons, AAV8-EGFP tracer and AAV8-EBFP or AAV8-KLF6 were delivered to the left cortex of adult female mice (18 animals per group). One week later, CST axons in the right cervical spinal cord were transected by a dorsal hemisection injury (Fig. 3A). Animals were sacrificed twelve weeks later, and immunohistochemistry readily confirmed expression of KLF6 protein in sites of AAV injections, but not in AAV8-EBFP animals (Fig. 3B,C). Sagittal sections of cervical spinal cord were prepared, and GFAP immunohistochemistry (blue) was used to visualize the site of injury. Three animals (two control, one KLF6) were lost to mortality, and three (two control, one KLF6) were excluded on the basis of insufficient injury depth and visibly spared axons. To assess axon regeneration, the number of EGFP-labeled CST axons, normalized to the total number EGFP-labeled CST axons in the medullary pyramid, was quantified at set distances distal to the injury. The total number of EGFP+ CST axons detected in the medullary pyramids was similar between groups (EBFP: 2681.9 ± 155.4 SEM, KLF6: 3228.5 ± 240.4 SEM, p > 0.05, t-test). As expected, EBFP control animals displayed very little CST growth distal to the injury site (Fig. 3D) In contrast, KLF6-treated animals showed significantly elevated CST growth up to 3 mm beyond the injury site (N = 14 EBFP control, 16 KLF6, p < 0.05 RM ANOVA with Sidak' s multiple comparison, Fig. 3E,F). Similar to the pyramidotomy experiment, however, KLF6 treatment had no effect on forelimb placement in the horizontal ladder task at any time point post-injury (Fig. 3G). Overall, these data indicate that forced re-expression of KLF6 in CST neurons promotes axon sprouting and regeneration, but is not sufficient to improve forelimb placement.

RNA-seq analysis of transcriptional effects of KLF6 overexpression. Although KLF transcription
factors have been functionally linked to regenerative axon growth, the underlying molecular mechanisms remain unclear. We therefore adopted an RNA-seq approach to examine the transcriptional consequences of forced KLF6 Scientific RepoRts | (2018) 8:12565 | DOI:10.1038/s41598-018-31101-5 overexpression. Postnatal cortical neurons were treated with AAV8-mCherry control or AAV8-KLF6 virus, using titers that resulted in >95% efficiency of transduction. In addition, FUDR was added to remove contaminating glia, resulting in >95% neuronal purity as assessed by NeuN immunoreactivity. Three days later, RNA (RNA Quality Score RIN >9.0) was collected from three experimental replicates from three separate litters, and deep sequencing performed with Illumina HiSeq 2000, with an average of 40 million 100 bp reads per sample. After extensive quality control in accordance with ENCODE guidelines (RNA integrity, library quality, sequence read quality and distribution) and alignment to the rat genome (see methods and Fig. S3), differential gene expression was assessed using Ballgown software. KLF6 overexpression significantly altered the expression of 454 genes (p-value < 0.05, FDR < 0.05), 55% of which were upregulated and 45% down-regulated. The genome-wide fpkm values for all samples and the complete set of differentially expressed transcripts are available in (Sup. Table 1).
To validate the RNA-seq data, newly prepared cortical neurons were transfected with KLF6 or EBFP control, and qPCR was used to assess changes in genes predicted by RNA-seq to be most affected by KLF6 expression. As an additional control, qPCR for putative KLF6-downstream genes was also performed after transfection with mutant forms of KLF6. Transcriptional activity of KLF6 and -7 has been previously shown to be sensitive to acetylation of two lysines (K) located in the first zinc finger of the DNA binding domain 48,49 . Accordingly, lysines 209 and 213 of KLF6 were mutated to either arginine (R) or glutamine (Q) to block or mimic the effects of acetylation, respectively (Fig. S4A). In assays of neurite outgrowth, mutation to glutamine (KLF6-K-to-Q) completely blocked growth promotion by KLF6, whereas transfection with the arginine mutant (KLF6-K-to-R) elevated growth promotion modestly but significantly above the level of wildtype KLF6 (Fig. S4B). qPCR analysis showed that transfection with either wildtype KLF6 or KLF6-K-to-R significantly upregulated s100a10, ccnd1, Igals3, and s100a11, with fold changes similar to those in the RNA-seq dataset (Fig. S4C). In contrast, transfection with KLF6-K-to-Q showed no significant elevation in putative KLF6 target genes (Fig. S4C). qPCR for KLF6 confirmed successful transfection with all constructs (fold increases of 189.2 ± 43.05 SEM, 631.28 ± 101 SEM, 434.42 ± 276 SEM, for wildtype, K-to-Q, and K-to-R mutants respectively). Overall these data are consistent with the previously established importance of these zinc finger lysines for KLF6 activity 48 , and serve to independently validate the RNA-seq-based identification of differentially expressed genes. As a first step in clarifying the mechanisms by which KLF6 affects neurite outgrowth we performed gene network analysis of the genes most up-and down-regulated by KLF6. Network analysis identified 5 distinct subnetworks within the upregulated cohort. Interestingly, genes within each subnetwork were highly enriched for distinct functions associated with axon extension, including regulation of cytoskeleton remodeling, regulation of cellular component size, lipid synthesis and transport, membrane raft assembly, and bioenergetics (Fig. 4A). In contrast, we identified 3 distinct sub-networks within the downregulated gene cohort, all enriched for functions related to synaptic transmission (Fig. 4B). Overall these data indicate that KLF6 expression may affect axon growth primarily by altering the abundance of actin regulatory genes, while also engaging complementary pro-growth mechanisms. To test whether forced expression of cytoskeletal-interacting genes might phenocopy KLF6 growth promotion, in follow-up experiments the cytoskeletal-interacting genes most affected by KLF6 were overexpressed in cortical neurons. Unlike KLF6, however, expression of these putative effector genes did not significantly alter neurite lengths (Anxa1 107.5 ± 7.5%, Anxa2 100.5 ± 2.7%, S100a4 97.4 ± 10.0%, S100a10 95.8 ± 3.5%, N > 200 cells in three experiments, p > 0.05 ANOVA with post-hoc Sidak's). One likely explanation is that growth promotion by KLF6 may require changes in the abundance of multiple, functionally interacting genes. If so, single or even multi-expression of putative target genes may not be a feasible strategy to mimic the KLF6 effect by engaging downstream effectors.
Functional Synergy between KLF6 and STAT3. As an alternative to clarifying downstream cellular effects, analysis of differential gene expression can also shed light on KLF6 function by identifying additional transcription factors that may functionally interact. We therefore performed motif analysis of the promoter regions of genes upregulated by KLF6 expression, in order to identify factors that could potentially synergize by co-occupying promoters of KLF6-responsive genes. Consistent with regulation by KLF6, the motif most enriched in these promoter regions was the recognition sequence for KLF7/KLF6 itself. Also over-represented were recognition motifs for additional factors including KLF4, SP1, CTCF, and STAT3 (Fig. 5A). To test potential functional interactions, we overexpressed these factors singly or in combination with KLF6 in assays of neurite outgrowth. Co-expression of KLF6 with KLF4 has been previously tested 14 , and we therefore focused here on SP1, CTCF, and STAT3. EGFP reporter was co-expressed with the factors to identify transfected cells, using 1:4:4 plasmid ratios that resulted in approximately 80% co-expression of both test plasmids in EGFP+ neurons (Fig. 5B). As expected, single overexpression of KLF6 significantly increased neurite lengths (123.9% of control, p < 0.01, ANOVA with post-hoc Sidak's) (Fig. 5C). SP1 and CTCF had no effect on neurite outgrowth when expressed singly, and when combined with KLF6 did not differ significantly from the effect of KLF6 alone (p > 0.05, ANOVA with post-hoc Sidak's). To test STAT3 function we used two constitutively active forms, based on prior work indicating that active STAT3, but not wildtype, affects neurite growth in this culture system and others 50,51 . Consistent with prior reports, single overexpression of VP16-STAT3 significantly (albeit modestly) increased neurite lengths to 115.6 ± 2.7% of control, whereas caSTAT3 had no significant effect (106.9 ± 2.8%, p < 0.05 ANOVA with post-hoc Sidak's). Remarkably, co-expression of caSTAT3 and VP16-caSTAT3 with KLF6 increased neurite lengths to 161.2.6 ± 5.6% and 175.4 ± 9.6% of control, respectively (p < 0.01 compared to KLF6 alone, ANOVA with post-hoc Sidak's) (Fig. 5B,C). This degree of neurite outgrowth is well above that predicted by the additive effects of cSTAT3 and KLF6, indicating functional synergy. Regulatory network analysis of genes upregulated after KLF6 overexpression revealed sub-networks enriched for distinct functional categories highly relevant to axon growth. Nodes correspond to target genes and edges to multiple modes of interaction (physical, shared upstream regulators, shared signaling pathways and inter-regulation). Node color represents their corresponding functions as denoted in the legends and node size is based on significance for enrichment in functional category. Only significantly enriched GO categories were included in the network analysis (p < 0.05, Right-sided hypergeometric test with Bonferroni correction, Cytoscape).
As a further test of synergy between KLF6 and VP16-STAT3, we tested co-expression in a second growth assay. Postnatal cortical neurons dissociated and transfected as previously, but then re-aggregated in drops of matrigel to create artificial explants. These explants were placed adjacent to lanes of matrigel that contained no cells. The key advantage of this assay is that unlike individually plated cells, aggregated neurons display sustained, long distance axon growth that extends several millimeters. This provides a better model of long distance axon regeneration in vivo (Fig. 6A-L). The assay detected an age-dependent decline in axon extension, with 104.7 ± 16 axons reaching 2 mm from P3 explants, 24.4.7 ± 11 from P5, and 0.1 ± 0.1 from P7 (p < 0.01, ANOVA with post-hoc Sidak's) (Fig. 6M-O). Forced expression of KLF6 increased axon growth at all ages: 175 ± 11.4, 78.1 ± 23.3, and 8.0 ± 2.5 axons reached 2 mm from P3, P5, and P7 explants (Fig. 6M-O). In contrast, VP16-STAT3 treated neurons showed no increased axon extension in this assay. When VP16-STAT3 and KLF6 were combined they again produced a synergistic increase in axon extension in P5 and P7 explants (138.8 ± 25.0 and and 60.73 ± 12.14 axons reached 2 mm) (Fig. 6N,O). The combined KLF6/VP16-STAT3 was the only condition in which P7 axons ever extended as far as 3 mm (Fig. 6O). Thus combined expression of VP16-STAT3 and KLF6 produces synergistic enhancement of axon elongation that partially blocks the age-dependent reduction in axon extension in postnatal cortical neurons. (B,C). Postnatal cortical neurons were transfected with KLF6, caSTAT3, VP16-STAT3, or controls for equal plasmid load, along with EGFP reporter. After two days in culture, neurite lengths were measured by automated image analysis (Cellomics) in transfected neurons (EGFP+, green, arrows in B). As expected, KLF6 and VP16-caSTAT3 significantly increased neurite lengths (*p < 0.05, **p < 0.001, p < 0.0001, ANOVA with Sidak's multiple comparisons test, N > 200 cells from three replicate experiments). Combined expression of KLF6 with either caSTAT3 or with VP16-caSTAT3 increased neurite length above the level of KLF6, caSTAT3, or VP16-caSTAT3 alone, and above the sum of their individual effects, demonstrating synergistic increases in neurite length. Scale bar is 50 µm (B).
Scientific RepoRts | (2018) 8:12565 | DOI:10.1038/s41598-018-31101-5 Co-occupancy by KLF6 and STAT3 in axon growth gene networks. Transcription factors can functionally synergize through a variety of mechanisms, including physical binding, transcriptional activation of upstream regulators, and co-occupancy of regulatory DNA 7 . We considered physical binding, which has been shown between STAT3 and another KLF family member, KLF4 13 . KLF6 lacks the STAT3 binding region identified in KLF4, however, and no KLF6/STAT3 interactions are reported in any of the several interaction databases we queried (STRING, IntAct, MIPS, MINT, DIP, BioGrid). Regarding cross-activation, recent profiling in oligodendrocytes indicated that KLF6 activity may influence STAT3 by increasing expression of upstream regulators including gp-130 and JAKs 52 . Examination of our RNA-seq analysis did not detect upregulation of these genes by KLF6, arguing against a similar mechanism in neurons. Similarly, forced expression of STAT3 in cultured cortical neurons produced no significant change in KLF6 transcript levels (Fig. S6).
Without strong support for physical binding or transcriptional cross-regulation in our system, we therefore focused on potential co-occupancy of DNA, which was indicated by the prior motif analysis. We performed footprinting analysis of ATAC-seq data from early postnatal cortex, which enables genome-wide identification of loci that are stably bound by TFs 38,53 (Fig. 7A,B). This analysis found 692 genes whose promoter or intragenic enhancer regions contained deep footprints that correspond to the consensus binding sequence for KLF6/7 (see methods for a discussion of enhancer identification) (Fig. 7A,B). We performed network analysis of the corresponding genes, which revealed ten subnetworks. Again consistent with the role of KLF6 in axon growth, three prominent sub-networks were enriched for axon-relevant terms including axon guidance, axon projection, and cytoskeletal organization (Fig. 7D). We next performed footprint analysis to identify sites of STAT3 binding that occurred within the same promoter/enhancer regions as KLF6 footprints. Overall, 6.5% of regulatory regions with KLF6 footprints were also co-occupied by STAT3. For reference, the maximal rate of co-occupancy for any two TFs in this dataset was 12% (CTCF and YY1). To assess statistical significance, we analyzed a background set of 10 additional TFs that were effectively footprinted in the dataset. The average co-occupancy of KLF6 with TFs in the background set was 0.4% (±0.13% SD), and thus KLF6's co-occupancy with STAT3 is more than 16-fold higher (>50 standard deviations) than the average background rate. Thus, in early postnatal tissue, KLF6 and STAT3 display substantial co-occupancy across the genome at a rate higher than chance. In adult tissue, co-occupancy analysis revealed two main changes. First, as predicted by a reduced level of KLF6 expression, the absolute number of genes with KLF6 footprints declined by 60% with age, to 280. Second, the rate of STAT3 co-occupancy in these residual genes declined to 2%. Combined, the reduction in total amount of KLF6 binding and the lower rate of STAT2 co-binding yielded an 87% drop in KLF6/STAT3 co-occupied genes across age. Interestingly, closer analysis of co-occupancy in the postnatal tissue revealed that co-occupied genes were not spread evenly across the entire KLF6 network, but instead occurred in only the axon-relevant sub-networks (guidance, projection, and cytoskeleton) (Fig. 7D). Taken together, these data strongly indicate that during early postnatal periods, but not in the adult, STAT3 and KLF6 influence axon growth by co-occupying regulatory DNA in growth-relevant gene networks 7 .

Discussion
Transcription factors are emerging as powerful tools to enhance axon regeneration, but likely must be supplied in appropriate combinations for maximal effect. Efforts to identify optimal combinations are hampered by limited understanding of the transcriptional mechanisms that underlie TF-triggered growth. Here we identified KLF6 as an effective pro-regenerative TF, and then integrated transcriptional profiling with recent advances in gene Figure 7. TF footprinting confirmed genome-wide co-occupancy of KLF6/STAT3 in regulatory regions relevant to axon growth. ENCODE consortia ATAC-Seq datasets generated from PO mouse cortex were used for genome-wide TF footprinting to identify regulatory regions bound in vivo by KLF6 and STAT3. (A) Representative UCSC browser image of gene promoter region with clear notches within peaks of accessibility that correspond to KLF6 and STAT3 motifs confirming co-occupancy. (B) Aggregate ATAC-Seq footprint for KLF6/7 (shown motif) generated over binding sites within the genome. (C) %co-occupancy for KLF6 with STAT3 is 50 SDs above background level at P0 and 15 SDs above background in adult neurons (D) Regulatory network analysis of genes with deep KLF6/7 footprints revealed 10 sub-networks (a-i) enriched for distinct functional categories highly relevant to axon growth (see legend). Genes showing co-occupancy by both KLF6 and STAT3 were functionally clustered in sub-networks a,b,c (indicated by *red asterisk) Nodes correspond to target genes and edges to multiple modes of interaction (physical, shared upstream regulators, shared signaling pathways and inter-regulation). Node color represents their corresponding functions as denoted in the legends and node size is based on significance for enrichment in functional category. Only significantly enriched GO categories were included in the network analysis (p < 0.05, enrichment -right-sided hypergeometric test with Bonferroni correction). network analyses to identify KLF6-responsive gene modules. We then leveraged these insights to identify STAT3 as a TF that functionally synergizes with KLF6 in promoting axon extension and identified co-occupancy of the two factors in regulatory DNA associated with specific pro-regenerative gene networks. These results offer new insight into transcriptional networks that underlie axon extension and suggest combined KLF6/STAT3 activity as a first step in their reconstruction.
KLFs, a 17-member family of zinc finger transcription factors, have emerged as important regulators of the neuron-intrinsic capacity for axon growth 14,54 . We previously demonstrated that KLF7 increases corticospinal regeneration in vivo when delivered as a transcriptionally active and protein-stabilized mutant (VP16-KLF7), but is ineffective in its wildtype form 15 . Here we show that KLF6, like KLF7, is developmentally downregulated in cortical neurons, and that forced re-expression enhances CST regeneration after spinal injury. Unlike KLF7, however, wildtype KLF6 protein is readily detectable after viral overexpression, indicating that KLF6 protein is not subject to the same regulation as KLF7. Importantly, KLF6 and -7 very likely act through similar transcriptional targets. They share nearly identical DNA-binding domains, indicating similar genomic targeting, and functionally compensate for one another in promoting axon growth 54 . Indeed, we showed previously that forced co-expression of KLF7 and KLF6 produces no gain in axon growth above either alone, again indicating a shared mechanism 14 . Thus, because KLF6's effects match KLF7's but do not depend on an artificial VP16 transcription domain, it is likely a better candidate for eventual translation and is a preferable tool to clarify the endogenous transcriptional mechanisms of pro-regenerative TFs.
In recent years, transcriptome profiling approaches have revealed large numbers of genes associated with various forms of endogenous or TF-stimulated axon growth 9,55-57 . Extracting biological insights from lists of differentially expressed genes, however, remains a non-trivial challenge. One standard approach is to test for overall enrichment of ontological terms, which can provide a first-pass indication of the cellular functions that are impacted by the gene set. One the other hand, this approach leaves unanswered the question of how the gene products might interact, and risks missing functions that exist in subnetworks but which are diluted by whole-set analysis. Accordingly, here we applied a network-based approach in which genes were first clustered into sub-networks according to mutual interactions. This analysis revealed that KLF6-responsive genes broke into five prominent sub-networks. Interestingly, when ontological enrichment analysis was applied to each subnetwork, all five showed significant enrichment for functions that play important but complementary roles in axon growth 10,[58][59][60][61][62][63] . Similarly, network analysis of genes downregulated upon KLF6 over-expression revealed subnetworks that were highly enriched for synaptic plasticity and neurotransmission. This is consistent with the growing body of evidence that the developmental decline in axon growth capacity is due in part to a neuronal transition from axon extension to synaptogenesis 3,55,64 . These results highlight the ability of transcriptional profiling and network analysis to yield new insight into the complementary cellular networks activated by pro-regenerative treatments.
It is increasingly appreciated that manipulation of multiple TFs will likely be necessary to evoke a full regenerative response, and recent efforts have sought to identify optimal combinations 7,9,17,65 . TFs can interact through a variety of mechanisms, and the details of these interactions determine the optimal strategies for detecting them 7 . For example, prior efforts to identify regenerative TF combinations have focused on direct protein-protein interactions in order to build networks and identify hub TFs 9 . In the case of KLF6 and STAT3, however, existing protein-protein interaction databases do not indicate direct association. Although this does not preclude the possibility of physical interactions, it led us to explore alternative mechanisms. TFs commonly synergize by targeting overlapping sets of regulatory DNA, without direct physical TF-TF interaction 66 . Thus a molecular signature of functionally synergizing factors is frequent co-occupancy of promoter regions. Detection of co-occupancy has been leveraged to clarify regulatory TF networks in other biological processes, but has not previously been examined for axon growth [67][68][69][70][71][72] .
Therefore, starting with the KLF6-responsive gene modules, we analyzed upstream regulatory DNA to identify relevant TFs that may synergize with KLF6 to drive growth. One such factor was SP1, which was previously implicated by gene network analysis as a potential contributor to regeneration in peripheral neurons 9,73 . We found, however, that forced expression of SP1 had no effect on neurite outgrowth in postnatal cortical neurons, perhaps owing to differences in cell type. Notably, the recognition motif for STAT3 is also highly over-represented in the promoter region of genes that are upregulated by KLF6 overexpression, suggesting a potential interaction between the two factors. STAT3 has been previously shown to have pro-regenerative effects in retinal ganglion cells and CST neurons 51,[74][75][76] . Interestingly, transcriptionally active VP16-STAT3, but not wildtype, was recently shown to promote axon growth in postnatal cortical neurons. Consistent with this, we found that VP16-STAT3 produced a modest increase in axon growth in one of the two assays used here. When combined with KLF6, however, VP16-STAT3 produced a dramatic increase in neurite extension in two growth assays, including an explant-based assay of long distance axon growth. Thus, as predicted by bioinformatics modeling, KLF6 and STAT3 synergize to enhance axon growth in developing cortical neurons.
Finally, we took advantage of recent advances in genome-wide digital footprinting to substantiate the prediction that KLF6 and STAT3 co-occupy regulatory DNA of pro-regenerative genes. To our knowledge, this analysis created the first comprehensive in vivo binding map for pro-regenerative TFs during periods of robust axon extension. In line with the RNA-seq based network analysis, we found that KLF6-footprinted genes separated into several sub-networks, some of which were highly enriched in functions relevant to axon growth including neuron projection extension, motor neuron axon guidance, and regulation of neuron projection extension. Importantly, instances of co-occupancy by STAT3 were not distributed randomly across the entire network, but instead occurred exclusively in the growth-relevant subnetworks. This finding lends strong support to the model that KLF6 and STAT3 synergistically enhance axon growth through a mechanism of co-occupancy of regulatory DNA. More generally, this analysis highlights the utility of a comprehensive systems approach to identify cooperative TF regulators of axon growth. This approach offers a rational framework to prioritize factors for future tests of combined gene function in vivo, for example experiments to test whether combined expression of KLF6/ STAT3 in CST neurons produces improvements in axon regeneration or functional recovery above the level of KLF6 alone. Starting from genome-wide expression datasets, integrated network and footprinting analyses offer enormous potential to further elucidate TF regulatory circuits that underlie axon growth. Data availability. All data generated or analyzed during this study are included in this published article (and its Supplementary Information files). RNA-Seq raw files are available from the corresponding author on request (currently undergoing submission in NCBI GEO-SUB2954799).