Transcriptome Analysis Reveals Neuroprotective aspects of Human Reactive Astrocytes induced by Interleukin 1β

Reactive astrogliosis is a critical process in neuropathological conditions and neurotrauma. Although it has been suggested that it confers neuroprotective effects, the exact genomic mechanism has not been explored. The prevailing dogma of the role of astrogliosis in inhibition of axonal regeneration has been challenged by recent findings in rodent model’s spinal cord injury, demonstrating its neuroprotection and axonal regeneration properties. We examined whether their neuroprotective and axonal regeneration potentials can be identify in human spinal cord reactive astrocytes in vitro. Here, reactive astrogliosis was induced with IL1β. Within 24 hours of IL1β induction, astrocytes acquired reactive characteristics. Transcriptome analysis of over 40000 transcripts of genes and analysis with PFSnet subnetwork revealed upregulation of chemokines and axonal permissive factors including FGF2, BDNF, and NGF. In addition, most genes regulating axonal inhibitory molecules, including ROBO1 and ROBO2 were downregulated. There was no increase in the gene expression of “Chondroitin Sulfate Proteoglycans” (CSPGs’) clusters. This suggests that reactive astrocytes may not be the main CSPG contributory factor in glial scar. PFSnet analysis also indicated an upregulation of “Axonal Guidance Signaling” pathway. Our result suggests that human spinal cord reactive astrocytes is potentially neuroprotective at an early onset of reactive astrogliosis.

Reactive astrogliosis is a common response of astrocytes to most Central Nervous System (CNS) injury 1,2 . The prevailing dogma of the involvement of reactive astrocytes in inhibition of axonal regeneration 3 , has been challenged by recent findings of Anderson et al. work, which demonstrated the critical role of reactive astrocytes in aiding axonal regeneration in rodents 4 . While some of the previous reports have identified reactive astrocytes as neuroprotective agents 1,5,6 , their direct neuroprotective effects were not well-documented in the human CNS astrocytes. As astrocytes' functions differ among regions in the CNS 7 , we investigated the presence of axonal regenerating and neuroprotective properties of reactive astrocytes 4 in human spinal cord derived astrocytes. Elucidating these properties at an early onset of reactive astrogliosis, will enable physicians and scientists to design novel therapeutic strategies exploiting the potential of reactive astrocyte mediated endogenous recovery.
It has been reported that among inflammatory mediators such as IL-1, TNF, IFN, and TGF, which are known to be active post-CNS injury, cytokine interleukin-1β (IL1β) is specifically critical for induction of reactive-astrocyte phenotype 8 . IL1β is a prominent inflammatory cytokine and is an early regulator of astrogliosis 9 . IL1β mRNA is The surface area was significantly reduced (p = 0.00863) in reactive astrocytes as compared to control (1200 astrocytes counted for control and treated groups each). (E) The number of processes to cell ratio was significantly increased in reactive astrocytes (p = 9.09 × 10 −6 , 500 cells were counted each in control and treated groups) after induced by IL1β. (F) The processes length was significantly increased in reactive astrocyte as compared to control (p = 0.009, 400 processes were counted for experiment and control each). Scale bar 50 μm.
in IL1β treated astrocytes (Fig. 1D). This change in the surface area was due to the fact that astrocytes acquired a more polarized morphology with extensive processes from the cell bodies. As reported in Fig. 1E, the number of processes to cell ratio for reactive astrocytes (0.25 ± 0) was increased in comparison to control group (0.16 ± 0). Although, a small fraction of control astrocytes displayed extensive processes, their lengths (84.6 ± 5 µm; p = 0.009) were shorter as compared to IL1β treated astrocytes (111.9 ± 5 µm, 400 total processes counted for treated and experiment each) (Fig. 1F).
Axonal guidance molecules and neurotrophic factors genes involved in neuroprotection and axonal regeneration. One of the aims of this study was to identify genes regulating axonal guidance molecules and neurotrophic factors expression that promote neuroprotection, neurogenesis and axonal regeneration in human spinal cord reactive astrocytes. Figure 3 shows the fold changes of these specific genes 4 . It is noteworthy that 9 out of 25 axonal growth permissive genes were upregulated in reactive astrocytes [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41] . Of those, Fibroblast Growth Factor 2 (FGF2) is the most upregulated (3.24 folds), while matrilin2 (MATN2) (−2.37 folds) is the most downregulated axonal permissive genes. On the other hand, Slit Guidance Ligand (SLIT2) (2.54 folds) and Dorsal Inhibitory Axon Guidance Protein (DRAXIN) (2.52 folds) were the most upregulated genes involves in axonal growth inhibitory molecules. 8 out of 13 genes regulating axonal growth inhibitory molecules were down regulated, with Roundabout Guidance Receptor 2 (ROBO2) being the most downregulated genes in reactive astrocytes (−2.5 folds). Additionally, various matrix metallopeptidase and hyaluronan synthases were upregulated  (Supp. File 1) as well. Interestingly, CSPGs clusters were unchanged in reactive astrocyte transcriptome as compared to control (Supp. File 1) [42][43][44][45] . Differential Expressed Genes (DEGs) were then analyzed with KEGG human pathways database to reveal affected biological pathways in human spinal cord reactive astrocytes. Using subnetwork analysis (PFSnet) within each pathway, it is revealed that reactive astrocytes differentially regulates subnetworks of various pathways, as compared to nascent astrocytes ( Fig. 4 and Tables 3-4) 46 . Consistent with the observable changes in morphology of the in vitro reactive astrogliosis (Figs 1-2), PFSnet analysis revealed that actin cytoskeleton signaling pathway is one of the most altered pathway. As numerous pathways were differentially regulated in reactive astrocyte, our aim was to focus on main pathways that have critical role in axonal growth and development. We found that reactive astrocytes affect two subnetworks involved in axonal attraction and repulsion; FYN & RASGAP (Ras GTPase-activating protein 1) and SFK (Src family tyrosine kinase). In the "Axonal Guidance Signaling", the axonal attraction and outgrowth seems to be regulated mainly by FYN and RASGAP (Supp. Figure 3). Another differential component of axonal attraction appeared to be directly regulated by SFK. In addition, a closely linked pathway to axonal development is the "ECM-receptor interaction" (Supp. Figure 4). Overall, the "ECM-Receptor Interaction" is downregulated in reactive astrocytes. Our result indicates that, fibronectin and laminin were downregulated, while collagen, Thrombospondin (THBS), and tenascin were upregulated.

