Retinal astrocytes transcriptome reveals Cyp1b1 regulates the expression of genes involved in cell adhesion and migration

Astrocytes (AC) are the most abundant cells in the central nervous system. In the retina, astrocytes play important roles in the development and integrity of the retinal neurovasculature. Astrocytes dysfunction contributes to pathogenesis of a variety of neurovascular diseases including diabetic retinopathy. Recent studies have demonstrated the expression of Cyp1b1 in the neurovascular cells of the central nervous system including AC. We recently showed retinal AC constitutively express Cyp1b1, and global Cyp1b1-deficiency (Cyp1b1-/-) attenuates retinal ischemia-mediated neovascularization in vivo and the pro-angiogenic activity of retinal vascular cells in vitro. We also demonstrated that Cyp1b1 expression is a key regulator of retinal AC function. However, the underlying mechanisms involved need further investigation. Here we determined changes in the transcriptome profiles of Cyp1b1+/+ and Cyp1b1-/- retinal AC by RNA sequencing. We identified 585 differentially expressed genes, whose pathway enrichment analysis revealed the most significant pathways impacted in Cyp1b1-/- AC. These genes included those of axon guidance, extracellular matrix proteins and their receptors, cancer, cell adhesion molecules, TGF-β signaling, and the focal adhesion modulation. The expression of a selected set of differentially expressed genes was confirmed by RT-qPCR analysis. To our knowledge, this is the first report of RNAseq investigation of the retinal AC transcriptome and the molecular pathways impacted by Cyp1b1 expression. These results demonstrated an important role for Cyp1b1 expression in the regulation of various retinal AC functions, which are important in neurovascular development and integrity.


Introduction
The cytochrome P450 superfamily consists of many heme-containing monooxygenases. They are best known for their roles in drug metabolism. CYP1B1 is involved in many processes in the body, such as assisting with reactions that break down drugs and produce certain fats (lipids). It is expressed in both adult and fetal human extrahepatic tissues, including most of the parenchymal and stromal tissues from brain, kidney, prostate, breast, cervix, uterus, ovary, lymph nodes [1], and ocular tissues [2,3]. Mutations in this enzyme are a risk factor for the development of primary congenital glaucoma in humans [4]. However, the underlying cellular and molecular mechanisms are not fully revealed. We previously showed expression of Cyp1b1 is essential for ischemia-mediated retinal neovascularization as occurs in retinopathy of prematurity, and the proangiogenic function of retinal vascular cells in culture [5][6][7]. However, how Cyp1b1 expression impacts these processes remained largely unknown. We showed Cyp1b1 is constitutively expressed in vascular endothelial cells and perivascular supporting cells from vascular beds of many organs including retina [5,7]. Recently, we also demonstrated that Cyp1b1 is expressed in retinal astrocytes (AC), and Cyp1b1-/-retinal AC are more proliferative and migratory [8]. These cells produced increased amounts of fibronectin and its receptors αvβ3 and α5β1 integrins. However, production of inflammatory mediators such as BMP-7 and MCP-1 were decreased in Cyp1b1-/-AC. In addition, we observed a significant increase in CD38 expression when Cyp1b1-/-AC were challenged with H 2 O 2 compared with Cyp1b1+/+ cells. Cyp1b1-/-AC also showed enhanced connexin 43 phosphorylation compared with Cyp1b1+/+ cells [8]. Thus, Cyp1b1-deficiency in AC was associated with increased resistance towards oxidative stress.
Astrocytes are the major cell type in the optic nerve head and are vital to the development and maintenance of the retinal astrocytic network and angiogenesis [9,10]. Under pathological conditions, AC become reactive and contribute to various ocular pathologies including glaucoma and diabetic retinopathy [11]. However, there still much more to delineate regarding Cyp1b1 expression and function in AC. The few studies available to date have confirmed Cyp1b1 expression in AC and neurons, and its upregulation in a variety of gliomas [1,12]. To our knowledge, we were first to report the impact of Cyp1b1 expression on retinal AC function. However, the intracellular pathways that mediate these activities of Cyp1b1 in AC remain elusive.
RNAseq is one transcriptomic approach where the total complement of RNAs from a given sample is isolated and sequenced using high-throughput technologies [13]. RNAseq technology has the potential to provide very useful, detailed information on the intracellular pathways impacted by Cyp1b1 expression, and identify the networks of genes involved. The purpose of the current study was to utilize this powerful technique to delineate the detailed molecular mechanisms of Cyp1b1 action in retinal AC by determining the changes in patterns of gene expression networks impacted by Cyp1b1 expression. The identification of genes whose expression is affected by the presence or absence of Cyp1b1 will provide additional clues to the intracellular mechanisms of Cyp1b1 action and function in retinal AC. This knowledge should lead to the discovery of new targets for modulation of Cyp1b1 activity and their potential therapeutic use. Expectation Maximization) [15]. Gene expression was normalized by the method of trimmed mean of M-values (TMM), where the product of this factor and the library sizes defines the effective library size [16]. For each sample contrast, simple boxplots per sample expression distributions were constructed before and after TMM normalization (Fig 1A). Overall TMM values were similar in each sample, indicated by the uniform distributions.
Analysis of differentially expressed genes was performed with a glm using the edgeR package [17]. In order to decide which genes are differentially expressed (DEG), the adjusted p-value-not the raw p-value-was defined to be 0.05. To control the false discovery rate (FDR), a Benjamini-Hochberg correction was applied [18]. Fig 1B represents a Volcano plot showing DEG as red and blue dots denoting up-and down-regulated expression, respectively, at an adjusted p-value (FDR) significance threshold of 0.05. The gray dots reflect those genes with no evidence of statistically significant changes in expression. The two solid gray lines denote the boundary of a two-fold change. We also conducted a hierarchical clustering analysis of DEG from all samples with Ward's method of Euclidean distances [19], and created a heatmap with the heatmap function from Heatmapper: webenabled heat mapping for all [20]. The results indicated that gene expression was similar in each group (Fig 1C).
Pathway enrichment analysis, conducted using KEGG (Kyoto Encyclopedia of Genes and Genomes) [28][29][30] as a mapping database, 57 pathways were identified with significance level of 0.05. To cut down on the list the P< 0.002 was applied. The pathways with multiple points of commonality and overlap among the Axon guidance, extracellular matrix-receptor interactions, pathways in cancer, cell adhesion molecules, TGF-β signaling pathways, and the focal adhesion were identified as most significantly enriched pathways ( Table 2).

