Hydrogel scaffolds promote neural gene expression and structural reorganization in human astrocyte cultures

Cell culture

The Normal Human Astrocyte (NHA; Lonza, CC-2565) cell line was implemented in this study because its neuroglial lineage and species origin present a practical in vitro system for studying the effects of hydrogel biomaterials on the human brain. Furthermore, because NHA are a commercially available cell line, they can be used in follow-up and confirmatory experiments by our research team as well as those of other investigators. NHA from two different biological donors (lots #000080982, #000022529; referred to hereafter as Donor A and Donor B) were cultured according to the manufacturer’s specifications. All cytokines, growth factors, and supplements from the SingleQuots™ kit (Lonza; CC-4123) were added to Astrocyte Basal Medium (Lonza; CC-3187). We omitted the manufacturer-recommended gentamicin from cell media because the aseptic techniques used in our laboratory enable NHA culture in an antibiotic-free environment. Frozen ampules of cells were thawed and plated into four T-25 flasks (passage 0) and incubated at 37°C, 5% CO₂, 275 mOsm. Media were replenished within 24 h of thawing cells, and every 48 h that followed. NHA were subcultured by partial digestion with ReagentPack™ subculturing reagents (Lonza; CC-5034) when cultures reached 80% confluence, five days after plating (passage 1).

Matrix preparation

We selected a self-assembling 16-mer peptide hydrogel matrix (PuraMatrix™; Ac-RADARADARADARADA-CONH₂; 3DMatrix Medical Technology, formerly BD Biosciences) for our experiments because it is chemically defined, commercially available, and cells are readily dissociated from the matrix for RNA isolation. The limited biological activity of PuraMatrix™ facilitates the evaluation of effects that occur primarily in response to the microenvironment structure. In particular, experiments have shown that PuraMatrix™ is not immunogenic, cytotoxic, pyrogenic, or hemolytic, and does not bind to cells via integrin receptors (Holmes et al., 2000; Zhang, Ellis-behnke & Zhao, 2005). Moreover, physiological conditions are retained in PuraMatrix™ because the osmolarity of the culture medium is not altered by addition of the hydrogel (matrix = 275 ± 2 mOsm; monolayer = 276 ± 4 mOsm; mean + S.D., n = 3).

PuraMatrix™ (1%) was vortexed for 30 s, diluted in sterile water to 0.35%, and vortexed again for 30 s. The matrix was centrifuged at 210 G for 5 min to remove bubbles. A total of 125 µl of peptide hydrogel per cm² was added to each culture vessel prior to the induction of gelation. Self-assembly of the microenvironment was initiated by adding an equivalent volume of cell culture media and equilibrating for one hour. In order to bring the environment to a neutral pH, the media were replenished twice over the following hour and allowed to incubate overnight prior to the addition of cells.

Experimental design

Two biological replicates (Donor A, B) were cultured simultaneously under two experimental conditions (TCPS and PuraMatrix™) for RNA isolation and sequencing (Fig. 1). The parallel culture and subsequent RNA isolation from donor cultures maintained in both microenvironments was designed to minimize the false discovery of differences due to technical variation. T-25 flasks for each NHA donor and each experimental condition were seeded with 5,000 cells per cm². As recommended by the vendor, the media were replenished on the 3rd day after plating, and RNA was isolated when monolayer cultures reached approximate confluence (day 5).

Live cell imaging with phase contrast microscopy

Prior to RNA isolation, phase contrast images of all four live samples were captured with Metavue image capture software (Molecular Devices), which controlled a Coolsnap HQ CCD camera (Photometrics) attached to an inverted Nikon TE-2000 microscope (20X Magnification).

RNA isolation

RNA isolation was undertaken according to the instructions in the TRIzol^® Plus RNA Purification Kit (Ambion). DNA was removed using the DNA-free™ kit (Ambion) according to the instructions provided by the manufacturer. RNA quality was assessed with an Agilent 2100 Bioanalyzer.

Library preparation and sequencing

RNA samples with RIN Values greater than 8.8 were submitted to the BioMicro Center at the Massachusetts Institute of Technology for library preparation and sequencing (Fig. S1). A 500 ng aliquot of total RNA (500 ng) from each sample was poly-A purified and converted to cDNA using the manufacturer’s instructions for the Illumina TruSeq RNA Sample Preparation Kit. Samples were fragmented on the SPRI-works system using BioMicro Center adapters, barcoded for multiplexing, enriched using BioMicro Center PCR primers, and assessed for fragment size and distribution on an Agilent 2100 Bioanalyzer. Samples were multiplexed and sequenced on an Illumina HiSeq 2000 (Illumina) in accordance with the protocol for 50 base pair (bp) paired end (PE) reads. Phred scores and nucleotide composition were assessed with FastQC to evaluate base call accuracy (Figs. S2 and S3). All four Illumina libraries were multiplexed into one lane (to minimize lane variation) and sequenced on two separate runs of the HiSeq 2000. RNA sequencing (RNA-seq) data collected from all four samples on these two separate occasions are considered technical replicates and labelled as Lane 1 and Lane 2 in the manuscript and data files.

Statistical and ontological analyses

RNA sequencing data analysis was undertaken with freely available online tools that included open source software and the statistical programming language, R. Illumina reads were aligned to the UCSC human genome (hg19) using TopHat software (Version 2.0.8b) which is made freely available on GenePattern by the Broad Institute (Langmead et al., 2009). The countOverlaps function in the R package Genomic Ranges (version 1.16.4) was used to assemble transcripts and estimate abundance (Lawrence et al., 2013). The abundance estimates produced by Genomic Ranges are not normalized and therefore referred to as raw data throughout this manuscript. Cuffdiff software (Version 2.1.1), made freely available on GenePattern by the Broad Institute, was used to normalize transcript abundance by the effective library size and transcript length. Prior to differential expression analysis, the raw read counts for technical replicates were summed to evaluate only the differences between true biological replicates (n = 2). The R package DeSeq (version 1.16) was used to evaluate differential expression because of the package’s ability to identify differentially expressed genes (DEGs) even among samples with low numbers of biological replicates (Anders & Huber, 2010). The first generation ‘DeSeq’ algorithm was chosen instead of the newer ‘DeSeq2’ package because ‘DeSeq’ estimates differential expression using a more conservative algorithm that is less likely to reveal false differentially expressed genes (Love, Huber & Anders, 2014). Coding transcript sequences (CTSs) for which DeSeq reported a Benjamini–Hochberg adjusted p-value of ≤0.1 and a log₂ fold change ≥2 met the criteria for differential gene expression. Cufflinks (version 2.2.1) was used to normalize raw read counts by fragments per kilobase of exon per million reads mapped (FPKM) to compensate for bias introduced by transcript length variance (Trapnell et al., 2012).

Functional significance was imparted to highly and differentially expressed genes using DAVID, the Database for Annotation, Visualization, and Integrated Discovery, with medium stringency settings [29,30]. The Gene Ontology Consortium (GO) terms, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, as well as keywords and sequence features from Uniprot were analyzed in detail for genes assigned to clusters with an enrichment score > 1.3 (Langmead et al., 2009; Sur et al., 2012; Lawrence et al., 2013). Protein networks were identified with STRING using the high stringency filter and a maximum of 10 interacting nodes (Szklarczyk et al., 2014). The role of DEGs in nervous system function, if any, was established through manual curation by PubMed query of each DEG using logical “AND” with terms “brain” or “neural,” and by searching the text in GeneCards (Safran et al., 2010).

Immunocytochemistry

Chambered slides with monolayer- and matrix-cultured NHA were rinsed twice in 1 ml/chamber of phosphate buffered saline (PBS) before fixing in an ice cold mixture of methanol and acetone (1:1) for 15 min (1 ml/chamber). Fixed cells were rinsed twice in ice cold PBS, then samples were incubated for one hour at ambient temperature (∼26°C) with 1 ml/chamber of blocking solution of PBS comprising 10% normal goat serum, 1% bovine serum albumin (BSA), 0.3 M glycine, and 0.1% TWEEN^® 20. Primary and secondary antibodies were diluted in PBS containing 1% BSA and 0.1% TWEEN^® 20 (PBS-ab). Rabbit anti-TUBB3 antibody (Abcam; catalog # ab52623, lot # GR14955-15) was diluted at 1:300 in PBS-ab and 250 µl of antibody solution was incubated with samples overnight at 4°C in a sealed, humidified chamber. Primary antibody was omitted from negative control samples to assess nonspecific binding of the secondary antibody. Cells were rinsed 3 × 5 min in 1 ml/chamber PBS before incubating with 250 µl of Alexa Fluor^® 488-conjugated goat anti-rabbit IgG (Abcam; catalogue # 150081, lot # GR282725-1) diluted at 1:500 in PBS-ab for one hour at ambient temperature (∼26°C) in the dark. Cell nuclei were counterstained with the nucleic acid stain Hoescht 33342 (0.2 µg/ ml) diluted in 18 MΩ water. Cells were rinsed 3 × 5 min in PBS before mounting in glycerol and PBS (1:1). Confocal images were acquired with a 10X/0.3 objective mounted on an inverted LSM 700 microscope (Zeiss). Hoescht 33342 and Alexa Fluor^® 488 were excited with laser excitation at λ = 405 and λ = 488, respectively. Main beam splitter was MBS 405/488/555/639 with the SP 555 filter used to collect emission from both fluorophores. Images were captured with settings which excluded fluorescence from control samples incubated with only secondary antibody. Immunolabelling and image capture were repeated in triplicate on the second and third passages of NHA (lot #0000402839). A blind, qualitative image assessment was undertaken by three individuals asked to combine images into groups, if applicable, and describe the staining patterns.

Figure preparation

Figures were prepared with Wordle (Feinberg, 2014), Excel (Microsoft, 2010), and the ‘ggplot2’ (Wickham, 2009) and ‘pheatmap’ (Kolde, 2015) data visualization packages for the R statistical programming language. Additional image assembly, labelling, and alignment were completed with Photoshop (Adobe CS6). The ZenLightEdition software package (Zeiss 2009) was used to optimize the relative contribution of the fluorescence signal from the emission channels to the merged confocal image, and to export confocal images as JPEG files.

Responsible conduct and reproducibility

Experimental procedures were developed in accordance with established guidelines for preclinical research (Landis et al., 2012; The National Institute of Health, 2016). NHA were produced by the vendor (Lonza) in compliance with national ethics standards (document available upon request from vendor). Laboratory culture procedures adhered to vendor specifications; cells were used within the maximum 10 population doublings (three passages). RNA-seq results (see below) further confirmed the human origin of the cell line.

The use of the commercial (human) cell lines was approved by the NMSU Institutional Biosafety Committee, approval # 1401SE2F0103, “Gene Expression in the Nervous System.” The protocol was exempt from review by the Institutional Review Board because the commercial cell lines were developed before 2015 and were de-identified by the vendor, Lonza, who retains a signed record of informed consent from human donors. The NIH Extramural Institutional Certification form for human cell lines created before January 25, 2015 was signed by the institutional offices and the principal investigator (EE Serrano) and submitted to NIH with the GEO data files. In accordance with data sharing policies set forth by the NIH, both the raw and processed data are available for download from GEO using the accession number GSE81995.

Official gene symbols were used in accordance with the standards set forth by the HUGO Gene Nomenclature Committee (Gray et al., 2015). Analysis was undertaken with two biological replicates because the vendor provides primary cell lines from two different human donors. A balanced and blocked design of two biological and two technical replicates was implemented to minimize the impact of confounding variables (Auer & Doerge, 2010). For a complete description of the statistical test that was used to evaluate differential expression, see Anders & Huber (2010). Sequence data were collected with a single blind protocol at the MIT BioMicro center, without prior knowledge of the nature of the biological samples. Researchers at NMSU assigned groups and assessed the outcomes using a single blind analysis.

Rabbit anti-TUBB3 monoclonal antibody (Abcam; catalog # ab52623, lot # GR149555-15, RRID:AB_869991) was used for immunocytochemistry with human astrocytes, as validated by the vendor with western blot analysis and confirmed in several references listed on the vendor’s website.

Live cell imaging with phase contrast microscopy

Phase contrast images of cells on TCPS and with peptide hydrogel demonstrate structural differences between NHA grown in TCPS or hydrogel microenvironments (Fig. 2). NHA cultured on TCPS surfaces (Figs. 2A and 2C) appear uniformly flat and phase dark, as is typical of astrocyte cultures, in comparison to NHA cultured with PuraMatrix™ (Figs. 2B and 2D). In the hydrogel environment, NHA reorganize and form aggregate structures that resemble neuron-glia cocultures. Clusters of phase-bright cell regions that are characteristic of neuronal somata are visible interspersed in the matrix among phase dark astrocytes (Serrano & Schimke, 1990; Abd-el-basset, 2013).

Evaluation of sequence quality and read alignments

Sequence quality assessment of PE RNA-seq reads from NHA revealed that, for all samples and read positions, the average Phred scores for nucleotides met the criteria for 99% base call accuracy (Phred score > 20; Fig. S2). The percent nucleotide composition values fluctuate during the initial reads (≤13), then converge to a more constant value (25%) as the reads lengthen (Fig. S3). On average, over a hundred million reads reads per replicate sample aligned to the reference human genome (hg19) for a total of 856,790,008 RNA-seq read alignments (Table 1). Approximately 80% of the reads aligned uniquely. Both paired ends aligned to the reference genome (hg19) for 45% of the PE reads, while 2% of alignments were discordant and did not meet the expected distance and/ or orientation constraints established by TopHat.

Normalization and establishment of threshold for expression

Box plots of log₁₀-transformed raw RNA-seq alignment reads and log₁₀-transformed FPKM-normalized coding sequences demonstrated consistency between the technical replicates, with slightly more variation between the biological replicates (Fig. 3). The interquartile ranges (25th–75th percentiles) for the raw read data were larger, and count values were greater, as compared to FPKM normalized data. In Fig. 4, the effects of normalization are portrayed in heat maps of Euclidean distances grouped by hierarchical clustering. Inspection of clusters and their internodal distances illustrates that, in both normalized and raw data, the most similar expression patterns are observed among the technical replicates. Cluster analysis of the raw data (Fig. 4A) linked the monolayer sample from donor B closest to the matrix sample from donor A. In contrast, the normalized data (Fig. 4B) illustrate the closer relation of expression between monolayer and matrix samples from the same donor. In Fig. 5, the effects of normalization are reiterated in principal component analyses. Inspection of groups in normalized and raw data revealed that the most similar expression patterns are observed among the technical replicates. Principal component analysis of the raw data (Fig. 5A) grouped the monolayer sample from donor B closest to the matrix sample from donor A. In contrast, the normalized data (Fig. 5B) illustrate the closer relation of monolayer and matrix samples from the same donor.

In order to evaluate genes with estimated mRNA copy numbers greater than one per cell, we set a detectability threshold of FPKM ≥ 1 (Hebenstreit et al., 2011; Nagaraj et al., 2011) and only CTS that met this criterion in at least one sample were eligible for inclusion in our detailed analysis. These criteria excluded 24 CTSs with a status of “NA” (File S1), 2,183 CTSs with all eight FPKM values of 0 or NA (File S2), and 6,447 CTSs with FPKM values from 0.1–1 (File S3). Moreover, 1,806 genes representing noncoding RNA sequences were removed from subsequent analyses (File S4). Application of these filter criteria identified 12,822 CTSs for downstream analyses (File S5).

Characterization of the most abundant transcripts

First, FPKM values of coding transcript sequences for cells grown in monolayer and matrix conditions were ranked from highest to lowest. Next, we identified the top 10% (n = 128) of the 12,822 CTSs with the highest average FPKM for cells grown in monolayer and matrix conditions. Approximately 90% (n = 115) of the most highly transcribed CTSs were common to cells grown in monolayer and matrix conditions, while 26 were not common to both conditions. The most highly transcribed common CTSs (n = 115) were used to construct a heat map, which depicts the FPKM for the CTSs (Fig. 6). DAVID analysis clustered the most expressed genes common to NHA grown in TCPS and hydrogel environments (n = 115) into 22 significant a clusters (enrichment scores > 1.3; File S6). The GO terms corresponding to the most highly transcribed CTSs were arranged into a Wordle image that represents higher frequency with larger font size (Fig. 7; Feinberg, 2014).

A total of 24 of the 26 CTSs that were not common to both conditions were found in the top 15% of the most transcribed genes in the other condition. Two genes, ACTA2 and HTRA1, ranked much lower in the other condition. When CTSs were ranked by highest average FPKM, the CTS for ACTA2 was ranked #61 for cells grown on TCPS; this ranking fell to #476 (37%) of the highest average FPKM for cells grown in PuraMatrix™. The ACTA2 gene corresponds to the alpha actin 2 protein, which is involved in the formation of tight junctions in the blood brain barrier (Safran et al., 2010). ACTA2 was not considered to be differentially expressed between the two conditions using our DEG criteria; however, the log₂FC value was −2.5. For cells grown in PuraMatrix™, the CTS for HTRA1 was ranked #122 of the most transcribed 12,822 genes; this ranking fell to #739 (58%) of the highest average FPKM for cells grown on TCPS. HTRA1 is a serine protease that plays a crucial role in neuronal maturation, possibly through downregulation of the TGF-β signalling pathway (Launay et al., 2008). HTRA1 met the differential expression criteria with a log₂FC = 2.2.

Identification of CNS biomarkers in monolayer- and matrix-cultured NHA transcriptomes

We compared the transcriptome of NHA in monolayer and matrix conditions with established profiles for distinct CNS cell types (Fig. 8). As per Cahoy et al. (2008), a suite of marker genes specify several classes of murine neural cells, including astrocytes, oligodendrocytes, and neurons, (Cahoy et al., 2008). Of the genes that Cahoy et al. (2008) determined are significantly upregulated in astrocytes, 60 % and 61 % met our detectable threshold criteria in monolayer and matrix conditions, respectively (Fig. 8A) (Cahoy et al., 2008). Similarly, 56% and 58% of the genes that are designated as significantly upregulated in oligodendrocytes by Cahoy et al. (2008) met our detection criteria in in monolayer and matrix conditions, respectively (Fig. 8A) (Cahoy et al., 2008). For the genes found upregulated in neurons by Cahoy et al. (2008), 47% and 44% were detectable in monolayer and matrix-cultured NHA cells, respectively (Fig. 8A) (Cahoy et al., 2008).

We further probed for cellular heterogeneity in NHA cultures using marker genes for distinct CNS cell subtypes. Cahoy et al. (2008) developed subset profiles (80 genes each) of the most upregulated genes in specific classes of astrocytes (in vitro, in vivo, developing, and mature) and oligodendrocytes (precursors, myelinating) (Cahoy et al., 2008). When we undertook a more detailed analysis of the NHA transcriptome in monolayer and matrix conditions using these gene lists, we observed that markers for all CNS cell subtypes were present in similar amounts for both monolayer- and matrix- cultured NHA. However, the relative amount of NHA expressed genes varied between subtypes. For example the subtype in vitro astrocytes comprised the highest number of biomarker genes (monolayer, 78%; matrix, 78%; Fig. 8B), as compared to myelinating oligodendrocytes, where a fewer number of biomarkers were detectable (monolayer, 34%; matrix, 36%; Fig. 8B).

Ontological analysis of differentially expressed genes with DAVID

A total of 43 CTSs met our criteria for differentially expressed genes (DEGs; Table 2). Of these DEGs, 16% were upregulated in monolayer-cultured NHA, and 84% were upregulated in cells cultured with PuraMatrix™ (Fig. 9). DAVID analysis did not group the DEGs from TCPS-cultured cells into any clusters (Huang, Sherman & Lempicki, 2009a; Huang, Sherman & Lempicki, 2009b). In contrast, DAVID analysis assigned the DEGs upregulated in hydrogel-cultured cells into 9 significant clusters (enrichment score > 1.3) which we ranked 1–9 with the highest enrichment score corresponding to cluster 1 (Table 3 and File S7) (Huang, Sherman & Lempicki, 2009a; Huang, Sherman & Lempicki, 2009b).

Predicted protein interactions for DEGs

We used the STRING database to conduct a deeper exploration of the network of protein interactions for the DEGs. We configured STRING for ten additional interacting nodes beyond the 7 and 36 DEGs used as input for the monolayer-cultured and matrix-cultured cells, respectively (Table 2 and Fig. 10). For cells cultured on TCPS, the top 10 interacting proteins predicted by STRING were EGFR, NPR3, UBC, ERBB4, NEB, MYH11, TPM2, NPR2, TPM1, and TPM3 (Fig. 10A). The proteins encoded by the genes upregulated in cells cultured with PuraMatrix™ were predicted to interact with ITGB1, EFNB1, NOTCH1, MAPK8, EFNB2, USH1C, HBD, CD44, EGFR, and TP53 (Fig. 10B).

The predicted functional significance of the protein interactors was explored with GeneCards (Safran et al., 2010). The CNS emerged as a common theme, and was set as a goal for further exploration. Seven of the top ten interacting proteins for matrix-upregulated genes had an established role in the nervous system that was mentioned in GeneCards (ITGB1, migration and survival of primary oligodendrocytes; EFNB1, mitigates commissural axons and growth cones; NOTCH, differentiation of Bergmann glia, represses neuronal differentiation; TP53, implicated in NOTCH signaling cascade, response to oxidative stress; MAPK8, controls neurite elongation in cortical neurons; EFNB2, constrains longitudinally projecting axons, angiogenesis, neuronal development; USH1C, part of the network that mediates mechanotransduction in cochlear hair cells). In contrast, exploring data on GeneCards uncovered an established role in the nervous system for only one of the top ten interacting proteins for monolayer-upregulated genes (ERBB4, nervous system development) (Safran et al., 2010).

Significance of differentially expressed genes for CNS function

The relevance of the DEGs in the nervous system was probed through pooled information sources (DAVID, GeneCards) as well as primary literature queries (PubMed) of individual genes along with the terms “brain” and “neural.” One of the seven genes upregulated in monolayer-cultured cells had a previously established neural function. Of the 36 genes upregulated in hydrogel-cultured cells, 34 were assigned to CNS function either through manual curation or functional clustering by DAVID (Fig. 11) (Huang, Sherman & Lempicki, 2009a; Huang, Sherman & Lempicki, 2009b; Safran et al., 2010). Of the hydrogel-upregulated genes, ∼42% (n = 15) have an established role in response to neural insult and correspond to processes such as oxidative stress (and/or angiogenesis), ischemic stroke, and traumatic brain injury. Approximately 39% (n = 14) of the genes upregulated in matrix-cultured cells are involved in neuronal differentiation and growth, synapse formation, and neural plasticity.

Moreover, over a quarter (n = 10) of the genes upregulated in matrix-cultured NHA play a role in sensory perception. The manually curated associations are described in more detail below.

Immunocytochemical detection of class III β-tubulin

In order to characterize the morphological differences observed in TCPS- and hydrogel- cultured NHA, class III β-tubulin was labelled using immunocytochemical methods. A pervasive, positive stain for class III β-tubulin was observed throughout the cell body of every cell cultured on TCPS (Fig. 12A). This result is consistent with the upregulated expression of class III β-tubulin in human fetal brain tissue (Safran et al., 2010). In contrast, a more heterogeneous staining pattern for class III β-tubulin was seen in matrix-cultured NHA (Fig. 12B). Some cells appeared unlabelled, and there was a spatially variable pattern of fluorescence in the cells that stained positively. In particular, long protrusions from the cell body of some cells appeared to label more intensely than other regions of the cell. Three different individuals blindly grouped optical sections and maximum confocal projections into two conditions and verified the labelling patterns seen in the two groups.