Discussion
In this study, we report the genome wide transcriptome profile of IL1β induced human spinal cord reactive astrocytes. We have studied a comprehensive subnetwork analysis based transcriptome pattern, which reveals the genes and their associated molecular function along with the important signaling pathways that regulate reactive astrogliosis. This provides a molecular roadmap to facilitate intervention and to optimize the benefits of reactive astrogliosis by either augmenting subnetworks of axonal attraction in the "Axonal Guidance Signaling" pathway or inhibiting signaling cascades in the "Axonal Repulsion" subnetwork.
CXCL5 and NMES1 are among the most upregulated genes by 205 and 108 folds respectively. The association of NMES1 with human spinal cord reactive astrocytes is not well-known currently. We speculate that a very likely role of C15orf48 or commonly known as NMES1 may be in the regulation of the cell-cycle process, as its downregulation has been associated with potential tumorigenesis. Similarly, rat model of reactive-astrocyte using lipopolysaccharide (LPS, chemical injury) and middle carotid arterial occlusion (MCAO, ischemic injury) to induce injury did not identify CXCL5 expression in their clusters of differentially expressed genes. Hence, we deduce that reactive astrocytes are a heterogeneous cell population whose molecular identity is based on the type of injury and origin of host cells 13,14,[47][48][49] .
The changes in the "Actin Cytoskeleton Signaling" pathway can be attributed to the change of the cytoskeleton into the more diffuse and ring-structure actin filaments 50 , leading to morphology changes in reactive astrocytes 14 . The reorganization of the actin cytoskeleton is affected by inflammatory cytokines and is closely related to the "Focal Adhesion" and "ECM Receptor Interaction" pathways. A component of MAPK signaling pathway, ERK phosphorylation has been shown to be increased after SCI 51 . However, our PFSnet subnetwork analysis suggest that for the "MAPK signaling pathway", the "JNK and p38 MAPK" sub-networks is activated in reactive astrocytes, while the classical "MAPK kinase sub-pathway" is down regulated. This selective subnetwork activation has not been clearly defined previously in reactive astrogliosis. Although the TGF-β signaling pathway is known to    be involved in astrocyte reactivity, its exact signaling cascade has not been well-established in reactive astrocytes.
Our findings indicate that one subnetwork of the pathway, involving ERK is upregulated while the other subnetwork involving the RhoA and PP2A is down regulated. Subnetwork analysis offers the advantage of differential view within a biological pathway, rather than revealing only an overall general alteration of biological pathway. One intriguing observation was an increase in neurotrophic factors and "Axonal Guidance Signaling" pathway, particularly brain derived neurotrophic factor (BDNF) in human spinal cord reactive astrocytes. The other neurotrophic factor that was upregulated was nerve growth factor (NGF). These two factors are critical in development, survival and regeneration of axons. Our findings thus suggest that early onset of reactive astrocytes could potentially promote neuronal and axonal development. Reactive astrocytes in the glial scar thus have the potential to secrete multiple axon-growth-permissive molecules, thus improving the microenvironment at and around the epicenter of injury. This will facilitate neuronal repair and regeneration. Recent findings also suggest that BDNF regulates the development of oligodendrocyte precursor cells 4,52-56 . Our results suggest that the primary axonal attraction aspect of reactive astrocytes, could be attributed to RasGAP and FYN over expression. Additionally, we also report the upregulated of "VEGF signaling" pathway, which has neurotrophic and neuroprotective effects on neuronal cells 57 .
We also identified that there was an absence of up regulation of CSPGs genes. It has been known that CSPGs are major axonal growth inhibitor in glial scar. Our results present the possibility that reactive astrocytes may not be the main contributor of CPSGs at the epicenter and around the injury. Furthermore, we also reported the upregulation of hyaluronan synthases, which are enzymes for producing hyaluronan. The overexpression of hyaluronan is known to improve SCI recovery by reducing the lesion, and pro-inflammatory cytokines suggesting another neuroprotective properties of reactive astrocytes.
Modification of extracellular matrix is an important element in formation of glial scar and modulation of axonal growth, regeneration and neuronal development post-CNS injury. MMP3, MMP1, and MMP12 belonging to the cluster of matrix proteins were found to be upregulated in reactive astrocytes. These genes have been associated with recruitment and migration of nascent astrocytes to the site of reactive astrocyte. Furthermore, MMP3 and MMP12 overexpression has been associated with remyelination. However, an increase in expression of matrix metalloproteases is known to enhance brain blood barrier (BBB) permeability and immune cells infiltration, resulting in inflammation post-CNS injury 58 . Our results indicate the upregulation of collagen and tenascin as well as downregulation of laminin and fibronectin in the ECM-receptor interaction pathways. In particular, collagen is critical in the early phase of tissue repair 59 and the role of tenascin, in promoting neural outgrowth is controversial. Laminin and fibronectin are well-known for promoting neurite outgrowth. Collectively, we interpret that the ECM changes in reactive astrocytes may potentially induce both neuro-permissive and inhibitory effects to modulate axonal regeneration and growth [59][60][61][62][63][64][65][66] .
In our report, we presented a fetal human spinal cord reactive astrogliosis induced by IL1β. Although the genomic profiles of fetal astrocytes and adult astrocytes are comparable, notable differences in expression has been identified from previous literature 67 . This include higher expression of pro-inflammatory miRNAs expression in adult astrocytes as compared to fetal 68 , while lower expression of fetal germinal matrix miRNAs in adult astrocytes as compared to fetal 67 . The higher expression of matrix associated miRNA could reflect the migrating status of the developing astrocytes in fetal CNS. Furthermore, young astrocytes also expressed genes involved in "neuronal differentiation" 69 . As our study is based on reactive astrogliosis in fetal astrocytes, it remains to be explored if neuroprotective and axonal regeneration potentials can replicated in adult reactive astrocytes. Also, our analysis indicates that both neuroprotective and inhibitory genes are activated in reactive astrocytes post 24 hours of IL1β exposure, suggesting that a balance of these gene expressions may regulate the functionality of reactive astrocytes.