Discussion
In this study we used RNAseq analysis to determine the global transcriptome profile of Cyp1b1 +/+ and Cyp1b1-/-AC from mouse retina. We identified 585 DEG, whose pathway analysis revealed the most significant biological functions. These included the Axon guidance, extracellular matrix (ECM)-receptor interactions, pathways in cancer, cell adhesion molecules, TGF-β signaling, and the focal adhesion regulation. We also found that some of the top downregulated genes were involved in biological and molecular processes including actin-mediated cell contraction, cardiac muscle contraction, cell proliferation, apoptosis and metabolic processes.
Our findings here were consistent with the results of our previous studies showing increased proliferation and migration of Cyp1b1-/-AC. Activation of AC proliferation and migration is important in repair of injuries in the central nervous system (CNS) and scar formation [31,32]. AC migration is regulated by various factors, among which transforming growth factor-β (TGF-β) plays an important role [33]. In AC, TGF-β suppresses cell proliferation by inducing p15 INK4B expression in a Smad3-dependent manner [34]. This is consistent with our findings that showed the downregulation of TGF-β in Cyp1b1-/-AC. The upregulation of Smad genes in Cyp1b1 -/-AC was associated with the enhanced pathways in axon guidance and cell proliferation.
Our results showed changes in Cxcl12, (also known as SDF-1). Cxcl12 is one of the most studied chemokines that induce cell proliferation and migration by binding to its receptor. Under normal conditions, Cxcl12 expression in the CNS is relatively low. However, its expression is upregulated when the CNS is affected by trauma, ischemia, inflammation or infection [35]. Enhanced Cxcl12 expression promotes proliferation of radial glia like cells after traumatic brain injury in rats [36]. Others have shown its therapeutic value by promoting autophagy and migration via PI3K-AKT-mTOR pathway [37].
We recently showed that Cyp1b1 deficiency affects retinal AC ECM production and expression of integrin receptors [8]. We also showed upregulation of cadherins, laminin, and tenascin. These molecules are well known for their roles in cell adhesion [38], and were associated with changes in adhesion observed in Cyp1b1-/-retinal AC [8]. Cadherin-4 (also known as Rcadherin) is involved in retinal angiogenesis during development. Dorrell et al. [39] used antibodies or peptides to neutralize R-cadherin, which prevented the normal formation of the retinal vascular network in newborn mice. They also showed that R-cadherin plays a crucial role in the endothelial-astrocyte interactions [39,40]. Thus, Cyp1b1 expression may impact AC interactions with EC.
Oxidative stress is implicated in many neurodegenerative diseases. Cytochrome P450 activities are generally involved in ROS production due to their involvement in the metabolism of steroids, fat-soluble vitamins, fatty acids, eicosanoids, drugs, carcinogens, and other xenobiotic chemicals [41][42][43]. CYP enzymes can generate superoxide and hydrogen peroxide through uncoupling reactions, for more details see [44]. We have demonstrated an important role for Cyp1b1 as a modulator of cellular redox state [45]. Studies utilizing vascular cells derived from Cyp1b1-/-mice showed an increase in oxidative stress in vascular endothelial cells and perivascular supporting cells [5,7,46]. In contrast, Cyp1b1-/-retinal AC did not show an increase in oxidative stress compared to Cyp1b1+/+ cells under basal conditions [8]. This is consistent with minimal changes in expression of genes that affect cellular redox state in the absence of Cyp1b1. However, incubation of these cells with known inducers of oxidative stress could reveal changes in genes that modulate oxidative stress. This notion is supported by our studies demonstrating that Cyp1b1-/-AC elicit a significantly more robust response in expression of CD38 when challenged with H 2 O 2 compared with Cyp1b1+/+ cells [47]. Thus, the elucidation of the mechanisms behind AC resilience to oxidative stress, especially in the absence of Cyp1b1, needs further investigation.
Cytochrome P450 enzymes are involved in metabolism of drugs, and are major source of variability in drug pharmacokinetics and responses. However, in some cases they can also activate compounds consumed in food, converting pro-carcinogens to carcinogens [48]. CYP proteins, involved in steroid or retinoic acid metabolism, could promote or suppress tumors development through hormonal control [49,50]. Genetic variability could play a role if a polymorphism affected a CYP protein involved in such processes [51]. Our study found multiple genes in the cancer pathway that are associates with the progression or suppression of cancer such as laminin, Tgfβ3, Hhip, Arnt2 and Mmp2 [52][53][54]. Mmp2 was upregulated in this pathway and is implicated in cancer cell migration [55]. A previous study using a microarray analysis, demonstrated that the increased expression of Mmp2 is involved in invasiveness of malignant glioma [56,57], an observation that is consistent with our findings. Thus, Cyp1b1 expression has an important role in modulating AC migration by suppressing Mmp2 expression.
A limitation of our studies is the use of cells prepared from mix genders. Our initial in vivo and in vitro vascular cell culture studies did not demonstrate a gender bias in the noted phenotypes with Cyp1b1-deficiency. However, our gene expression studies here showed some of the differentially expressed genes in Cyp1b1 null cells are sex linked: Uty; Y-chromosome, Cysltr1; X-chromosome, and Kdm5d; Y-chromosome. Thus, Cyp1b1-deficiency impact on gene expression, and likely noted phenotypes, could be deferentially impacted by gender. Future studies will further address the gender contributions to various phenotypes noted with Cyp1b1-deficiecny.
In summary, RNAseq technology was used to investigate the transcriptome profiles of retinal AC and how Cyp1b1 expression modulates their cellular functions. A pathway analysis of DEG indicated the most significantly enriched pathways included Axon guidance, ECMreceptor interactions, as well as cancer and other pathways (cell proliferation, focal adhesion, and cell adhesion). Our transcriptomic approach in this investigation, which relied on  RNAseq, was powerful and effective way to allow us to obtain a global view of genes whose expression are impacted by Cyp1b1, likely in a gender dependent manner.

