Corticosterone-mediated regulation and functions of miR-218-5p in rat brain

Chronic stress is one of the key precipitating factors in major depressive disorder (MDD). Stress associated studies have underscored the mechanistic role of epigenetic master players like microRNAs (miRNAs) in depression pathophysiology at both preclinical and clinical levels. Previously, we had reported changes in miR-218-5p expression in response to corticosterone (CORT) induced chronic stress. MiR-218-5p was one of the most significantly induced miRNAs in the prefrontal cortex (PFC) of rats under chronic stress. In the present report, we have investigated how chronic CORT exposure mechanistically affected miR-218-5p expression in the rat brain and how miR-218 could trigger molecular changes on its downstream regulatory pathways. Elevated expression of miR-218-5p was found in the PFC of CORT-treated rats. A glucocorticoid receptor (GR) targeted Chromatin-Immunoprecipitation (ChIP) assay revealed high GR occupancy on the promoter region of Slit3 gene hosting miR-218-2 in its 3rd intron. RNA-sequencing data based on RNA Induced silencing Complex Immunoprecipitation (RISC-IP) with AGO2 in SH-SY5Y cells detected six consistent target genes of miR-218-5p (APOL4, DTWD1, BNIP1, METTL22, SNAPC1, and HDAC6). The expression of all five genes, except APOL4, was successfully validated with qPCR in CORT-treated rat PFC. Further, Hdac6-based ChIP-seq experiment helped in mapping major genomic loci enriched for intergenic regions in the PFC of CORT-treated rat. A proximity-based gene ontology (GO) analysis revealed a majority of the intergenic sites to be part of key genes implicated in central nervous system functions, notably synapse organization, neuron projection morphogenesis, and axonogenesis. Our results suggest that the upregulation of miR-218-5p in PFC of CORT-treated rats possibly resulted from GR biding in the promoter region of Slit3 gene. Interestingly, Hdac6 was one of the consistent target genes potentially found to regulate CNS related genes by chromatin modification. Collectively, these findings establish the role of miR-218-5p in chronic stress and the epigenetic function it plays to induce chromatin-based transcriptional changes of several CNS genes in triggering stress-induced disorders, including depression. This also opens up the scope to understand the role of miR-218-5p as a potential target for noncoding RNA therapeutics in clinical depression.