Materials and Methods
Cell Culture. Human astrocytes derived from the spinal cord (19 weeks old fetus) was purchased from ScienCell Research Laboratories (Cat. No: 1820) 3 . The cells were cultured in Astrocytes Basal Medium (ScienCell TM ) supplemented with 2% FBS (ScienCell TM ), 1% Astrocyte growth supplement (ScienCell TM ) and 1% of penicillin-streptomycin (ScienCell TM ). The human astrocytes cells were plated on Matrigel (Corning) coated polystyrene 35mm plates with cell density of 1 × 10 5 cells. 24 hours after plating, cells were washed with PBS (HyClone ™ ). Fresh medium containing IL1β (Invitrogen TM ) at 100 ng/ml concentration was added and incubated for 24 hours in 5% CO2 at 37 °C.   Microarray. (i) All RNA samples were quality checked prior to microarray analysis in accordance to Affymetrix recommended protocols. Affymetrix 3′ IVT PLUS reagent kit was used. All RNA samples were assessed with spectrophotometric measurements (BioSPEC-Mini, Shimadzu) and quality RNA Integrity Number (RIN) with Agilent Bioanalyzer, (ii) for target preparation, 100 ng of total RNA was reverse transcribed to generate cDNA, and later used as template to generate biotin-labeled amplified RNA (aRNA). aRNA was then fragmented and hybridized to Affymetrix Human U133 Plus 2.0 Arrays, for 16 hours (45 °C; 60 rpm rotation). Affymetrix Human U133 Plus 2.0 Arrays (HG-U133 Plus 2.0) comprised of 1,300,000 unique oligonucleotides, encompassing 47, 000 transcripts and variants, representing approximately 39, 000 of well characterized human genes, providing a complete coverage to the human genes. All samples were processed with similar reagent kit, were washed, stained and scanned with Affymetrix 300 7 G scanner. Scanned images were assessed for hybridization efficiency.

