The effect of Bacopa monnieri on gene expression levels in SH-SY5Y human neuroblastoma cells

Bacopa monnieri is a plant used as a nootropic in Ayurveda, a 5000-year-old system of traditional Indian medicine. Although both animal and clinical studies supported its role as a memory enhancer, the molecular and cellular mechanism underlying Bacopa’s nootropic action are not understood. In this study, we used deep sequencing (RNA-Seq) to identify the transcriptome changes upon Bacopa treatment on SH-SY5Y human neuroblastoma cells. We identified several genes whose expression levels were regulated by Bacopa. Biostatistical analysis of the RNA-Seq data identified biological pathways and molecular functions that were regulated by Bacopa, including regulation of mRNA translation and transmembrane transport, responses to oxidative stress and protein misfolding. Pathway analysis using the Ingenuity platform suggested that Bacopa may protect against brain damage and improve brain development. These newly identified molecular and cellular determinants may contribute to the nootropic action of Bacopa and open up a new direction of investigation into its mechanism of action.


Introduction
Bacopa monnieri (Bacopa), also known as Bacopa monniera, Herpestis monniera, water hyssop or Brahmi, has long been used in Indian traditional medicine (Ayurveda) as a brain tonic to enhance memory performance, learning and concentration [1,2]. These traditional claims have recently been supported by several animal and clinical studies. Animals treated chronically with Bacopa showed better acquisition and improved retention in learning tasks [3][4][5][6][7][8][9][10][11]. Similarly, clinical studies and a meta-analysis of randomized control trials also demonstrated that chronic oral administration of Bacopa (over a period of more than 12 weeks) to healthy subjects resulted in improvements in the subjects' information processing speed, free recall, verbal memory and learning. Bacopa treatment also resulted in a decrease in anxiety, which PLOS ONE | https://doi.org/10.1371/journal.pone.0182984 August 23, 2017 1 / 21 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 improved learning [12][13][14][15][16][17][18][19]. With its low toxicity risks and apparent beneficial effects as a nootropic [20][21][22], Bacopa has been extracted and marketed as a dietary supplement (KeenMind or CDRI08, Soho Flordis International) and its production and imports are tightly inspected by the Food and Drug Administration (FDA) (https://tinyurl.com/ybkmhfes). Notwithstanding its wide availability, the mechanisms of action of Bacopa have yet to be delineated. Several mechanisms have been proposed concerning the nootropic effects of Bacopa. These included alteration of the levels of several neurotransmitters, including serotonin (5-hydroxytryptamine, 5-HT) [6,23,24], acetylcholine [4,13,25] and dopamine [6,20,24,26]. Elevation of the neurotransmitter 5-HT resulted in activation of cAMP response element-binding protein (CREB) and subsequent changes in transcription, protein phosphorylation and histone modification [8,23,27]. Other research groups have suggested that Bacopa regulate pre-and postsynaptic proteins [9] and induce the formation of new dendrites [11,28]. These proposed mechanisms are not mutually exclusive and the discrepancies could be attributed to a selective bias in the targets that were investigated by each research group. In order to better define the molecular and cellular components of Bacopa action, we have applied a deep sequencing technique, RNA-Seq, to identify the changes in the transcriptome upon Bacopa treatment. We have performed both a transcript level and gene level analysis to investigate the changes in gene expression caused by Bacopa. These experiments have suggested underlying mechanisms of action for Bacopa, the biological pathways altered by this plant extract and the potential upstream mediators for these processes.

