Animals

For the RNA-Seq studies adult male Sprague-Dawley rats (270 g–370 g) obtained from Charles River Laboratories (Wilmington, MA) were maintained under normal laboratory conditions (temperature: 21–23°C, 12 h light-dark cycles with light on at 6:00 AM) with access to unlimited food and drinking water. All procedures were carried out in accordance with the National Institutes of Health (NIH) guidelines on the care and use of animals and an animal study protocol No. 1278 approved by the National Institute of Neurological Disorders and Stroke (NINDS) Animal Care and Use Committee.

For the microarray and qRT-PCR experiments all procedures were approved by the University of Bristol Ethical Review Committee and were carried out in accordance with the Animals (Scientific Procedures) Act, 1986. Adult male Sprague-Dawley rats (10–12 weeks old; Harlan Sera-lab, Loughborough, UK) were maintained in standardized temperature (22±1°C), humidity (50±5%) and diurnal conditions (10 hours light, 14 hours dark; lights on at 0700). All animals were allowed to acclimatize for 5 days prior to the start of the experiment at which point animals (n = 5) were given 2% saline (w/v) for 5 days. Control animals had normal access to drinking water and all animals had access to regular chow ad libitum. Sacrifice time was always between 1100 and 1300 each day.

Chronic Salt loading

To induce chronic hyperosmolality (hypernatremia), rats were given 2% NaCl solution ad libitum as their only drinking fluid for 5 d. Euhydrated (Control) rats had free access to drinking water over the same time period. Daily fluid intake was monitored. Rats were weighed, euthanized by decapitation and their brains removed and quickly frozen on dry ice (see below). The control rats were killed at the same time as the salt loaded rats. Rats tolerate salt loading very well and no special care is necessary. This method is standard in the field to produce systemic hyperosmolarity & hypernatremia. Some papers describing the use this method have previously been published [15,22,23,24,25,26,27]. One of the consequences of chronic salt loading is a nearly two-fold increase in cell volume of the MCNs [15], confirmed in this study) which is known to be accompanied by a global increase in gene expression in the MCNs [7]. Note also the doubling of the total

Tissue Preparation for Laser Capture Microdissection (LCM)