Immunohistochemistry. (i)
(iii) Prior to microarray analysis, quality control checks were also carried out on the array. The 3′/5′ ratio of housekeeping genes, presences of spike controls, background value, raw Q noise, scaling factor, and percent of genes presence were evaluated. Signal intensity ratio of 3′/5′ probe sets is important in establishing good cDNA synthesis, integrity of starting RNA and hybridization properties. In order to obtain reliable and accurate data comparison, the background intensities were measured and were made sure to be consistent and closed to each other. Additionally, the raw Q noise was also measured by performing pixel to pixel variations in background intensities. This is to make sure that the variations in the digitized signal observed by the scanners as it samples the probe array's surface is consistent and close to each other. The measure of brightness of the array which may vary between arrays to array was also normalized to a standard level. This normalization will be important in comparing array to array data. To assess the hybridization, washing and staining steps quality, bacteria spike controls were performed. These are probe sets hybridized by pre-labeled bacterial spike controls (BioC, BioC, BioD, and Cre in staggered concentration). Finally, Poly A control was also carried out to identify problems with target preparation. These are a set of poly-adenylated RNA spikes of Lys, Phe, Thr, and Dap, prepared in staggered concentration. These are prepared together with the samples throughout cDNA synthesis onwards. All RNA samples have RIN 10 and good quality absorbance (OD) ratio for 260/280 m, (>1.8) and 260/230 (>2.0) (Supp. Table 1, Supp Fig. 5), suggesting high purity of RNA, with minimal degradation. Furthermore, the 3′/5′ ratio of housekeeping gene GAPDH, approximated 1 onwards (Supp. Table 2). The array was also quality checked with consistent background intensity ranging from 30.781324 to 40.71225 (Supp. Table 2). The normalized scaling factors was done to standardize the measure of brightness of the array, and the result of the array falls within the acceptable range of below 3-fold (Supp. Table 2). In the bacteria spike controls, all samples also showed correct signal intensities of the control including BioB, suggesting good hybridization (Supp. Fig. 6). To rule out any problems with the target preparation, PolyA controls were also carried out, and all samples have shown consistent and staggered concentration of the poly A controls (Supp. Fig. 6).