Ethics statement
All animal experiments were performed in accordance to the Association for Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Vision Research and were approved by the Institutional Animal Care and Use Committee of the University of Wisconsin School of Medicine and Public Health (the assurance number A3368-01). Animals were sacrificed according to an approved protocol by CO2 asphyxiation.

Isolation and culture of Cyp1b1-/-retinal AC
Retinal AC were isolated from mouse retina by collecting retinas from one litter of 4-week-old (6 to 7-mix gender) mice using a dissecting microscope, as previously described by us with greater than 98% purity [8,58]. Briefly, retinas (12 to 14) were rinsed with serum-free Dulbecco's Modified Eagle's Medium (DMEM), pooled in a 60 mm dish, minced and digested for 45 min with collagenase Type I (1 mg/ml; Worthington, Lakewood, NJ) in serum-free DMEM at 37˚C. Cells were rinsed in DMEM containing 10% fetal bovine serum (FBS) and centrifuged for 5 min at 400 xg. Digested cells were rinsed again in DMEM containing 10% FBS and filtered through a double layer of sterile 40 μm nylon mesh (Sefar America Inc., Fisher Scientific, Hanover Park, IL). Cells were centrifuged for 5 min at 400 xg and medium was aspirated. Cells were washed twice with DMEM containing 10% FBS, resuspended in 1 ml of DMEM containing 10% FBS in a 1.5 ml microfuge tube and incubated with rat-anti-mouse CD31 (Mec13.3; BD Biosciences) coated with sheep anti-rat magnetic beads, and were gently rocked for 1 h at 4˚C. Using a Dynal magnetic tube holder, cells not bound to magnetic beads were collected and washed in DMEM containing 10% FBS. Cells were plated in growth medium in a single well of a 24 well plate coated with human fibronectin (2 μg/ml in serum-free DMEM; BD Biosciences, Bedford, MA), and incubated at 33˚C with 5% CO2. The cells bound to magnetic beads are generally used to remove retinal EC as we described [58]. Retinal AC were grown in DMEM containing 10% FBS, 2 mM L-glutamine, 2 mM sodium pyruvate, 20 mM HEPES, 1% nonessential amino acids, 100 μg/ml streptomycin, 100 U/ml penicillin, freshly added heparin at 55 U/ml (Sigma, St. Louis, MO), endothelial growth supplement 100 μg/ml (Sigma), and the murine recombinant interferon-γ (R&D, Minneapolis, MN) at 44 U/ml. Cells were maintained at 33˚C with 5% CO2. Cells were progressively passed to larger plates, maintained, and propagated in 1% gelatin-coated 60 mm dishes. For all experiments, cells were used at 70-80% confluence unless stated otherwise.