The rats were euthanized by decapitation, their brains were quickly removed from the skull, and the cerebellum and brain stem were removed and discarded. The frozen brain was then trimmed at the level of the anterior commissure prior to sectioning. The brain block containing the hypothalamus was then rapidly frozen on dry ice and stored at −80°C until sectioning was done. Details about this procedure were previously published [28,29]. Briefly, the frozen brain was mounted onto Tissue-Tek OCT compound (Sakura Finetek USA, Torrance, CA) on a chuck with the caudal portion of the brain placed down on the chuck, and the brain on the chuck was placed in the cryostat (Reichert-Jung 2800, Frigocut, Heidelberg, Germany) for twenty minutes to equilibrate the brain with a chamber temperature (CT) of −20°C and an object temperature (OT) of −16°C. Four 12 μm thick sections of the SON region were cut and placed on polyethylene naphthalate (PEN) Membrane Frame slides (Life Technologies, Grand Island, NY, Cat. # LCM0521) under RNAase-free conditions. Four brain sections were usually mounted on a single PEN membrane slide, and between 25 and 30 slides were made from a single brain specimen for up to 120 total sections. Slides were stored immediately in a slide box embedded in dry ice, and the slide box was placed in the −80°C freezer until the slides were used for LCM (slides were used within one month of being sectioned).

Laser Capture Microdissection (LCM)

Slides were kept in their slide box on dry ice until the dehydration step for visualization on the LCM instrument. Before LCM, the sections were fixed and dehydrated with ethanol and xylene under RNase free conditions, as described previously with some modifications [30]. The slides were thawed for 30 seconds in 75% ethanol and further dehydrated by sequential immersion in 95% ethanol for 5 seconds, twice in 100% ethanol for 5 seconds, and followed by xylene for 30 seconds. The slides were left in a hood to allow the xylene to evaporate for 3 minutes. All solutions were made with DEPC water. Three slides can be placed in the Arcturus XT LCM instrument at one time, and therefore, only three slides were dehydrated for microdissection at one time. The rest of the slides were stored on dry ice until they were needed. The ethanol fixation technique was found to provide the best visualization of the SON under the LCM microscope (S1 Fig) and the highest RNA yield and quality (Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Properties of Normosmotic and Salt-Loaded SON samples used in RNA-seq study^a.

https://doi.org/10.1371/journal.pone.0124523.t001

LCM was performed using the Arcturus XT (Life Technologies, Washington DC 20001) laser-microdissection system, immediately after fixation and dehydration of the slides. The SONs were microdissected using CapSure Macro LCM Caps (Life Technologies, Cat. # LCM0211). The PEN membrane frame slide was loaded on the instrument with the tissue side down and a clean glass slide underneath it acting as support [28]. The rat SON is completely visualized in the 20X objective view of the Arcturus XT LCM microscope which allowed for the microdissections of multiple SONs from many serial sections from the same brain specimen onto the same cap.

RNA Isolation and Quality Assessment

The RNA was isolated using the Picopure RNA Isolation Kit (Life Technologies, Cat. #KIT0204). The Macro cap with the microdissected SONs was immediately seated into a 0.5 mL tube containing 50 μL of PicoPure Extraction buffer and incubated for 30 minutes at 42^°C. When multiple caps for an individual rat were collected, the buffer from each cap containing the lysed cells was pooled in a 1.5 mL microcentrifuge tube before continuing with the isolation. After isolation, the RNA’s concentration was measured using a Nanodrop spectrophotometer (ND-1000, Thermo Scientific, Wilmington, DE). The RNA samples were analyzed using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA 95051) to obtain RNA integrity numbers (RIN numbers) as measures of their quality (Schroeder et al, 2006). Three RNA samples from Control SONs and three RNA samples from Salt-loaded SONs met our RIN quality criterion of >6.5, and these RNAs were taken for RNA-Seq analysis. See Table 1 for RNA yield and RIN numbers for these samples.

RNA-Seq Methods

Amplified cDNA libraries were prepared from the isolated RNA samples shown in Table 1 and sequenced by the NIH Intramural Sequencing Center (NISC), Bethesda, MD under the direction of Baishali Maskeri and Alice Young. Briefly, mRNA libraries were constructed from 50–750 ng mRNA using the Illumina TruSeq RNA Sample Prep Kits, version 2. Equal quantities of the the resulting cDNAs for the individual control and SL samples were fragmented using a Covaris E210. Library amplification was performed using 12–15 cycles to minimize the risk of over-amplification. Unique barcode adapters were applied to each library. Equal volumes of individual libraries were pooled and run on a MiSeq. The libraries were then repooled based on the MiSeq demultiplexing results and sequenced on a HiSeq 2000 with ver 3 flow cells and sequencing reagents. At least 4x40 million 100-base read pairs were generated for each individual library. Data was processed using RTA 1.13.48 and CASAVA 1.8.2, providing for four sets of. fastq files per library. Reads contained within these files were subjected to sixteen analysis steps (See S2 Fig and legend). For the first analysis step, 3' adaptor sequence was removed from a read if present using the FASTQ/A Clipper tool (http://hannonlab.cshl.edu). For the second analysis step, reads were import into the CLCbio Genomics Workbench (www.clcbio.com) and quality inspected. For the third analysis step, 15nt from the 5' end of each read was globally trimmed in order to remove nucleotide bias and the last occurring nucleotide at the 3' end was removed due to low quality. Nucleotides with a call accuracy < 95% were also removed. While, read pairs with more than two ambiguities in at least one read were discarded along with pairs having at least one read < 15nt in length. For the fourth analysis step, remaining read pairs were aligned to the rat genome (rn5) for each library using the Workbench’s "RNA-Seq" tool under default parameters. For the fifth analysis step, the number of aligned read pairs to each known gene for each library were exported from the Workbench and imported into R (http://www.r-project.org/). In R, the number of aligned read pairs per known gene per library were converted to "Reads Per Kilo base of transcript Million mapped reads" ("RPKM") values. These RPPKM values were then pedestalled by a value of two (RPKM+2) to control for low-end variance. For the sixth analysis step, pedestalled RPKM values were transformed (Log2(RPKM+2) to promote normality then used to identify and remove non-informative genes (i.e., those not having at least one library with a transformed RPKM value >1). For the seventh analysis step, transformed RPKM values were quantile normalized (Quantile(Log2(RPKM+2))) to correct for library-to-library differences in distribution location, spread, and skew. For the eighth analysis step, normalized RPKM values were explored across libraries via Tukey box plots, covariance-based principal component analysis (PCA) scatter plot, and Pearson correlation-based heat maps to ensure absence of outliers. For the ninth analysis step, normalized RPKM values were organized by library class (i.e., salt-loaded control) and the coefficient of variation (CV) and mean normalized RPKM were calculated for each gene. Locally weighted scatterplot smoothing (LOWESS) was then applied by library class (CV ~ mean normalized RPKM) and the resulting fits plotted. For the tenth analysis step, the LOWESS fits were visually inspected to identify the mean normalized RPKM at which the linear relationship between mean normalized RPKM (i.e., signal) and CV (i.e., noise) is grossly and concordantly lost. This mean normalized RPKM value was then defined as the “confidence criterion”. For the eleventh analysis step, genes not having at least one library with a normalized RPKM greater than the "confidence criterion" were discarded as noise-biased. Normalized RPKM values for remaining genes were floored to equal the "confidence criterion" value if less. For the twelfth analysis step, the Welch-modified t-test was applied to the post-floored RPKM data on a gene-by-gene basis under Benjamin–Hochberg (BH) False-Discovery Rate (FDR) Multiple Comparison Correction (MCC) condition. This provided for one corrected p-value per gene. This resulting p-value describes the probability that the mean difference occurring between the library classes is due to chance. For the thirteenth analysis step, genes were filtered to keep only those having a corrected p-value < 0.10 and an absolute difference of means between library classes >/ = 1.5X. Genes kept were deemed to be those having differential expression between salt-loaded and control. For the fourteenth analysis step, floored RPKM values for these differential expressed genes were used in the generation of a covariance-based principal component analysis (PCA) scatter plot and Pearson correlation-based heat map to confirm good intra and inter library grouping by class. For the fifteenth and sixteenth analysis steps, the Ingenuity Pathway Analysis (IPA) tool (www.ingenuity.com) was used. Specifically, symbols for the genes deemed to have differential expression between salt-loaded and control were import into IPA and the corresponding enriched pathways and functions identified (uncorrected Fishers Exact Test p-value < 0.05).

Fig 1 describes the quality of the RNA-Seq data. There is good separation of the control and SL samples and Fig 1A illustrates the uniformity of the distribution of expression which indicates no sample outliers. Fig 1B illustrates the good within and between sample level groupings indicating neither experiment-level nor cohort-level outliers based on magnitude of expression. Fig 1C corroborates conclusions in 1B but is based on pattern of expression, and Fig 1D shows a heat map based on clustered differential expression.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. RNA-Seq Analysis Results.

(A) Tukey box plot comparing the sample level distributions of gene expression data post normalization (Quantile(Log2(RPKM+2))) for 6 samples (3 controls, green-filled; 3 salt-loaded samples, red-filled). Distributions depicted consist of 26,313 gene expression measurements each. (B) covariance-based Principal Component Analysis (PCA) scatter-plot depicting relationships across 6 samples (3 controls, green-filled; 3 salt-loaded samples, red-filled) when normalized gene expression data (Quantile(Log2(RPKM+2))) is used for the 552 genes which are identified to have significant differential expression between salt-loaded (SL) condition and control. (C) correlation-based unclustered Heat-Map depicting relationships across 6 samples (3 controls, green-outlined; 3 salt-loaded samples, red-outlined) when normalized gene expression data (Quantile(Log2(RPKM+2))) is used for 552 genes identified to have significant differential expression between salt-loaded (SL) condition and control. (D) correlation-based clustered Heat-Map depicting relationships across 6 samples (3 controls, green-outlined; 3 salt-loaded samples, red-outlined) when normalized gene expression data (Quantile(Log2(RPKM+2))) is used for 552 genes identified to have significant differential expression between salt-loaded (SL) condition and control. Results reveal absence of outliers and excellent within and between sample class grouping.

https://doi.org/10.1371/journal.pone.0124523.g001

Raw fastq files for all libraries can be found in the Short Read Archive (SRA; http://www.ncbi.nlm.nih.gov/Traces/sra/; SRP049482).

Microarray Methods

CEL files were generated by the Affymetrix AGCC program for two separate microarray batches, one file per library, These files were subjected to twelve analysis steps (See S4 Fig and legend). For the first analysis step, all files were imported into the Affymetrix Expression Console (http://www.affymetrix.com) and Robust Multi-Array (RMA) expression summarization and normalization performed. This provided Log2 expression values for 31,099 gene fragments per library. For the second analysis step, expression data for all libraries was imported into R (http://cran.r-project.org/) and baseline subtraction performed to correct for microarray batch differences using the common library run between batches. For the third analysis step, batch-corrected expression values by library were explored via Tukey box plots, covariance-based principal component analysis (PCA) scatter plots, and Pearson correlation-based heat map to ensure absence of outliers. For the fourth analysis step, batch-corrected expression values were organized by library class (i.e., salt-loaded, control) and the coefficient of variation (CV) and mean expression calculated for each gene fragment. Locally weighted scatterplot smoothing (LOWESS) was then applied by library class (CV ~ mean expression) and the resulting fits plotted. For the fifth analysis step, the LOWESS fits were visually inspected to identify the expression value at which the linear relationship between mean expression (i.e., signal) and CV (i.e., noise) is grossly and concordantly lost. This expression value was then defined as the “confidence criterion”. For the sixth analysis step, gene fragments not having at least one library with an expression value greater than the "confidence criterion" were discarded as noise-biased. Expression values for remaining gene fragments were floored to equal the "confidence criterion" value if less. For the seventh analysis step, the Welch-modified t-test was applied to the post-floored expression values on a gene fragment-by-gene fragment basis under Benjamin–Hochberg (BH) False-Discovery Rate (FDR) Multiple Comparison Correction (MCC) condition. This provided for one corrected p-value per gene fragment. This resulting p-value describes the probability that the mean difference occurring between the library classes is due to chance. For the eighth analysis step, gene fragments were filtered to keep only those having a corrected p-value < 0.05 and an absolute difference of means between library classes >/ = 1.5X. Gene fragments kept were deemed to be those having differential expression between salt-loaded and control. For the ninth analysis step, floored expression values for these differential expressed genes were used in the generation of a covariance-based principal component analysis (PCA) scatter plot and Pearson correlation-based heat map to confirm good intra and inter library grouping by class. For the tenth, eleventh, and twelfth analysis steps, the Ingenuity Pathway Analysis (IPA) tool (www.ingenuity.com) was used. Specifically, the gene fragments deemed to have differential expression between salt-loaded and control were imported into IPA and the corresponding gene symbols identified along with the enriched pathways and enriched functions for those genes (uncorrected Fishers Exact Test p-value < 0.05). Comparison and intersection of results with those obtained by RNA-Seq was accomplished using the IPA-assigned gene symbol.

Raw. CEL files for all microarrays can be found in the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/; WD = GSE3110, SL = GSE65663).

Immunohistochemistry

In separate experiments involving immunohistochemical assays, animals were anesthetized with isoflurane and immediately perfused transcardially with 100 ml of 0.1 M PBS at 10ml/min followed by 200ml of 4% paraformaldehyde (PFA) in 0.1 M PBS, PH 7.4 at 10ml/min. Brains were removed and, if they were soft and pink following the perfusion, they were post-fixed in 1% PFA overnight. Brains were then cryoprotected in increasing concentrations of sucrose in 0.9% saline. Sucrose concentrations were 5%, 10%, and 15%.

Coronal sections (16μm) were made on a cryostat (Reichert-Jung 2800; Frigocut, Heidelberg, Germany) and mounted onto coated slides (Fisher Scientific 12-550-19, Pittsburgh, PA). The sections were rinsed with 0.3% Triton-X in 0.1M PBS for 5 minutes, followed by two rinses of PBS of 5 minutes each. The sections were then incubated overnight at 4°C with a pan-specific anti-neurophysin mouse monoclonal antibody, PS 45, [31,32] at a dilution ratio of 1:100 in 0.1M PBS. After incubation in primary antibody the slides were placed on a shaker table at room temperature for 2 hours, rinsed 3 times with 0.1M PBS, and then incubated in secondary antibody, Alexa Flour 594-conjugated goat anti-mouse (Molecular Probes, Eugene, OR) at 1:500 dilution for two hours. After completing the immunohistochemistry, slides were counterstained with DAPI to detect all the nuclei in the section, rinsed three times in 0.1M PBS, and then coverslipped using Prolong Gold Antifade reagent (Molecular Probes) as a mounting medium.

In order to visualize the stained sections an epifluorescent microscope (Nikon Eclipse E400) was used. Representative images were then recorded using the accompanying camera (Q-Imaging RETIGA EXi Fast, Cooled Mono 12-bit). Optimizing the exposure time of the photograph to best display the representative images was done using Image-Pro Plus 5.1 (MediaCybernetics, Inc) imaging software. An illustration of an SON section showing PS 45 positive MCN cells, and MCN and non-MCN positive DAPI stained cells is shown in S3 Fig, and quantitative comparisons of these phenotypes is shown in S1 Table.

Quantitative PCR validation of targets

Quantitative PCR methods were as described previously [33]. Brains from both the euhydrated and the 5-day salt loaded rats were carefully extracted and flash-frozen with powdered dry ice. The brain was mounted and 60μm sections stained with Toluidine blue (Sigma Aldrich; 0.1% in ethanol (70%, w/v) were made in order to identify the SON. Once localized, 1mm bilateral punches of the SON were taken from unstained tissue and stored in RNase-free tubes at -80°C.

When all samples were prepared, 1mL QIAzol lysis reagent (Qiagen) was added to each tube and left at room temperature for 5 minutes followed by centrifugation (10 minutes, 10,300rpm, 4 °C) in order to remove cellular debris following lysis. Chloroform was then added to the supernatant homogenate and the samples were centrifuged (15 minutes, 11,200 rpm, 4 °C) in order to extract the nucleic material. Total RNA was then precipitated using one volume of ethanol (70% w/v) and purified using the Qiagen RNeasy Mini Kit (Qiagen) according to manufacturers protocol. Using the QuantiTect Reverse Transcription Kit (Qiagen), cDNA was then prepared for qPCR assessment using primers according to S3 Table (constructed by MWG Eurofins unless otherwise stated) using the Applied Biosystems StepOnePlus Real-Time PCR System and fold-change assessed by establishing the ΔΔCt between the Rpl19 calibrator gene and the target gene.

Characterization of the laser microdissected SON

The SON is composed of a single neuronal phenotype, the MCNs, which are large (20–40 um) neurons easily visualized in the Arturus XT LCM microscope (see S1 Fig). The MCNs are further divided into two subtypes, the Oxt- and Avp-MCNs with the Avp neurons predominating in an approximately 2:1 ratio in the SON [2,34]. MCNs are the primary sources of the RNA isolated from the laser microdissected SON. However, there is also a non-neuronal small cell (e.g, glial, vascular, etc.) population in the SON that can be visualized by DAPI-staining of all nuclei in the SON (S3 Fig). By quantitatively comparing the small dark-stained nuclei in the small cell population to the larger, lighter DAPI stained nuclei in the immunohistochemically identified MCNs in the SON we estimate that there are twice as many non-neuronal cells than MCNs in the SON (S1 Table).

Table 2 shows the quantitative determination of all the phenotype-specific RNAs found by the RNA-Seq analysis of the microdissected SONs presented in both quantile normalized Log2(RPKM+2) units and in relative gene expression linear units (in parentheses). The most abundant transcripts expressed in the SON are the the Oxt-and Avp peptide RNAs, whose expression is orders of magnitude greater than the relatively abundant pan-neuronal specific beta tubulin 3 (TUBB3) gene. The Avp RNA expression is about two fold greater than the Oxt-RNA expression in the SON consistent with the greater number of Avp MCNs as compared to Oxt MCNs present in the SON. Other MCN-specific marker genes such as prodynorphin (Pdyn), tyrosine hydroxylase (TH), and the pan-neuron-specific genes such as TUBB3 and neurofilament M (NEFM) are also highly expressed in the SON.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Relative Expression of Cell-Type Specific Markers in SON Samples^a.

https://doi.org/10.1371/journal.pone.0124523.t002

As expected, various non-neuronal cell (e.g, astrocyte, microglial, vascular, etc.) RNAs are also found in the SON sample but at much lower levels. The presence of oligodendroglial-specific RNAs such as myelin basic protein (MBP) RNA could represent a slight tissue contamination from the optic chiasm adjacent to the SON, but is more likely due to the detection of sequences of a related but distinct neuronal protein known as golli-MBP [35,36], previously shown to be significantly expressed in the SON [29].

RNA-Seq Analysis of the SON in control (normosmotic) and SL conditions

A systematic description of our RNA-Seq analysis is presented in S2 Fig (and in Methods). A critical component in this analysis is our approach to set the noise threshold level of RNA-Seq gene expression (see S5A Fig) at or below which no data was accepted as reliable. This meant that most of the data obtained for the 26,313 genes that were detected by the RNA-Seq method before noise level correction were removed from consideration, and in this way probably increased the reliability of the positive data but also increased the likelihood of false negatives. Following noise correction the genes found to be present in the Control and SL SON transcriptomes are shown in S4 Table and S5 Table, respectively. S6 Table presents the total number of 9,709 genes that were expressed in both the control the SL SONs, and S7 Table provides the statistical testing data for these 9,709 genes. Thus, the data in S6 and S7 Tables represents the most reliable identification of the genes that are present in the SON, provides their ingenuity descriptions, and shows their fold changes in expression in response to SL. While these genes are all definitely present in the SON, few significantly change in expression as a result of SL. An example of this is shown in S8 Table which shows that of the 772 of the genes found expressed in the SON and involved in either transcriptional or translational processes (in S7 Table), only about 42 are significantly changed in gene expression by SL The transcriptional and translational regulator genes that do significantly change in gene expression by SL are shown in Tables 3 and 4. Table 3 lists those genes that increase and Table 4 lists those that decrease in expression after SL. Another functional class of genes related to transcriptional regulation are the ligand-dependent nuclear receptor family of genes. A list of the nuclear (hormone) receptor mRNAs that are found in the SON and their changes in gene expression is presented in Table 5. Here too, there were only modest, if any, changes in gene expression of the ligand-dependent nuclear receptors produced by SL. In contrast, Tables 6 and 7 list a wide variety of other genes with diverse functions that are significantly and strongly increased and decreased, respectively, in expression by SL (with greater than two-fold changes).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Transcriptional and Translational Regulator mRNAs increased in expression by SL.

https://doi.org/10.1371/journal.pone.0124523.t003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Transcriptional and Translational Regulator mRNAs decreased in expression by SL.

https://doi.org/10.1371/journal.pone.0124523.t004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Ligand-dependent nuclear receptor mRNAs expressed in SON.

https://doi.org/10.1371/journal.pone.0124523.t005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Other mRNAs in SON that are increased in expression by SL.

https://doi.org/10.1371/journal.pone.0124523.t006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 7. Other mRNAs in SON that are decreased in expression by SL.

https://doi.org/10.1371/journal.pone.0124523.t007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Pie Chart showing Locations and Types of 552 selected, expressed genes (See columns N and O in S9 Table).

https://doi.org/10.1371/journal.pone.0124523.g002

Microarray Analysis of the SON in control (normosmotic) and SL conditions

The strategy used to analyze the microarray data is shown in S4 Fig (and Methods), and the analysis of the noise threshold for the microarray data is illustrated in S5B Fig. Using this approach, 10,382 gene probes were identified as being expressed in control SONs S10 Table), 10,476 gene probes in SL SONs (S11 Table), and 11,293 gene probes in both control and SL SONs (S12 Table). After conversion of the identified gene probes derived from the microarray analysis and the RNA-Seq data to gene symbols, the intersections of these data, based on their unique symbols, were determined and illustrated in the Venn diagram shown in Fig 3. There were 6,466 genes present in the SON that were represented in both the microarray and RNA-Seq data, however 1,040 of the expressed genes were found only in the microarray data, and 2,762 of the expressed genes were selectively found in the RNA-Seq data and not in the microarray data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Venn diagram of intersections of RNA-Seq and Microarray data of post-noise filtered genes representing 7,506 microarray and 9,228 RNA-Seq genes.

Complete data are shown in S4 Table for RNA-Seq genes and S10 Table for microarray genes. Note that 6,466 unique genes are found in both platforms to intersect.

https://doi.org/10.1371/journal.pone.0124523.g003

Comparison of RNA-Seq and Microarray Analyses of specific gene expression changes in the SON as a result of SL

Fig 4 depicts a Venn diagram display of the intersection of RNA-Seq and microarray data after correcting for noise thresholds and statistical analyses of each data set. The RNA-Seq data shown in Fig 4 is derived from the 552 identified genes presented in S9 Table, and the microarray data used is from the 1,309 identified gene probes in S13 Table. The 146 unique genes illustrated in Fig 4 whose expression in the SON is altered by SL were found to intersect between the two different platforms. These correspond to the 170 gene probes in the microarray study (S14 Table) and are listed in S15 Table as the 146 intersecting gene candidates. Fig 4 also shows that in addition to the 146 genes present in the SON that are represented in both data sets, there are 406 expressed genes that are altered in expression by SL that are selectively found in the RNA-Seq data but not in the microarray data (see S16 Table). In addition, there are 1139 expressed gene probes (S17 Table) which convert to 884 unique genes (Fig 4), that are altered in expression in the SON by SL and are found only in the microarray data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Venn diagram of intersections of RNA-Seq and microarray data after correcting for noise threshold and statistical analyses (see S2 Fig, Steps 11 and 12, and S4 Fig steps 6 and 7 for microarray).

There are 552 genes selected from the RNA-Seq data (see S9 Table) and 1030 unique genes selected from the microarray analysis (S13 Table). Note that 146 unique genes are found from both platforms to be unequivocally expressed and altered by salt-loading (See S14 and S15 Tables). Lower panel: Shows X-Y scatter plot of the fold changes in the 146 genes observed between Con & SL samples that were found in both the Microarray (X-axis) and RNA-Seq (Y-axis) data.

https://doi.org/10.1371/journal.pone.0124523.g004

Ingenuity Pathway Analysis (IPA) of networks, enriched functions, enriched pathways for the 552 RNA-Seq genes selected as differentially expressed between SL and Control SONs

Applying Ingenuity Pathway Analysis (IPA) (www.ingenuity.com) to the data shown in S9 Table suggests that certain functional relationships and interactions between the genes in the SON have been altered as a result of SL treatment. S18 Table shows the top 25 Networks proposed by IPA for the 552 RNA-Seq genes that were selected as differentially expressed between SL and Control. SONs. The top two predicted Networks are illustrated in Fig 5 (with Atf4 and Akt as hub genes) and in S6 Fig (with NFkB as a hub gene). Other IPA core analyses suggest predicted enriched functions (S19 Table) and enriched pathways (S20 Table) for the 552 RNA-Seq genes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Network Analysis depicting gene products and known relationships between them for the second-ranked scoring network by Ingenuity (http://www.ingenuity.com/) when provided list of differentially expressed genes observed between salt-loaded and control by RNA-Seq (based on data in S9 Table).

Gene products are represented using circle-shaped symbols with connected edges drawn between them to describe interactions (solid edge = direct interaction, dashed edge = indirect interaction). Color-filled shapes indicate the direction of differential expression observed between salt-loaded and control (green = up, red = down). Circle-shaped symbols not color-filled represent gene products not observed differentially expressed between salt-loaded and control.

https://doi.org/10.1371/journal.pone.0124523.g005

qRT-PCR validations of selected genes in the SON

Table 8 summarizes the results of qRT-PCR data illustrated in S7 Fig for several genes that showed very large changes in expression between control and SL SONs in RNA-Seq (e.g, vgf and procr) and several genes that were identified in the network analysis (see Fig 5). Relevant RNA-Seq data are also shown for comparison in Table 8. Note that all but two genes (Insig1 and Eaf1) were confirmed by qPCR as being significantly regulated by salt loading treatment.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 8. Selected qPCR determinations of fold changes in gene expression due to SL in the SON.

https://doi.org/10.1371/journal.pone.0124523.t008

Alternatively Spliced Genes expressed in the SON under Control and SL conditions

The CLCbio Workbench (www.clcbio.com) was employed to search for evidence of alternative splicing in the RNA-Seq data that occurred during SL. RNA-Seq expression was generated for a total of 27,180 transcripts. After normalization, noise modeling, and filtering, there were data for 8,537 transcripts that underwent statistical testing. Expression levels differed significantly between the SL and control groups for 529 of the transcripts subjected to the IPA analysis. Of these 529 transcripts, 525 are represented by unique genes and 4 mapped to the same gene (2 were mapped to the Argn gene and 2 were mapped to the Nnat gene). This indicates a very modest extent of alternative splicing in the 529 evaluated SON genes. These data are shown in S21 Table.

Graphic representations of these data are shown in S8 Fig for Argn which has 38 exons and in S9 Fig for Nnat which has 3 exons. The only difference between the two Argn transcripts is that the final exon in this gene is missing from one of the sliced variants (see star in S8 Fig). S9 Fig shows that the middle exon (see star) is missing from one of the transcripts of Nnat. Both Nnat transcript forms are expressed in SL and control SONs: expression is higher under SL conditions, and the 2 exon form is more highly expressed than the 3 exon form.

The SON Transcriptome and Effects of SL

RNAs representative of all of the cellular phenotypes expected to be present in the SON are present in the RNA-Seq data obtained from the SON (Table 2), and the MCN- and astrocyte-specific RNAs are predominant in the SON samples. Both the MCN- and astrocyte-specific RNA expressions are significantly altered by SL treatment, but in opposite directions. All the MCN-specific and pan-neuronal RNA markers were specifically and substantially increased by SL, which is comparable to the two-fold increase in total RNA levels in the SON produced by SL (S2 Table). In contrast to the MCN-specific transcripts, all of the astrocyte RNA markers (e.g, GFAP and S100b) were significantly decreased in expression under SL conditions (Table 2). This increase in expression of neuronal markers together with a decrease in expression of astroglial markers such as GFAP are similar to that previously reported in the SON of dehydrated rats [51]. None of the other non-neuronal marker RNAs shown in Table 2 showed statistically significant changes in expression levels with SL. The extremely high levels of the Oxt- and Avp- RNAs in the RNA-Seq data compromised our ability to accurately quantify the changes in expression of these genes during SL. Such high levels of expression produce difficulties in RNA-Seq analysis which are described as transcriptional amplification [52]. However, it is well established that Oxt- and Avp gene expression approximately doubles in the SON during chronic SL [7] and this was confirmed by the qPCR data obtained in this study (see Table 8).

Our RNA-Seq analysis reliably detected 9709 genes in the SON samples. These data for the SON transcriptome are presented in S7 Table together with their fold changes in SL, their Ingenuity symbols and descriptions of the listed genes. Since the SON contains non-neuronal cell phenotypes in addition to the MCNs, it is important to note that the GFAP marker molecule for the most abundant non-neuronal phenotype, the astrocytes, contains one hundreth of the number of transcripts than that found for the Avp mRNAs in the SON (Table 2). As noted above, the MCN-specific markers and pan-neuronal mRNAs increased with SL, whereas the opposite, a decrease in mRNA was found for the astrocytic markers (e.g, GFAP, in Table 2). The abundance of the mRNAs from the other non-neuronal phenotypes (e.g. microglia, oligodendroglia, vascular) were very low and close to the noise level, and showed no significant changes as a result of SL. Therefore, it is likely that most of the increases in gene expression that are found with SL are due to changes within the MCNs. The 552 genes that show statistically significant changes in SL are listed in S9 Table. These 552 genes represent a variety of functions and subcellular locations (Fig 2). Some of these changes are in transcriptional and translational regulator genes (Table 3) and various other genes that were changed more than two-fold by SL treatment (Table 5). The fold changes found in a selected number of genes in this group were confirmed by qPCR (Table 8).

The microarray analysis identified 7,506 genes that were expressed in both control and SL SONs (S12 Table). Comparisons of the unique genes predicted by the microarray data and the genes predicted from the RNA-Seq data to be present in the SON is shown in the Venn diagram in Fig 3. There are 6,466 genes that intersect between these different platforms, and 1,040 genes were found only in the microarray data, and 2,762 genes in the RNA-Seq data were not present in the microarray data (Fig 3). Similar differences in comparisons between RNA-Seq and microarray data analyses of the same experiments have been reported in other studies [48,53]. These differences also may be due, in part, to the difference in sensitivity between the two methods, and the larger linear, dynamic range of RNA-Seq [48,49]. Other possible causes of these differences may be that the SON samples for RNA-Seq were obtained by laser microdissection and the microarray SON samples were obtained by manual dissection of the SON, and the duration of SL was for 5 days in the RNA-seq study and for 7 days in the microarray study. A Venn diagram of intersections of the 552 genes from the RNA-Seq study (S9 Table) and the 1,309 genes from the microarray data (S13 Table) that changed with SL is shown in Fig 4. Of these genes, only 146 were found to intersect and the list of these genes and their fold-changes and ingenuity descriptions are presented in S15 Table. Since these 146 genes identified as definitely changing were confirmed by these two different methods, this represents a type of validation of the reliability of these conclusions. Also in the Venn diagram in Fig 4 there are 406 genes that are solely identified by the RNA-Seq method (S16 Table) and 1139 genes identified only by the microarray method (S17 Table). Many of the latter genes are substantially changed in expression by SL (e.g. with > twofold changes) and with statistically significant p-values, and in some cases have been validated by qPCR. In addition, many of the genes identified in this study as being expressed in the SON have also been validated by other independent methods, for example for Neuronatin(Peg5) and Pspss2 by differential hybridization [54], and for Trpv2, Vgf, and opsin by in situ hybridization or immunohistochemistry [29,55,56].

Biological Significance of the SON Transcriptome Data

One can easily be overwhelmed by the abundance of data obtained from either the RNA-Seq or the microarray platforms. One clear value of these data is that it allows the investigator to determine whether the mRNA of interest is actually present in the tissue being studied (e.g., the SON), its relative abundance, and whether it responds to experimental perturbations (e.g., SL). One approach to organize such data is to analyze the complex transcriptomic and genomic information by using Ingenuity Pathway Analysis (IPA) (www.ingenuity.com). IPA analyzes data from expression profiling and proteomic experiments and generates large gene lists which suggest molecular interactions in the data and networks based on a large database in the literature. S18 Table shows a list of 25 networks predicted by IPA for the 552 RNA-Seq derived genes selected as being differentially expressed in the SON between SL and control rats. S6 Fig and Fig 5 illustrate the molecular interactions suggested for two of these networks.

The network shown in Fig 5 is of particular interest because it links the CREB family of genes which had been proposed as transcription factors that regulate transcription of the Avp gene in the Avp-MCNs through its CRE sites in its promoter that is activated by cyclic AMP [7,57,58]. It is known that the SON responds to hyperosmotic stimuli by increasing cAMP [39]. While the increase in Avp gene expression was initially thought to be caused by CREB activation, a study which injected a highly selective, dominant-negative CREB inhibitor (A-CREB) into the SON showed that this inhibitor which could block the induced expression of c-fos in response to acute hyperosmotic simulation, but did not inhibit the simultaneous increase in Avp gene transcription under the same circumstances [59]. In this regard, it is interesting that both the RNA-Seq and the microarray data identify three other CREB family members, Creb3l1, Atf4, and Atf5 [60], as being present in the SON. Each of these genes interact in the network shown in Fig 5, and all of them increase their expression in response to SL treatment (Table 3 and Table 8). None of these transcription factors’ actions would have been inhibited by the dominant-negative inhibitor, A-CREB [60,61], therefore all of these genes are candidates to bind the CRE sites in the Avp promoter and thereby regulate Avp gene expression. In this regard, a recent paper provides compelling evidence that one of these genes, CREB3L1, does indeed activate Avp gene expression by binding to the CRE-like and Gbox sequences between -170 and -120 bp in the Avp promoter [62]. This -170 to -120 bp region in the Avp promoter was also identified as a key domain involved in the cell-specific expression of the Avp gene in the Avp MCNs in in vivo experiments [63]. It would be interesting to determine if the other two CREB family candidates shown in Fig 5, Atf 4 and Atf5, also bind to these CRE-like and Gbox sequences in the Avp promoter and will also regulate the Avp gene’s expression in vivo, possibly in a synergistic manner. It is also notable that several other genes in the network in Fig 5, e.g., Cebpg, Psph and Tef also increase and some such as Dsp, and Agap2 decrease in expression during SL, but no member of the Akt gene family (a hub gene in this network) underwent a significant change in gene expression in response to the osmotic (SL) stimulus.

That there was also no change in expression of the ligand-dependent nuclear receptor RNAs in the SON in response to SL (Table 5) was surprising. This is because this family of transcription factors (nuclear receptors) have been implicated in the regulation of the Oxt gene, and have been hypothesized to act via a hormone responsive element found in the Oxt gene promoter [7,64] The latter data raise the question whether the expectation that a transcription factor gene that regulates a target gene must increase in expression when the target gene increases in expression under a specific physiological condition (such as for Oxt gene expression and SL). It is possible that such co-regulation need not be an absolute criterion for identifying a transcription factor candidate for any given target. Moreover, many of the genes in Table 5 that undergo substantial changes in gene expression in response to SL may be unrelated to the regulation of the Oxt and Avp genes in the SON, and are likely to be related to the other profound physiological and biochemical adaptive changes occurring in the MCNs during SL.

Concluding Remarks

A fundamental question in molecular neuroscience is how specific neuronal phenotypes are formed and maintained in the central nervous system. It is generally recognized that the expression of specific genes during development and in the fully differentiated state are important determinents of phenotype. MCNs in the SON are highly specialized to release large amounts of Avp or Oxt into the blood stream and play critical roles in the regulation of body fluid homeostasis. The MCNs are osmosensory neurons and are excited by exposure to hypertonic solutions and inhibited by hypotonic solutions. In addition, the MCNs respond to systemic hypertonic and hypotonic stimulation with large changes in the expression of their Avp and Oxt genes, and microarray studies have shown that these osmotic perturbations also cause large changes in global gene expression in the HNS. In this paper, we examine gene expression in the rat SON under normosmotic (control) and chronic salt-loading SL) conditions by the first time using “new-generation”, RNA sequencing (RNA-Seq) methods. We reliably detect 9,709 genes as being present in the SON by RNA-Seq. We also compare the SON transcriptome resolved by RNA-Seq methods with the SON transcriptome determined by Affymetrix microarray methods in rats under the same osmotic conditions. We find that there are 6,466 genes present in the SON that are represented in both data sets, although 1,040 of the expressed genes were found only in the microarray data, and 2,762 of the expressed genes are selectively found in the RNA-Seq data and not the microarray data. We find 552 of these genes were changed in expression as a result of chronic SL. These genes reflect diverse functions, and 42 of these are involved in either transcriptional or translational processes. These data provide the research community a comprehensive view of the transcriptome in the SON under normosmotic conditions and the changes that occur in specific gene expression evoked by salt loading. The RNA-Seq database we present in this paper is a valuable resource for the further study of SON functions during osmoregulation, and should drive new research by identifying new molecular targets to be interrogated with regard to osmotic regulation in the SON. In this regard, our IPA analysis suggests that a key network involving several Creb family transcription factors, e.g, Creb3l1, Atf4, and Atf5 (see Fig 5) should be worthy of further study under SL conditions. This RNA-Seq database also provides a baseline for the determination of the molecular differences between the Oxt and Avp phenotypes. As noted above, the SON is heterogeneous with respect to its cellular phenotypes. Hence, it is not possible from these RNA-Seq data of the whole SON to computationally deconvolute the individual contributions of the Oxt and Avp MCNs. However, transgenic rats containing green fluorescent (eGFP) Avp MCNs [65] or red fluorescent (mRFP1) Oxt MCNs [66] now exist and methods have been developed to dissect single, fluorescent neurons [28,67] for high throughput RNA analyses [68,69,70]. Therefore, it should now be possible to do RNA-Seq studies of the specific Oxt and Avp MCN phenotypes in the SON and to identify the unique RNAs that are present and are changed by SL in the transcriptomes of each of these individual phenotypes.