Cell culture, differentiation of SH-SY5Y cells and Bacopa treatment
The human neuroblastoma cell line, SH-SY5Y, was purchased from the American Type Culture Collection (ATCC CRL-2266). Cells were cultured in DMEM/F12 (Sigma) supplemented in 10% (v/v) fetal bovine serum (FBS, Gibco) and 1% (v/v) penicillin/streptomycin (P/S, JR Scientific). This medium will be referred to as Complete Medium henceforth. Subculturing was performed as per manufacturer's instructions (ATCC). In brief, as the SH-SY5Y cells grow as a mixture of floating and adherent cells, care was taken to ensure the floating cells in the medium were collected and recovered by centrifugation. These collated floating cells would be combined with trypsinized adherent cells and subcultured. Cells were also passaged less than three times to ensure that the cells remained neuroblast-like [29] (Fig 1A and 1C). For experiments involving undifferentiated SH-SY5Y cells, the plating density was 0.4 x 10 6 cells/cm 2 . To differentiate SH-SY5Y cells, they were plated at a density of 0.5 x 10 6 /cm 2 on culture surfaces coated with 10 μg/ml laminin (Sigma) and maintained in Complete Medium for 18 h. After which, they were maintained in serum-free Complete Medium. 50 nM of human insulin-like growth factor-I (IGF-1) (Sigma) was added to promote differentiation [30]. 48 h after the switch to serum-free Complete Medium and the addition of IGF-1, the medium was replenished. Bacopa treatment was carried out 72 h after the start of differentiation. Undifferentiated and differentiated cells were treated with 3 μg/ml Bacopa for 24 h or 10 μg/ml Bacopa for 4 h, or with vehicle controls. For all experiments, we used a standardized extract of Bacopa (CDRI-08), containing no less than 55% bacoside A and bacoside B as its bioactive components that was extracted by ethanol extraction (Laila Impex, Vijaywada, India) [31,32].
High (10 μg/ml) and low concentration (1 μg/ml) of Bacopa were supplemented during the addition of H 2 O 2 for investigation of neuroprotective effects of Bacopa.
RNA sample preparation, library construction and RNA sequencing RNA was harvested from treated SH-SY5Y cells using the NucleoSpin RNA/Protein kit (Macherey-Nagel GmbH). Library construction and RNA sequencing were performed by the Duke-NUS Genome Biology Facility (DGBF). Briefly, 2.2 μg of RNA was sent for library construction. Quality of RNA was analyzed with the Agilent Bioanalyzer and RNA with RIN > 9 was used. Following poly-A enrichment, recovered RNA was processed using the Illumina TruSeq RNA Sample Preparation Kit v2 protocol (non-stranded) to generate adaptor-ligated libraries. A total of six samples were sequenced, obtained from undifferentiated and differentiated SH-SY5Y cells treated with (i) vehicle-treated control, (ii) 3 μg/ml Bacopa for 24 h, and (iii) 10 μg/ml Bacopa for 4 h. Samples were sequenced using two lanes of an Illumina HiSeq2000 sequencer using 76-bp paired-end reads.
Computational analyses of RNA-sequencing data RNA-Seq data was mapped to the human genome using Partek Flow (version 4.0.15.0406) and Partek Genomics Suite (version 6.15.0327). After the adaptor sequences were trimmed away, reads were mapped to the Homo sapiens genome (hg38) with TopHat2 (version 2.0.8). Local alignment was performed on the unaligned reads from TopHat2 to the human genome (hg38) with Bowtie2 (version 2.1.0). Aligned reads from the TopHat2 and Bowtie2 alignment were combined in Partek Flow. Post-alignment QA/QC was performed after each alignment step and aligned reads had an average quality Phred score above 30. The unique paired reads were used for gene expression quantification. Reads were assigned to individual transcripts of a gene based on the Expectation/Maximization (E/M) algorithm [33]. In the Partek Genomics Suite software, the E/M algorithm was modified to accept paired-end reads, junction aligned reads, and multiple aligned reads if these are present in the data. RNA expression was calculated as fragments per kilobase of transcript per million mapped reads (FPKM) values of the human RefSeq genes for paired-end sequencing. To identify differentially expressed genes, Partek's Gene Specific Analysis (GSA) algorithm was used. Read counts between samples were normalized with the Upper Quantile method and analysis was performed at the transcript level. A cutoff value of multimodal P < 0.05 and fold change > 2 or < -2 were set. A gene ontology analysis was conducted using Partek Genomics Suite.