RNA purification
Total RNA from Cyp1b1+/+ and Cyp1b1-/-retinal AC was purified using RNeasy Mini kit according to manufacturer's protocol with the DNAse treatment step to eliminate traces of genomic DNA (Qiagen, Germantown, MD). The quality and quantity of the total RNA were measured using an Agilent Model 2100 Bioanalyzer, and samples showing a RIN >8 were selected for further analysis. Samples were stored in RNase-free water and kept at -80˚C until further processing.

RNA sequencing
Eight samples of purified RNA (4 from Cyp1b1+/+ retinal AC and 4 from Cyp1b1 -/-retinal AC) were subsequently subjected to a double round of poly-A mRNA purification, fragmented, and primed for cDNA library synthesis using the TruSeq RNA sample preparation kit (RS-122-9004). All procedures were carried out according to the manufacturer's instructions (Illumina, San Diego, CA). Following validation (Agilent 2100 Bioanalyzer, DNA 1000) and normalization, samples were clustered (TruSeq pairedend cluster kit v3-cBot-HS, PE-401-3001) followed by paired-end sequencing (100 bp; TruSeq SBS kit v3-HS 200 cycles, FC-401-3001) on a HiSeq2500. The following quality control statistics were used to evaluate the technical quality of the experiments; (1) combined per cycle base quality, (2) per cycle base frequencies, (3) per cycle average base quality, (4) relative 3k-mer diversity, (5) Phred quality distribution, (6) mean quality distribution, (7) read length distribution and (8) read occurrence distribution using the trimming software skewer [59]. Low-abundance genes with a read count below a threshold of 1.0 in two or more samples were excluded. To compare gene expression between Cyp1b1+/+AC and Cyp1b1-/-AC, samples were normalized by trimmed mean of M-values (TMM) using Edge R (version 2.5 of Bioconductor) software [16]. After the inspection of preliminary data, transcript reads were aligned to the preassembled selected reference genome sequence using STAR (Spliced Transcripts Alignment to a Reference) using the default settings [14]. Transcript abundance were performed by using RSEM (RNASeq by Expectation Maximization) [15]. Subsequently, differential analysis of significant changes in gene expression was performed with a glm using the edgeR package [17] in the different genotype pairs (e.g. Cyp1b1+/+ vs. Cyp1b1 -/-). All sequence data have been deposited in the Gene Expression Omnibus with accession number GSE145103. . Standard curves were generated from known quantities for each of the target gene of linearized plasmid DNA. Ten times dilution series were used for each known target, which were amplified using SYBR-Green qPCR. The linear regression line for ng of DNA was determined from relative fluorescent units (RFU) at a threshold fluorescence value (Ct) to quantify gene targets from cell extracts by comparing the RFU at the Ct to the standard curve, normalized by the simultaneous amplification of RpL13a, a housekeeping gene. All primer sequence used are listed in Table 4.

Biological interaction network and KEGG pathway enrichment analysis
Pathway Enrichment Analysis was performed with WEBGESTALT web analysis software (http://bioinfo.vanderbilt.edu/wg2/) by mapping significantly changed genes to corresponding KEGG enrichment pathways, Gene Ontology enrichment and conducting a hypergeometric statistical test with significant level <0.05 after multiple testing correction [60][61][62]. Significance for pathway level enrichment was defined as having an enrichment score False Discovery Rate (FDR) corrected p-value < 0.05.

Statistical analysis
RNA-seq data were analyzed and gene expressions were normalized by the method of trimmed mean of M-values (TMM) and glm using the edgeR package, respectively. This yielded a total of 13,575 genes. KEGG pathway enrichment analysis cut-off criteria of FDR<0.05 with a | FC|>2 was apply. RT-PCR data were analyzed with student's unpaired t-test (2-tailed) or oneway ANOVA with post-Bonferroni's test for multiple comparisons. P� 0.05 was considered significant. Data are presented as Mean ± SEM from cells with n� 6 (as indicated in figure legends). All data analysis was done in GraphPad Prism or Microsoft Excel.