Methods
Animals. All the experiments were carried out in accordance with the National Institutes of Health guide for the care and use of laboratory animals. The Animal Care Committee of the University of Alabama at Birmingham approved this study. The study is reported in accordance with ARRIVE guidelines (https:// arriv eguid elines. org).
Virus-free Sprague-Dawley male rats, initially weighing 220-250 g, were used. Rats were housed in groups of three under standard laboratory conditions (temperature 21 ± 1 °C, humidity 55 ± 5%, 12-h light/dark cycle). Animals were provided free access to food. Rats were acclimatized for 1 week before the start of experiments. Six animals in each group were used for miRNA expression analysis, 5-6 animals in each group were used for chromatin immunoprecipitation (ChIP) based qPCR experiments and target gene expression profiling were done in 13-22 animal in each group.
Corticosterone treatment. The procedure for corticosterone treatment has been described in our earlier publication 22 . Rats, under light halothane anesthesia, were subcutaneously implanted with corticosterone pellets containing 100 mg of corticosterone in a cholesterol base (Innovative Research of America, Sarasota, FL, USA). These corticosterone pellets can maintain a physiological serum concentration of corticosterone for 21 days. The release of CORT after implantation of a 100-mg corticosterone pellet is 4.76 mg per day. Control rats underwent an identical surgery procedure with implantation of a cholesterol base pellet. Rats were decapitated 21 days after pellet implantation. This procedure has been published in our earlier study 22 .
The trunk blood was collected on ice and centrifuged; then, the serum was stored at − 80 °C until the assays were performed. Serum corticosterone levels were measured by a commercially available radioimmunoassay kit (ICN Biomedical, Inc., Cleveland, OH, USA). Brains were removed quickly after the blood was taken. Prefrontal cortex (PFC) was dissected out and immediately stored at − 80 °C until analyzed. Rats were decapitated between 09:00 h and 11:00 h, corresponding to 3-5 h after lights on. All experiments were performed in 6-12 rats per group.

RNA-immunoprecipitation assay (RNP-IP).
Cell lysates collected in RIP lysis buffer were used for RNA-immunoprecipitation assay following a protocol by 18,23 . Briefly, prewashed protein A/G beads (Santa Cruz Biotechnology, Dallas, TX, USA) were conjugated with Anti-Ago2 monoclonal antibody (MABE253, Millipore, Billerica, MA, USA) and were added with cell lysate. Following overnight incubation, the immunoprecipitated RNA-induced silencing complex (RISC) was used to extract overexpressed miR-218-5p and its bound targets mRNAs using TRIzol reagent as RNA-IP samples. Raw data files in FASTQ format were generated from the Illumina sequencer. To examine the sequencing quality, the quality score plot of each sample was plotted. Sequence quality was examined using the FastQC software. After quality control, the fragments were 5′, 3′-adaptor trimmed and filtered ≤ 20 bp reads with Cutadapt software. The trimmed reads were aligned to reference genome with Hisat 2 software. Based on alignment statistical analysis (mapping ratio, rRNA/mtRNA content, fragment sequence bias), it was determined whether the results could be used for subsequent data analysis. The expression levels (FPKM value) of known genes and transcripts were calculated using ballgown through the transcript abundances estimated with StringTie. The number of identified genes and transcripts per group was calculated based on the mean of FPKM in group (≥ 0.5). Differentially expressed gene and transcript analyses were performed with R package ballgown. The volcano plot was created based on FPKM values with web-based tools (https:// paolo. shiny apps. io/ Shiny Volca noPlot/). The heat map was also created based on FPKM value by ClustVis 25 .

Identification of consistent target genes of miR-218-5p based on input RNA-seq and IP
RNA-seq results. The prominent target genes of miR-218-5p were determined by the following criteria: (1) significantly downregulated genes in input RNA-seq, and (2) significantly upregulated genes in RISC-IPed RNA-seq.
ChIP-sequencing (ChIP-Seq) and ChIP-qPCR (qChIP). Briefly, 30 mg of frozen PFC tissue was cut into small pieces. The 1% formaldehyde in ice-cold phosphate buffered saline (PBS) containing proteasomal inhibitors (1X complete protease inhibitor, 1 mM of PMSF [phenylmethylsulfonyl fluoride] and 25 μm MG-132) was added to brain pieces. The brain pieces were incubated for 10 min at 37 °C for cross-linking reaction. The cross-linking reaction was stopped by adding 125 mM glycine. Crosslinked tissue pieces were washed twice with ice-cold 1XPBS and suspended with SDS lysis buffer (1% SDS, 10 mM EDTA, and 50 mM Tris-HCl [pH 8.1]) containing protease and proteasome inhibitors for 10 min on ice bath. Brain lysates were sonicated using Active Motif Q800R sonicator (Amplitude: 25%, Pulse Duration: 30 s, pulse on time: 10 min, number of pulses: 20). The resulting homogenates were centrifuged for 10 min at 13,000×g at 4 °C to remove debris. The supernatants were diluted in ChIP dilution buffer with protease and proteasome inhibitors (0.01% SDS, 1.1% Triton X-100, 1.  ChIP with HDAC6 antibody and corresponding input DNA samples were subjected to ChIP-seq. ChIP-seq libraries were sequenced with the Illumina HiSeq2000 sequencer (50-nucleotide pair-ended read) at the Heflin Center Genomics Core, University of Alabama at Birmingham. On the other hand, immunoprecipitated ChIP DNA with GR antibody and corresponding input DNA samples were subjected to relative quantification with EvaGreen dye-based chemistry. The control qPCR experiment was done using DNA from IgG antibody-based ChIP-IP and corresponding ChIP-input samples. All the primers are shown in Table S1. The input normalized fold enrichment was calculated according to the method described earlier 26 .
Functional annotation of nearest gene from peak change in ChIP-seq. The nearest genes from the significantly altered peaks in ChIP-seq in each sample were used in standalone Cytoscape program to perform gene ontology (GO) analysis using ClueGO plugin. Display pathways were selected with p ≤ 0.05. Clustering was done based on the common functionality of genes enriched for specific term. The kappa score was set at 0.5 to define the term-term relationship (i.e., edges between the nodes).
qPCR based mRNA and miRNA-specific gene expression. mRNA-specific cDNA was synthesized with the Oligo (dT) 18 priming method. A mixture of 0.5 µg of total RNA, oligo (dT) 18 , dNTP, and double distilled water was incubated at 65 °C for 5 min and quickly chilled on the ice. Subsequently, the reaction mixture was added with 1× 1st strand synthesis buffer, 0.01 mM DTT, 2 U RNaseOut and 200 U M-MLV reverse transcriptase (Invitrogen, Grand Island, NY, USA) and incubated at 37 °C for 50 min. Finally, the reaction was inactivated at 70 °C. miRNA-specific cDNA synthesis was performed using poly(A)-tailing method and priming with oligo dT adaptor primer. qPCR-based relative expression was measured with EvaGreen chemistry (Applied Biological Materials, Canada). All the primer sequences are shown in Table S1. Gapdh was used as a normalizer for mRNA expression, whereas U6 was used for miRNA expression. Fold expression changes were calculated by Livak's ΔΔCt method 27 . qPCR based miR-218 expression assay using in vitro cellular model with GR gain of function mutation. We used HEK-293 cell line to transfect GR overexpression clone (SinoBiologicals Inc., Wayne, PA, USA). In parallel, control group was prepared by transfecting the cell line with scramble plasmid. Fortyeight hours post-transfection, cells were harvested, and RNA was extracted as mentioned above. The extracted RNA was used in miRNA-specific cDNA preparation and subsequently used in qPCR with primers to determine miR-218 expression. The expression was normalized by U6 gene as reference. The fold-change was determined as mentioned above.
Statistical analysis. The SPSS 25.0 software (IBM, USA) was used for statistical analyses. Shapiro-Wilk test was used to assess the normality of the data. The average differences of mRNA and miRNA expressions were assessed by the student's t test. The difference of fold enrichment in ChIP experiment was also assessed by the student's t test. Statistical significance was set at the 95% level (p ≤ 0.05).

Results
Serum corticosterone levels. Serum corticosterone levels were measured in sham and corticosterone-treated rats. The serum corticosterone levels in these rats were as follows: sham, 40.3 ± 0.10. 3; CORT, 109.8 ± 27.1 ng/ml, which was significantly higher in the CORT-treated group. As with our earlier studies 22,28 , corticosterone level was significantly higher (p < 0.001) in corticosterone pellet-implanted rats.

miR-218-5p expression in CORT-treated rats.
To check miR-218-5p expression in CORT-treated rats, qPCR was performed in the PFC. Expression of U6 was not significantly different between control and CORTtreated groups. The miR-218a-5p expression was significantly elevated in the CORT-treated group (p = 0.031, Fig. 1A).
Input samples-based RNA-seq of coding transcripts and core analysis of downregulated genes. A total of 12,222 genes were detected in RNA-seq for input samples. The volcano plot was created using all the differentially regulated genes based on input RNA-seq and Ago-RISC IP-RNA-seq results (Fig. 2B,C). Based on input RNA-seq results (Fig. 2B), 378 genes were found significantly upregulated and 379 were significantly downregulated in the mimic group compared to the vehicle group (downregulated genes are shown in Table S2). The 379 downregulated genes were further used for subsequent analysis since miRNAs generally regulate mRNA expression inversely. In Fig. 2C, volcano plot showing differential enrichment (up) and depletion (down) of target genes in RISC pull down (IP) complex, 560 genes were found up and 125 genes were down in IP group.
RNA-seq for IP samples to identify consistent target genes based on core analysis with the input RNA-seq and IP RNA-seq data. A total of 12,522 genes were detected in RNA-seq for IP sam- www.nature.com/scientificreports/ ples. The volcano plot was created using all the genes (Fig. 2C). Out of 12,522 genes, 560 were significantly upregulated and 125 were significantly downregulated in the mimic group compared to the vehicle group (all upregulated genes are shown in Table S3). We focused on 560 upregulated genes for subsequent analysis because those genes should be packed in RISC complex under the overexpressed miR-218-5p condition and will be subsequently downregulated. The result of the putative target genes, that were downregulated in input RNA-seq and upregulated in IP RNA-seq experiments, was used to determine an overlapped set of genes. The overlapped six genes (APOL4, DTWD1, BNIP1, METTL22, SNAPC1, and HDAC6) were predicted to be experimentally consistent target genes. In the chord diagram, the degree of interactions between input and IP groups on the six targets are shown with respective fold-change values (Fig. 2D). Two expression heat maps (Fig. 2E,F) were created based on six genes as they were originally identified in IP-and input-RNA-seq from miR-218 transfected human neuroblastoma cell line (SHSY-5Y). Subsequently, five genes, except APOL4 were tested in the qPCR experiments in rat PFC because APOL4 gene doesn't exist in rats (https:// www. ncbi. nlm. nih. gov/ gene/).

Validation of consistent target genes in CORT-treated rats based on qPCR.
The qPCR of five consistent target genes was done to investigate whether these genes are altered in PFC of CORT-treated rats. The level of Gapdh (control vs CORT = 22 vs 13, t = -2.005, df = 33, p = 0.053) was not significantly different between control and CORT-treated groups. As shown in Fig. 3, all 4 consistent target genes were significantly downregu-

Mapping the genomic region bound by Hdac6 revealed by ChIP-seq. ChIP-seq experiments were
performed to reveal the genome-wide Hdac6 binding regions in PFC of sham and CORT-treated groups. As shown in Fig. 4A, the ChIP-seq experiment revealed genome-wide binding of Hdac6 largely located in the intergenic region (25/455, 94.5%); rest of them were in the intronic region (25/455, 5.5%). In Fig. 4B, Hdac6 associated peaks and their distribution, proximal to respective genes, have been shown with ridgeline distribution plot. Considering the p-values, the statistically highest known binding motif was common to ZNF264 (Fig. 4C).

GO analysis based on nearest genes from the region of significant peak change. All 455 genes
identified from ChIP-Seq results were subjected to Cytoscape-ClueGO for biological processing (BP) analysis and creating network based on significant BP terms. As shown in Table S4, multiple CNS related BP terms (e.g., synapse organization, dendrite development, morphogenesis, axonogenesis, and neuron projection morphogenesis) were revealed as significant GO terms. These terms had high connectivity among themselves (Fig. 4D).  www.nature.com/scientificreports/ onstrated consistence increase in miR-218 expression (20%) induced by GR as compared to scramble control (Fig. S1).

Discussion
In this study, we found that CORT enhanced the gene transcription of miR-218-5p, which appeared to be through the biding of GRs on the promoter region of Slit3 gene. Using RNA-IP and input RNA-seq data from miR-218-5p overexpressed SH-SY5Y cells, six genes (APOL4, DTWD1, BNIP1, METLL22, SNAPC1, and HDAC6) were identified as consistent target genes of miR-218-5p. Five genes except APOL4 (not expressed in rats) were downregulated in PFC of CORT-treated rats. Of these, ChIP-seq data in PFC of CORT-treated rats revealed that the nearest genes from Hdac6-associated peaks were associated with several CNS related biological processes terms. Glucocorticoid receptors are known to activate the transcription of early response genes associated with stress 5,6 . Our aim was to understand the mechanism of upregulated transcriptional activity of miR-218a by CORT. Our results show that mir-218 could possibly be co-activated by CORT, given that miR-218 gene is located in the intron of Slit3 gene. ChIP-seq experiment using GR antibody indicated significant GR binding sites on Slit3 promoter region of CORT-treated rats. We found a trend towards the increased binding of GR on both miR-218a-1 and -2; however, it did not reach statistical significance. This is not surprising given that there are three predicted GR biding sites on Slit3 promoter regions compared to two predicted GR biding sites on other Slit genes. Our in vitro GR overexpression model also supplemented this finding where we found an increased miR-218 expression under exogenous gain of function manipulation of NR3C1 gene following a transfection of cDNA overexpression clone in HEK-293 cell line. It is intriguing that GR mediated gene promoter activity largely depends on interaction with coregulatory molecules. In fact, depending on the nature of the coregulator, GR can either activate or inhibit promoter specific gene transcription. Studies in the past have reported allosteric www.nature.com/scientificreports/ modulation of GR by transcriptional coregulators, leading to chromatin destabilization and subsequent activation of target gene promoter 29 . Although our ChIP experiment results are not supplemented with coactivator experiments, the present findings related to GR-induced miR-218 expression can be well articulated with the allosteric modulation theory of GR. Despite the mixed findings on stress-induced miR-218 expression changes, our report for the first time elaborated the direct role of GR in transcriptional activation of miR-218 by Slit3 gene promoter. This has also, for the first time, helped in understanding the co-transcriptional activation of an intronic miRNA via a shared promoter usage. Given the role of Slit3 gene in axonal guidance and genetic predisposition to major depressive disorder due its duplicated locus on 5q35.1, it is possible that co-transcriptional activation of Slit3 and miR-218 may be necessary to develop depressive phenotypes under chronic stress 30 .
Since induced transcriptional activity of miRNAs is directly associated with posttranscriptional target gene silencing, we investigated the direct repression of targetome of miR-218-5p using Ago2 mediated RISC-sequencing experiment. There are eight AGO family members, and AGO1-4 are capable of loading miRNA and endonuclease activity. However, RNA interference-dependent gene silencing exclusively belongs to AGO2, which is one of the essential components of RISC. Functionally, mature miRNAs are incorporated into RISC, which helps them regulate the fate of target gene expression by pairing primarily to the 3'-UTR of protein-coding mRNAs. In this context, high throughput RISC profiling is critical to ascertain the role of specific miRNA, in this case miR-218, to determine its regulatory potential on cellular targets critical in stress responsive pathways. To find consistent target genes of miR-218-5p, we conducted RNA-seq with RNA-input samples and found six genes as miR-218-5p consistent target genes. Subsequently, we confirmed the downregulation of five genes (Dtwd1, Bnip1, Mettl22, Snapc1, and Hdac6) by qPCR experiment in PFC of CORT-treated rats. DTWD1 (DTW domain containing 1) is found in STAR*D study as a candidate gene associated with the side effect of antidepressant 31 . Budde et al. reported that rs586758 located in DTWD1 is related to the Global Assessment of Functioning (GAF; DSM-IV Axis V) 32 . Another target gene, BNIP1, is implicated in the functions of the endoplasmic reticulum (ER), which is primarily responsible for cellular stress responses 33 . Under a stressful condition, protein folding is disrupted 34 . ER-mediated unfolded protein response (UPR) is activated to re-establish the protein homeostasis. We have earlier reported the dysregulated expressions of UPR genes in postmortem PFC of MDD subjects 35 and rat brain which showed depression phenotype 36,37 . The abnormality of UPR genes has also reported in PFC of schizophrenia subjects 38 . METTL22, a member of methyltransferase like protein family, is responsible for the methylation of DNA-RNA-binding protein Kin17 and is implicated in DNA repair and replication and mRNA processing 39 . SNAPC1, one of the subunits of small nuclear RNA (snRNA) activating protein complex, acts as a transcription factor to mediate transcription of snRNAs 40,41 . Interestingly, SNAPC1 acts as a transcriptional coactivator facilitating RNA polymerase II elongating the transcripts 42 .
One of the miR-218-5p target genes is HDAC6, which belongs to epigenetic family of chromatin modifiers. HDAC (histone deacetylase) is a key enzyme that regulates the reversible acetylation of histones [43][44][45] . Out of four subtypes in HDAC, HDAC6 belongs to Class II. As a direct target of miR-218-5p, Hdac6 can potentially alter genome-wide chromatin conformation of many stress-responsive genes under chronic stress. It would be interesting to know if miR-218 as a posttranscriptional gene regulator can provide an additional layer of epigenetic masking on HDAC6 expression. Our ChIP-seq data obtained from in vitro cell culture model showed that HDAC6 peaks were associated with genes related to key CNS functions including axonal morphogenesis and synaptic functions. Since glucocorticoid signaling system is partly regulated by reversible acetylation-deacetylation level of HSP90 chaperonin molecule 46 , there could be a possibility to see abnormalities in stress responsiveness mediated by altered HDAC6 activity. We strongly anticipate that reduced level of HDAC6 in PFC of rats (as we found in our CORT treated rats) could be part of an epigenetic compensatory mechanism. At the molecular level, this compensation could be attributed to the induced expression of miR-218-5p to directly check the deacetylation status of HSP90 and facilitating other HDAC6 target proteins to achieve a hyperacetylation level under chronic stress 47 . In our ChIP-seq results, we found that the majority of Hdac6 peaks were localized in the intergenic regions (94.5%), as previously reported for other class II HDACs such as HDAC4 48 , HDAC5 49 , HDAC7 50,51 , and HDAC9 48 . The gene ontology results indicated that nearest genes from Hdac6 peak were associated with several CNS related biological processes terms. It has earlier been reported that HDAC6 regulates neural migration and dendrite development 52 , as well as modulates synaptic plasticity and neuronal differentiation in human stem-cell-derived neurons 53 . Especially in PFC of acute stressed rats, HDAC6 regulate the synaptic efficacy 54 . The decreased expression of HDAC6 has also been found in the peripheral blood of bipolar patients 55 . Importantly, HDAC6 inhibitor potentially has antidepressant effect and Hdac6-deficient mice exhibit less anxiety and antidepressant-like behavior 56 .
Altogether, we found that miR-218-5p, upregulated in brain of CORT-treated rats, was GR-mediated. This upregulation, and consequent downregulation of target genes, could potentially be involved in the development of depressive phenotypes in stressed rats. We found 6 consistent target genes of miR-218-5p that were detected by in-vitro experiments in which miR-218-5p was overexpressed followed by RNA-seq. Of them, the nearest genes from HDAC6 peak, identified via ChIP-seq data, were associated with key CNS functions including axonal morphogenesis and synaptic functions. Our ChIP-seq results also identified several HDAC6-associated peaks across the genomic loci, which largely show the involvement of this intragenic miR-218-5p in the overall modification of chromatin sites. Even though the subcellular distribution of HDAC6 is cytoplasmic, it can also shuttle back to the nucleus and may interact with many nuclear proteins 57 . This has been seen as a tight control on chromatin accessibility outlining the importance of reversible chromatin modifications of key stress-associated genes. Together, these findings provide regulatory insights of miR-218-5p, which could possibly be associated with stress-induced depression. It also provides possible feature therapeutic target for the treatment of depression. Further studies in animals overexpressing miR-218, and examining consequent behavioral phenotype, will be key to several of these answers.