Functional class scoring using gene-set enrichment and overrepresentation analysis
To identify biological functions affected by differential gene expression, evidence for functional class enrichment involving biological pathways or gene-sets was sought via two independent methods. First, pathway enrichment analysis was conducted via the Gene Set Enrichment Analysis tool using the "pre-ranked" option [34]. For this analysis, a total of 10744 genes with a FPKM ! 5 in at least 1 sample were included, and ranked by their fold-change in expression between the control and Bacopa treated samples. The ranked gene-list was used to query pathways from Gene Ontology-Biological Process and Gene Ontology-Molecular Function, (downloaded from the Molecular Signatures Database, MSigDB, [35], as well as custom pathways (Reactome pathways plus user-defined gene-sets for brain-specific functions) for statistically significant enrichment of higher or lower ranked genes. Significance estimates were adjusted for multiple testing via the false discovery rate (FDR).
Over-representation analysis (ORA) of biological functions and putative upstream regulators was carried out by subjecting a pre-filtered list of 576 differentially expressed genes (FPKM 5 in at least 1 sample and absolute fold-change ! 1.5) to QIAGEN's Ingenuity Pathway Analysis tool (IPA, QIAGEN Redwood City, https://www.qiagenbioinformatics.com/ product-login/). First, reference gene-sets corresponding to 'biological functions' (as defined in the Ingenuity Knowledge Base), were analyzed via Fisher's exact test for statistically significant over-representation in the list of differentially expressed genes. Additionally, predictions of changes in the activity status of 'upstream regulators' (specifically, transcription factors), that could putatively explain the observed gene expression changes due to Bacopa treatments, was also carried out. An ORA was first performed to determine whether an upstream regulator was enriched for differential expression of its target genes (the list of regulators and their target genes were, again, obtained from the Ingenuity Knowledge Base). The overall activation or inhibition status of the regulator was then inferred from the degree of consistency (up-or down-regulation) in the expression patterns of its target genes, expressed as a z-score. Regulators with z ! 2 or z ! −2 were considered to be activated or inhibited, respectively. cDNA synthesis and quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) RNA was harvested using the NucleoSpin RNA/Protein kit (Macherey-Nagel GmbH). The amount of RNA was measured spectrophotometrically using a Nanodrop ND-1000 Spectrophotometer. cDNA was generated from 2 μg of RNA using the iScript cDNA synthesis kit (Bio-Rad Laboratories). Random hexanucleotides were annealed for 5 min at 25˚C. cDNA synthesis was performed for 30 min at 42˚C, followed by an enzyme inactivation step for 5 min at 85˚C. cDNA was stored at -20˚C until use. 1 μl of the cDNA reaction mix was used for qRT-PCR, which was performed using iQ SYBR green reagents (Bio-Rad Laboratories) on the iQ5 Multicolor Real-Time PCR Detection System (Bio-Rad Laboratories) with the following PCR profile: 95˚C for 3 min, 40 cycles of 95˚C for 10s and 55˚C for 30 s. After the completion of the PCR, melt curve analysis was performed using the following paradigm: 95˚C for 1 min, 55˚C for 1 min followed by ramping up the temperature from 55˚C to 95˚C. RPL19 was used as a control.

Characterization of SH-SY5Y cells
The goal of these experiments was to characterize the effects of Bacopa on functional properties and gene expression profiles of SH-SY5Y human neuroblastoma cells. In the presence of serum, undifferentiated SH-SY5Y cells can be propagated indefinitely, while they can be made to terminally differentiate by withdrawing serum, plating on a properly coated surface and addition of differentiating factors. SH-SY5Y cells (passage number 27, www.atcc.org) were differentiated by plating them on laminin-coated coverslips, using a protocol optimized by Dwane et al [30], which consisted of removing fetal bovine serum (FBS) and adding 50 nM Insulin Growth Factor 1 (IGF-1). Cells showed a differentiated neuronal phenotype (formation of extensive neurites) after 3-5 days (Fig 1B). In order to obtain a better visualization of the individual cells, we transfected the cultures with a cDNA encoding Green Fluorescent Protein (GFP) and used fluorescence microscopy to document the differentiation (Fig 1D and  1E).