Analysis of array intensities. For the changes in differential genes expression level, Affymetrix
Transcriptome Console was used to analyze the level of expression for each individual gene. One-way between subject ANOVA, was used to assess the statistical significance (p < 0.05).

R e v -G G G T C AG G G G T G G T TAT T G C ) , G FA P ( F w d -C AG AT T C G AG G G G G C A A A AG C , R e v-AGGCTCACTC C CTGTCAAGC), and C x cl 6 (Fwd-TGC GT TGCACT TGT T TAC GC, R e v -C T T C C C G T T C T T C A G G G A G G ) N M E S 1 ( F w d -G C C C A C C A G G C G AT C A ATA C ,
Rev-ACACAGCGAAAGATGAGGCT) were detected with the fluorescent double-stranded DNA binding dye, SYBR Green. qRT-PCR amplification was performed in triplicates for each sample and the results were replicated in four independent experiments. Gel electrophoresis and melting curve analyses were performed to validate PCR product sizes. The expression level of each gene was normalized against β-actin using the comparative CT method 70 .
Subnetwork analysis. Candidate subnetworks were generated by inducing connected components on known biological pathways with highly expressed genes in each phenotype. Two scores are computed for each subnetwork; these scores denote the level of expression of the subnetwork in majority of the samples in each phenotype. Finally, the difference of the two scores is tested for statistical significance. The theoretical t-distribution is used as the null distribution for estimating the statistical significance of subnetworks scored in PFSNet. An additional criteria set was that any subnetwork tested for statistical significance, has to be highly expressed in all the samples in the corresponding group (control/IL1β). Changed in subnetworks were analyzed with pathway information from the PathwayAPI database 34 which contains the aggregation of human pathways from KEGG 71 and Ingenuity. Expression data is preprocessed using the espresso function in R affy package 336 . Pathways maps were generated based on KEGG pathways database tool 3 . Each subnetwork was ranked by their effect size. Effect size is a quantitative measure of difference between two groups. In our study, the effect size for each subnetwork was computed as the standardized mean of paired difference "(mean of paired differences)/(standard deviation of paired differences)" between PFSNet scores, corresponding to the control and reactive astrocyte groups 72,73 . As the sample size is small, multiple-testing correction was not performed since the p-value of the same gene fluctuates in a small range, with a large portion in the insignificant part even though the gene is differentially expressed by construction. Multiple-testing correction approaches would simply shift the null-hypothesis rejection threshold left-wards and thus would be insensitive against the wide fluctuation range in the p-value of this gene. We have also discuss similar approach previously 62,63 .

Subnetworks in pathways
Effect size Genes  Imaging. The samples were imaged using ZEISS LSM 810 confocal microscopy instrument. Image J (Fiji) was used for image analysis and quantification. Tile scan and stitching were performed with ZEN blue software to image the entire spinal cord slice. Quantification for spinal cord GFAP intensity was done over 5 random Region of Interest (ROIs), covering the grey and white matter. The ROI is a 200 × 200 μm 2 square visualized with (Plan-APOCHROMAT Zeiss NA = 0.45) 10X objective lens. 3 slices with 400 um between each of them were used for each rat. For in vitro cell length and processes quantification, 20X Plan-APOCHROMAT (NA = 0.8) objective lens was used. Cells that displayed processes were selected for quantification since not all the cells displayed a morphology that included processes. The main primary process of the cell was defined as any extension protruding from the cell body. The length of the primary process was defined as the distance from the center of the nuclei to the tip of the extended process. Only the main primary processes were used for quantification since very few cells displayed branching processes after exposure to IL-1B. For morphological quantification, the grayscale images were converted to binary type using ImageJ and their threshold was adjusted to clearly demarcate the cell boundary. The wand tool was used to outline the cellular perimeter in order to measure the cell surface area. Measuring was done using the line tool in Image J. Imaging settings were fixed for all groups within each experiments.
Statistical analysis. All results were expressed as mean ± SEM, unless stated. Experiments were repeated independently for four times (n = 4) unless stated. Statistical significance was evaluated with unpaired student t-test for significance (*p < 0.05; **p < 0.01).