RNA-Seq analysis
In order to understand the changes in gene expression levels caused by Bacopa, we performed a comprehensive RNA-Seq analysis for both undifferentiated and differentiated SH-SY5Y neuroblastoma cells. Two different concentrations and treatment times were used: 10 μg/ml for 4 hours and 3 μg/ml for 24 hours. These concentrations were determined in pilot experiments to be the highest concentrations that were not deleterious to the cells. Vehicle (DMSO) was used as a control for the Bacopa effect. In addition, we evaluated the effect of differentiation itself. For each condition, the sequence of 50 million mRNA fragments ("reads") was obtained which were then mapped to the human genome, as described in the Methods section, resulting in a relative abundance for all mRNAs expressed at a reasonable level. For each "treatment", a list of differentially expressed mRNAs was generated, using a P-value smaller than 0.05 and an absolute fold-change of at least 2. Table 1 summarizes the results.

Effect of differentiation
Differentiation of SH-SY5Y cells was accomplished by growing them on laminin and replacing serum with IGF-1. Cells showed a differentiated phenotype (formation of extensive neurites) after 3-5 days (Fig 1B, 1D and 1E). The RNA-Seq analysis detected 31,500 different transcripts in SH-SY5Y cells. Due to alternative splicing, the number of distinct mRNAs in a cell is typically larger than the number of genes (~23,000) in the human genome. The differentiation protocol resulted in a change of 502 mRNA levels, of which 78 could be classified as transcribed from 'neuronal' genes (S1 Table). The results indicated that (1) our differentiation protocol was effective in altering the gene expression profile towards a more neuronal phenotype, and (2) the RNA-Seq approach can be used to identify changes in mRNA levels in an un-biased manner and at a genome wide scale. Table 1 summarizes the effect of four different Bacopa treatment protocols (2 concentrations/ durations in undifferentiated and differentiated SH-SY5Y cells) on the mRNA profile of the neuroblastoma cells. The 4-hour treatment with Bacopa altered the expression of more mRNAs than the 24-hour Bacopa treatment: 57 and 66 vs. 37 and 29 altered mRNAs in undifferentiated and differentiated cells, respectively (Table 1). This indicated that the effect of Bacopa on gene transcription seen after 4 hours subsides after one day for many of the affected transcripts. 4 hours of Bacopa on differentiated cells was most effective, with 66 mRNAs being altered at least 2-fold (Table 1). This data set was therefore selected for further analysis. A gene ontology analysis was conducted using Partek Genomics Suite. Fig 2 illustrates  We considered mRNA levels to be altered if the P-value was smaller than 0.05 and the absolute value of the fold change was larger than 2. Four hours of Bacopa on differentiated cells was most effective (highlighted in bold and italics). listing the biological processes, cellular components and molecular functions affected by Bacopa treatment. Many of the affected biological processes referred to functions that were of critical importance in the central nervous system. Table 2 summarizes the mRNAs altered by Bacopa treatment encoded by genes with a neuronal function. The most striking finding was the Neuroplastin gene (NPTN), because a single nucleotide polymorphism (SNP) in the Neuroplastin locus associates with cortical thickness and intellectual ability in adolescents [36]. Neuroplastin is a synaptic glycoprotein involved in long-term potentiation (LTP) in hippocampal CA1 synapses that modulates neuritogenesis and neuronal plasticity [37,38].

Effect of Bacopa treatment: Gene level analysis
The analysis described above was based on mapping the reads to individual transcripts of the human genome (see Materials and methods). Because of the extensive alternative splicing seen for many genes, a process that is most prominent in the brain, mapping the reads is a difficult process and bound to be error-prone. We have therefore used an additional analysis, which is more robust, in which the reads mapped to all exons belonging to a gene were combined to provide a gene level statistic. Several additional genes regulated by Bacopa were identified this way, as summarized in Table 3. Effect of Bacopa on gene expression in human neuroblastoma cells Validation of the RNA-sequencing results: qRT-PCR Next, we validated the RNA-Sequencing results using qRT-PCR experiments for a subset of genes. Fig 3 summarizes the results by showing the absolute value of the fold change produced by Bacopa treatment. In the case of genes with alternatively spliced transcripts, PCR primer pairs were designed that were unique for the mRNA isoform that was altered by Bacopa. In some cases, this proved to be hard, since the exon that was unique for the altered transcript was very short. A case in point is Neuroplastin, which has four mRNA variants, one of which is affected by Bacopa (NPTN_D). The distinguishing feature of the affected variant is the absence of exon 2 and that exon 7 is shortened by 12 nucleotides. We were able to validate the effects of Bacopa on the transcripts identified by the gene-level analysis (ANKRD1, SLC7A11, CHAC1, TXNIP, YPEL4 and STON1) ( Table 3). However, for the mRNAs for which only one or a few of the alternatively spliced variants were affected by Bacopa, the validation was less successful: the effects were either much smaller than those seen in the RNA-Seq analysis (HBA1 and HBA2) ( Table 3) or there was no measurable effect (NPTN_D) ( Table 2). Effect of Bacopa on gene expression in human neuroblastoma cells

Gene Set Enrichment Analysis (GSEA)
The deep sequencing experiments (RNA-Seq) we have performed yielded a set of genes that were differentially affected by Bacopa treatment. In order to better understand the underlying biological processes affected by Bacopa, we compared the genes affected by Bacopa to the gene sets belonging to each of the entries in the Gene Ontology database. A statistical test was performed for each GO term to investigate if it is enriched for the Bacopa-regulated genes. The Broad Institute has developed a set of tools to query the Molecular Signatures Database (MSigDB), a collection of annotated gene sets. Table 4 shows the results of a GSEA of three databases: Reactome, Gene Ontology (Molecular Function: GOMF, Biological Process: GOBP) and Canonical. Twenty-five gene sets were identified with a significant enrichment score (Pval < 0.05). Twenty-two of these fell into three broad functional categories: (1) oxidative stress response, (2) translation regulation, and (3) membrane transport. Fig 4 provide more detailed results for the 'Oxidative Stress Response' gene set, including the Enrichment plot (Fig 4A), the genes enriched in this pathway ( Fig 4B) and the Mean-Average (MA) plot (Fig 4C).

Ingenuity Pathway Analysis (IPA)
With the current results, we were interested to understand the possible upstream biological causes and potential downstream effects. As such, we performed IPA on the RNA-Seq data for Bacopa treatment of differentiated SH-SY5Y cells. IPA infers statistically, based on the altered gene expression in our dataset, the possible downstream effects on biological functions and association with diseases (Table 5) and Upstream Regulators (Table 6). Table 5 shows four Biological Functions identified by IPA that were over-represented with differentially expressed genes in Bacopa-treated vs. Control samples: (1) carbohydrate metabolism, (2)   to be inhibited by Bacopa treatment, including (1) organismal death, (2) damage of brain, (3) apoptosis and (4) synthesis of reactive oxygen species (ROS), and (5) growth failure. Note that the activation states of these disease-related categories were decreased by Bacopa treatment. Many of these Bacopa effects would therefore be expected to contribute to increased brain health and improved brain development. Table 6 shows the three Upstream Regulators (transcription factors) that the IPA analysis predicted as differentially regulated based on the observed effects of Bacopa treatment on gene transcription. The transcription factors, Activating Transcription Factor 4 (ATF4) and Nuclear Factor, Erythroid 2 Like 2 (NFE2L2) were predicted to be activated, while Forkhead Box O3 (FOXO3) was inhibited. ATF4 is a transcription factor involved in the endoplasmic reticulum (ER) stress response [39]. It is known to collaborate with NFE2L2 (aka NRF2), a master redox switch that turns on cellular signaling involved in the induction of cytoprotective genes [40][41][42][43]. FOXO3 is a 'forkhead box' transcription factor that regulates autophagy and neural oxidative stress-mediated stem cell homeostasis and has been associated with longevity in worms [44][45][46][47][48]. Therefore, all three upstream regulators identified by IPA are transcription factors involved in regulating cellular stress pathways. regulator IPA identified three transcription factors that each regulated the transcription of genes whose mRNA levels were affected by Bacopa treatment. This suggests that Bacopa exert its effect by inhibiting FOXO3 and activating both NFE2L2 (NRF2) and ATF4 (CREB2). The transcription level of FOXO3 and NFE2L2 was not significantly (n.s.) changed by Bacopa, but they were predicted to be functionally inhibited and activated, respectively, by the treatment, To validate the findings obtained from the RNA-Seq data, we have performed a functional oxidative stress assay of Bacopa and studied the neuroprotective properties of Bacopa. SH-SY5Y cells were first challenged with hydrogen peroxide (H 2 O 2 ) and the effect was evaluated using a fluorescence-based live-death assay. Fig 5A and 5B illustrate the concentration-survival curve for H 2 O 2 . H 2 O 2 was quite efficacious in causing cell death. Bacopa was able to partially protect against H 2 O 2 toxicity, at a concentration of 10 μg/ml, while the lower dose of 1 μg/ml showed no efficacy in this neuroprotection assay (Fig 5C). Protection was seen at a lower dose than previously reported for Bacopa in SH-SY5Y cells, in which the lowest dose tested was 25 μg/ml [49]. This neuroprotective property of Bacopa seen in these functional oxidative stress assays validated the results and analysis obtained from the RNA-Seq, which suggested an effect of Bacopa on the oxidative stress response pathway.

Discussion
By identifying changes in the gene transcription profile in human neuroblastoma cells (SH-SY5Y) using an RNA-Seq approach, we have identified some of the molecular/cellular components and pathways that underlie the mechanism of action of Bacopa Monnieri (Bacopa). Changes in mRNA levels produced by Bacopa treatment were validated using qRT-PCR. A gene ontology analysis of the mRNAs affected by Bacopa indicated effects on brain development, Ca 2+ and K + ion homeostasis, synaptic function and long-term potentiation (Fig 2), providing the first mechanistic support for the well-published nootropic effects of Bacopa on cognitive function and memory performance. The Gene Set Enrichment Analysis (GSEA) performed for the effect of 4 hours of Bacopa treatment on differentiated SH-SY5Y neuroblastoma cells identified three biological functions/pathways affected by Bacopa: 1) oxidative stress response, 2) transmembrane transport and 3) translation regulation ( Table 4). The first two pathways were upregulated by Bacopa and the third was downregulated, suggesting an inhibition of general translation. Interestingly, aminoacylation of tRNAs was upregulated by Bacopa, which would indicate that the levels of tRNAs will increase, thereby stimulating translation rates. It is possible that the overall effect of Bacopa treatment is an introduction of a selective bias in the translation of a specific class of mRNAs. This would have to be investigated further.
A second interesting biological function uncovered by the GSEA was the up-regulation by Bacopa of a select group of members of the SLC family of solute carriers, which consists of over 300 proteins functionally grouped into 52 sub-families, including facilitative transporters, primary and secondary active transporters, ion channels, and the aquaporins. Bacopa increased the expression of 31 (out of 300) SLC family members. Table 7 summarizes the function of some of these transmembrane carriers by listing the molecules they transport. Many have critical functions in the central nervous system. For instance, Glutamate, L-DOPA, norepinephrine and monoamines are neurotransmitters, or their precursors. The Na + /Ca 2+ exchanger and K-Cl co-transporter are critical for intracellular ion homeostasis and neuronal excitability. Zinc is a co-factor that regulates NMDA receptor function. Finally, there are transporters for glucose, fatty acids, phosphate and many amino acids, which are required for normal metabolism. These transporters are heteromeric complexes assembled from multiple subunits. Bacopa altered the mRNA levels for two and four subunits for the complexes that transport glucose and Zinc respectively.
Moreover, Bacopa was found to upregulate many genes that respond to oxidative stress, thereby likely improving the capacity of the cell to properly deal with such insults. Oxidative damage resulted from reactive oxygen species caused pathological manifestations of aging resulting in cognitive dysfunction [50]. Overexpression of superoxide dismutase in aged mice exhibited enhanced hippocampal LTP, better cerebellum-dependent motor learning and better hippocampus-dependent spatial learning [51]. Similarly, rats exposed to another antioxidant, Curcumin, also had improved memory retention [52]. Identification of the exact subset of genes of the oxidative stress response pathway that are controlled by Bacopa (Fig 4) now enables more mechanistic studies of the neuroprotective effect of Bacopa. This important Effect of Bacopa on gene expression in human neuroblastoma cells finding goes a long way towards understanding the mechanism underlying Bacopa's neuroprotection capabilities that others and we have characterized (Fig 5) [24, 25, [53][54][55][56][57][58].
Supporting information S1