Quantitative transcript and protein profiling

P deficiency is a critical driver of phytoplankton dynamics in the ocean, but the molecular level drivers of cellular responses to changes in P supply are not well described in many marine phytoplankton. Herein, deep sequencing of transcript tags and quantitative proteomic analysis were performed on pooled biological replicates of P-replete and P-deficient cultures of T. pseudonana CCMP 1335 (Figure S1), and coupled with physiological measurements to elucidate biochemical responses to P deficiency.

To examine coverage of the 11,242 predicted protein coding genes [30], the unique tags were unambiguously mapped with 100% identity to the T. pseudonana CCMP 1335 genome, and those falling within a gene model were tabulated. Tag-seq, detects both sense and antisense transcripts. Here, tags frequently mapped to both the forward and reverse strands of a given NlaIII site (data not shown). It is unclear whether this is a function of true antisense transcription or a methodological issue [31], but the observation of antisense transcripts and their frequency was similar in a tag-seq study of C. reinhardtii [29]. For the purposes of the analysis presented herein, antisense tags were excluded. The data validate transcription of >85% (9,572) of the modeled gene set (Table S1; Table S2). Past transcriptional profiling in this diatom using a tiling array detected transcription of 4,653 (43%) of the predicted protein coding genes [19]. Tag-seq compares favorably with RNA-seq [29], and is of equivalent or better coverage than other tiling array, DGE, or SAGE studies [19], [32]–[34] with non-model organisms, like phytoplankton, where the transcriptome size is uncertain.

Differential expression patterns were determined on the sequence tags using Analysis of Sequence Counts (ASC) [25]. A total of 1382 differentially regulated sequence tags were detected in this dataset with an estimated false discovery rate <0.05 using ASC. Analysis of the 50 most significantly upregulated sequence tags in this set, identified 27 that could be assigned to annotated gene models (Table S3). Of the remaining tags, 15 mapped to the genome, and 8 could not be assigned to the genome (Table S3). Three of the 8 unassigned tags mapped to ESTs in the JGI genome portal and thus, may represent tags that map to unpredicted splice sites (Table S3). The remaining 5 tags map to neither the genome nor the ESTs, and likely represent tags that map to splice sites in genes with no EST support, or possibly tags that cover SNPs and do not map because of the requirement for a 100% identity match. Tags that did not align with gene models but aligned to the genome often matched intergenic spaces proximal to gene models associated with P metabolism, suggesting this is a function of tags mapping to the 3′ untranslated regions (UTRs). For example, the percentage of tags that map to gene models increases by roughly 24% when each gene model is extended 200 bp into the 3′ UTR. In several cases, tags align to the genome in a relatively isolated region of intergenic space (Table S3). These data suggest the presence of novel P-regulated genes that were missed by the in silico gene modeling and highlight the value of surveying tags that fall outside of the modeled gene set.

Differential expression patterns were also determined on transcripts with one or more sequence tags mapping to that gene model in the genome. Using ASC [25], a total of 318 differentially regulated transcripts were detected (Figure 1A; Table S1) with an estimated false discovery rate <0.05. The number is lower than for the raw sequence tags because there are multiple tags mapping to secondary NlaIII sites in the same gene model. Although the tag-seq method is designed to sample the most 3′ NlaIII site [22], mapping of tags to secondary sites within a gene model is common [29], [32], and likely related to incomplete digestion during library construction.

Download:

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

Figure 1. Phosphorus deficiency signals in the transcriptome and proteome.

Global quantitative transcriptomic (A) and proteomic (B) analyses of T. pseudonana comparing P-replete (+P) and P-deficient (-P) treatments. Each data point represents a unique gene model (A) or protein (B), with those that are significantly differentially regulated noted in color. For the transcriptome data, tags mapped to a gene model were summed and then the total count was normalized to library size in tags per million (tpm). FDR is the false discovery rate p<0.05.

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

Quantitative global shotgun proteomic analyses were conducted in parallel to the transcriptome profiling, essentially providing independent replication of the regulation patterns observed in the transcriptome. Focusing on quantitative precision of the most abundant proteins, a total of 9512 unique peptides and 1264 unique proteins were detected with a false positive rate of 0.96% [35] (Table S4). Comparison of technical triplicate injections on the LC-MS demonstrated robust relative quantitation of the abundant proteins, with repeat injections showing similar spectral count values for each protein (Figure S2). The extent of protein coverage of gene products is similar to a study of the picoeukaryote Ostreococcus tauri [36] and a previous characterization of the T. pseudonana proteome [37], with the additional identification of certain P-responsive proteins (e.g. PID:23858) (Table S4). Changes in the relative abundance of proteins here were measured by spectral counting, where the number of spectra associated with a protein were counted and normalized. Using Fisher's Exact test (p<0.05), 136 differentially expressed proteins were identified (Figure 1B; Table S1). Of those, 79 were more abundant in the P-deficient treatment and 57 were more abundant in the P-replete treatment (Figure 1B; Table S1). On the order of 40 of the proteins were ≥2-fold more abundant in one of the conditions (Figure 1B). These observations show an increase in abundance of a number of transcripts and proteins that respond to an environmental stressor, and these targets may be useful for the development of biomarkers of diatom P deficiency in the oceans.

Coordination in the transcriptome and proteome reveals diverse aspects of the P stress response

A comparison of the transcriptome and proteome reveals a coordinated response for a number of targets, especially for those targets with a clear role in P metabolism (Figure 2; Table S5). In the combined dataset (where both the transcript and protein were identified for a given gene), 23 proteins were significantly more abundant in the P-deficient treatment, and 26 proteins were significantly less abundant in the P-deficient treatment (Table S5). Of the 23 proteins more abundant in the P-deficient treatment, no genes showed opposing patterns between the protein and the transcript. Coordination was also seen in the combined dataset for the 26 significantly downregulated proteins, with no transcripts showing opposing patterns (Table S5). Some transcripts were not significantly different between treatments, even though the corresponding protein was more or less abundant under P deficiency. This may be the result of post-transcriptional or post-translational regulation, or because they were expressed at low enough levels as to make it to difficult to resolve differences between samples. Of the few studies with model systems that have examined parallel responses of the proteome and the transcriptome in this way the extent of coordination between proteins and transcripts is similar. In yeast, 53% of proteins observed to be statistically more abundant had corresponding transcripts that were also statistically upregulated by rampamycin exposure [28]. In the green alga C. reinhardtii, 25% of transcripts responsive to copper exposure had a detectable protein with a coordinated pattern [29]. These studies have similar values to the diatom data examined here (60% for those that increased in abundance, 30% for those targets that were less abundant). With differences in the expression and turnover times for transcripts relative to proteins, and variability in post-transcriptional and post-translational regulation, there are many reasons why snapshot analyses of the transcriptome and proteome may not display similar patterns. Despite these possibilities, the patterns observed here for the diatom are tightly linked, particularly for the upregulated genes.

Download:

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

Figure 2. Comparison of transcript and protein signals.

Comparison of proteome and transcriptome changes in response to P deficiency. Fold-change presented as the log₂ of the ratio of deficient∶replete conditions. Unity lines are shown in grey solid (fold-change = 1), while a linear regression (log₂[proteins] = 0.49*log₂[transcripts]-0.25) of proteins that are >2-fold in abundance in either treatment against their corresponding transcripts is shown in yellow (r² = 0.53). The dashed line is the 1∶1 line denoting equal fold change between the deficient and replete conditions for the transcriptome and the proteome. Proteins and transcripts of interest that correspond to P-metabolism, glycolysis and ribosomes/translation are highlighted.

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

Analysis of the choreographed patterns, in combination with physiological data demonstrate that T. pseudonana has evolved a sophisticated, and multi-faceted, response to P deficiency involving five major strategies. These strategies include 1) change in cellular P allocation, 2) increased P transport, 3) a switch to utilization of DOP, 4) a remodeling of the cell surface, and 5) modulation of glycolysis and translation.

Changes in cellular phosphorus allocation.

Polyphosphate is comprised of variable length chains of phosphate, and it can serve as a storage compound that accumulates during luxury uptake in P-replete environments. Polyphosphate metabolism in yeast and other eukaryotes involves the vacuolar transporter chaperone (Vtc) 1–4 family [38]. In yeast, the Vtc4 protein interacts with the vacuole membrane and generates polyphosphate from ATP in a phosphotransfer reaction to form polyphosphate chains [39]. Two putative T. pseudonana Vtc genes were differentially regulated in this study. A Vtc4 polyphosphate polymerase homolog (PID: 43150) and an additional Vtc gene (PID: 38190) were significantly upregulated in the P-deficient transcriptome (Figure 3A). The VTC4 protein was also more abundant in the P-deficient proteome (Figure 3A). The other Vtc gene (PID: 38190) was not detected in the proteome, either because it is in low abundance, or because Vtc proteins are membrane associated and difficult to extract. Modulation of polyphosphate stores can occur in both prokaryotic and eukaryotic microbes experiencing general stress among other factors [40]. Although T. pseudonana is known to increase P allocation to polyphosphate in response to nitrogen limitation in chemostat cultures [10], polyphosphate in marine systems is typically thought to be the product of luxury uptake and storage of P in phytoplankton [8], not a stress response. Herein, the gene expression and protein abundance patterns are consistent with an increased allocation of P to polyphosphate detected in P-deficient T. pseudonana with solid state ³¹P NMR (Figure 3A). These data emphasize that not all diatom polyphosphate allocation is driven by luxury uptake in high P coastal systems, and thus supports the hypothesis that polyphosphate is present and its cycling could be important in low P systems [41], and in scenarios where diatoms are subject to environmental stressors. With these novel regulation data for a diatom polyphosphate polymerase, and their link to cellular polyphosphate changes, the transcript or protein may be used to track polyphosphate dynamics in field populations and could serve as a tool for better constraining patterns in P allocation, polyphosphate export, and concomitant apatite formation in marine sediments.

Download:

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

Figure 3. Transcript, protein, and physiological parameters linked to phosphorus deficiency.

Normalized transcript and protein abundance for significantly differentially regulated signatures and their associated physiological patterns, for polyphosphate metabolism (A), phosphate transport (B), alkaline phosphatase (C, D), phosphodiesterase (E), and sulfolipid synthesis (F) across P-replete (+P) and P-deficient (-P) conditions. Protein data are distinguished with a “p” next to the PID and by the hatched pattern. Polyphosphate abundances as measured by solid state ³¹P NMR (A) and enzyme activities (C, E), were assayed and are reported below each graph. Cell-associated alkaline phosphatase activity (green color) was detected using an enzyme labeled fluorescence substrate. The green fluorescence indicating enzyme activity is present in -P cells (panels 1, 2, 4, 5, and 6) and not present in +P cells (panel 3) (D). Chlorophyll autofluorescence (red) in also visible. Panels 4, 5, and 6 are a Z series through a labeled –P cell. The SQDG:PC ratio is reported from Van Mooy et al. (2009) from replete and P-deficient T. pseudonana cultures [51], which were grown similarly to those in this study. SQDG: sulphoquinovosyldiacylglyerol; PG: phosphatidylglycerol.

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

Conclusions

The diatom response to P deficiency involves a specific and multi-faceted remodeling of the transcriptome that is tightly coupled to the response of the proteome. There are also clear physiological changes in cellular P allocation patterns, enzyme activity, and lipid composition that are consistent with the transcriptome and proteome patterns. We identified differential regulation of a polyphosphate polymerase leading to a P deficiency-related increase in allocation of cellular P to polyphosphate. These data suggest that cycling and export of P as polyphosphate could be important in low or variable P systems where there is diatom P deficiency, in addition to high P environments where there is luxury uptake. This broader distribution of polyphosphate containing cells will influence the polyphosphate fueled formation and distribution of apatite containing sediments [8]. Although it is known that algae will induce phosphohydrolytic enzymes to acquire P from DOP in low P environments, to our knowledge this is the first work to demonstrate that diatoms have a cell-surface associated alkaline phosphatase, and may tightly couple hydrolysis and uptake. These results have implications for the competitive success of eukaryotic phytoplankton in low P systems and biogeochemical modeling efforts. Past work on DOP utilization has primarily focused on phosphomonoester hydrolysis and the enzyme alkaline phosphatase, but the strong induction of a diesterase suggests that recovery of phosphate from phosphodiester either from cell membranes, or exogenously, is a new critically important part of the diatom P stress response. Modeling the extent to which DOP supports primary production based on alkaline phosphatase activities alone, may thus be missing a large component of DOP hydrolysis. These observations, as well as the P regulation of a number of genes of unknown function, and the regulation of genes related to glycolysis and translation, provide new molecular-level insight into P stress responses in eukaryotic marine phytoplankton. In summary, these molecular data reveal that T. pseudonana has evolved a sophisticated response to P deficiency involving multiple biochemical strategies that are likely critical to its ability to rapidly respond to variations in environmental P availability, which should inform future field programs and biogeochemical modeling efforts.

Culture conditions and physiology measurements

T. pseudonana (Strain 1335 from the Provosoli-Guillard National Center for the Culture of Marine Phytoplankton (CCMP)) was grown in a modified f/2 medium [64] made from Sargasso Sea water. Macronutrients and vitamin B₁₂, biotin, and thiamine solutions were treated with prepared Chelex-100 resin to remove trace metal contaminants followed by trace-metal clean syringe sterilization [65] to yield final nutrient concentrations of 882 µM NaNO₃ and 106 µM Na₂SiO₃, and vitamin concentrations of 75 pM B₁₂, 400 pM biotin, and 60 nM thiamine. The Fe concentration was also modified from f/2 to 400 nM. All conditions were run in triplicate at 14°C, in constant light (120 µmol photons m⁻² s⁻¹). Cells were grown with f/2 phosphorus concentrations (P-replete; 36 µM PO₄) and with low phosphorus concentrations (P-deficient; 0.4 µM PO₄). Growth was monitored daily with cell counts. Duplicate P-replete treatments were pooled and harvested in mid log phase, and triplicate P-deficient treatments were pooled, and also quickly harvested onto 2 µm filters at the onset of P depletion (Figure S1A). All RNA samples were snap frozen in liquid nitrogen. Replicate P-deficient cultures were refed to 36 µM phosphate at the onset of P depletion, and subsequently resumed growth (Figure S1A). Additional growth studies were performed as described above substituting glycerophosphate and adenosine monophosphate at 36 µM, for the phosphate in the medium (Figure S1B).

Alkaline phosphatase activity was assayed on cells using a fluorometric method described elsewhere [12]. Labeling of alkaline phosphatase activity on whole cells was performed using enzyme labeled fluorescence (ELF) with previously reported methods [13]. Confocal microscopy of ELF-labeled cells was performed at the Marine Biological Laboratories using a Zeiss Inverted LSM 710, and a spectral detector to isolate the chlorophyll (680 nm) and ELF-97 product (535 nm) signals after Dyhrman and Palenik (1999). To assay diesterase activity cells were harvested via gentle filtration and resuspended into 2 mL of Buffer A [43] with 5.0 mM bis-4(nitrophenyl)phosphate. Absorbance was measured at 410 nm relative to substrate only and cell only controls on a spectrophotometer. Polyphosphate was detected using solid state ³¹P nuclear magnetic resonance spectroscopy on cells collected onto 0.2 µm polycarbonate filters and analyzed as described elsewhere [66], [67].

Transcriptome analysis

Transcriptome profiling was completed using a tag-seq profiling approach by Illumina Inc. Briefly, equal volumes from duplicate P-replete cultures, and triplicate P-deficient cultures were pooled prior to filtration. RNA was extracted from 5.40×10⁸ cells and 3.51×10⁸ cells, from the P-replete and P-deficient samples, respectively, using Qiagen's RNeasy Midi kit according to the manufacturer's instructions with the following exceptions; cells were lysed using 0.5 mm zirconia/silica beads (BioSpec) along with the lysis buffer solution and vortexed until the solution looked homogenous (2–3 min). The lysis solution was then put over Qiashredder columns (Qiagen) to aid in the removal of any large cell wall material that could clog the spin columns. To aid in the removal of DNA, an on-column DNase digestion was performed during the extraction using the RNase-free DNase Set (Qiagen) according to the manufacturer's instructions. Roughly 11 µg of P-replete and 8 µg of P-deficient total RNA was sent to Illumina and converted to cDNA. Illumina library construction included a polyA selection and digestion with NlaIII. This resulted in 21 bp sequence reads called tags that should map to the most 3′ NlaIII site in each transcript, but often map to alternative sites likely due to incomplete digestion during library construction.

Herein, tag-seq detected a total of 781,210 unique tags across both libraries generated from P-replete and P-deficient cultures (Table S1). The data are deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE28134. The 34.5 Mb T. pseudonana genome sequence [30] was used for tag annotation. Tag to gene mapping was performed through a pipeline designed by Genesifter Inc. In brief, the tag including the NlaIII site was mapped to the genome sequence data available at NCBI (AAFD00000000), requiring a 100% identity. Tags that mapped to more than one location were removed from the analysis. Of the roughly 12 million total sequence reads for each treatment, approximately 75% could be mapped to the genome (Table S1), a high number in part because of the presence of many unique tags with high counts. This percentage is similar to that observed for a DGE study in Chlamydomonas reinhardtii [29]. Approximately 10% of the unique tags could be mapped. This percentage of tag to genome mapping is somewhat lower than other studies with lower coverage [32]. For organisms with the vast majority of genes in their genome unmapped and uncharacterized, such as T. pseudonana, it is common to have relatively low percentages of tags mapping to the genome. Unique tags may not map because of single nucleotide polymorphisms (biological or introduced by sequencing error), because of incorrectly assigned intron/exon boundaries, and because tags map to 3′ UTRs that were not included in the predicted gene models among other possibilities.

Tag-seq, like Long-SAGE, detects both sense and antisense transcripts. Here, tags frequently mapped to both the forward and reverse strands of a given NlaIII site. For the differential expression analysis, all sense tags mapping within a given gene model were pooled, and antisense tags were excluded. Tag frequency was normalized to the total tags sequenced in a given library and the frequencies are reported in tags per million (tpm). Dynamic range was robust with gene model counts ranging from 1 to >5600 tags per million (tpm). Statistical significance in differential expression was determined by computing the posterior mean of the log ratio of proportions over the two conditions using an empirical Bayes method (ASC) [25]. In brief, the expectation of log fold change given the observed counts for each tag is estimated, and a standardized score is computed by dividing the estimated log ratio by a shrunken standard error estimate representing the biological variation. Transcripts, or tags, with a local false discovery rate (FDR)<0.05 are reported as differentially expressed (Table S2). The ASC statistic is designed for high density, high coverage datasets of transcriptional patterns generated with high throughput sequencing and is less biased by high tag counts than the Z score and Fisher test [25].

Proteome analysis

Biomass from frozen filter samples was scraped into tubes with a clean spatula then resuspended in 600 µL B-PER reagent (Thermo Scientific, Rockford, IL) supplemented with 0.5 mM EDTA and protease inhibitor 1mMphenylmethanesulfonylfluoride. Samples were vortexed for 1 minute and incubated at room temperature for 20 minutes then chilled on ice for 10 minutes. The cells were then sonicated with a microtip on ice for 1 minute at constant duty cycle, with the temperature not exceeding 28°C. Samples were centrifuged for 30 minutes at 14,100 RCF and 4°C, and supernatants were transferred and precipitated overnight in 50% acetone 50% methanol with 0.5 mM HCl at −20°C. Precipitated protein was collected by centrifugation at 14,100 RCF for 30 minutes at 4°C and dried by speed vacuum at room temperature. Protein was resuspended in B-Per/EDTA/PMSF extraction buffer, gently mixed and incubated for 30 minutes at room temperature and splits were combined. Aliquots were taken for protein determination by DC assay using bovine serum albumin as a protein standard (BioRad Inc., Hercules CA). Proteins were stored at −80°C until digestion.

Protein was then digested following the tube gel digestion procedure [68] with some modifications. Briefly, samples were immobilized in 15% acrylamide in pH 7.5 Tris buffer, incubated twice with 10% acetic acid and 50% ethanol for 20 min and 1 h, then with 10% acetic acid and 50% methanol for 2 h at room temperature and mixing at 350 RPM, decanting between. Gel samples were cut into ∼1 mm³ pieces and incubated twice with 50% acetonitrile 50% 25 mM ammonium bicarbonate solution for 1 h and again overnight, shaking at 350 RPM at 16°C and decanting between. Proteins immobilized within the gel were successively reduced with 10 mM dithiothreitol (DTT) at 56°C for 1 h, decanted, and alkyated with 30 mM iodoacetamide for 1 h at room temperature, washed in 25 mM ammonium bicarbonate for 20 minutes and twice with 100% acetonitrile for 10 minutes, and dried for 20 minutes by speed vacuum. Gel pieces were rehydrated and digested with a solution of 10 µg trypsin in 25 mM ammonium bicarbonate for 16 h at 37°C (1∶20 ratio trypsin to total protein, Promega Gold Mass Spectrometry Grade, Promega Inc., Madison WI). The peptides were extracted by successive additions of 50 percent acetonitrile (Fisher Optima) with 5% formic acid (Michrom Ultra Pure). The extracted peptides were combined and concentrated by speed vacuum for about three hours to less than 20 µL, diluted with 2 percent acetonitrile and 0.1 percent formic acid in water (Fisher Optima) and stored at −80°C.

The protein digestions were analyzed using a peptide Cap Trap in-line with a reversed phase Magic C18 AQ column (0.2×150 mm, 3 µm particle size, 200 Å pore size, Michrom Bioresources Inc. Auburn CA) on a Paradigm MS4 HPLC system (Michrom Bioresources Inc.) at a flow rate of 2 µl minute⁻¹. A LTQ linear ion trap mass spectrometer (Thermo Scientific Inc. San Jose CA) was used with an ADVANCE electrospray source (Michrom Bioresources Inc.). The chromatography consisted of a hyperbolic gradient from 5% buffer A to 95% buffer B for 300 minutes, where A was 0.1% formic acid (Michrom Ultra Pure) in water (Fisher Optima) and B was 0.1% formic acid in acetonitrile (Fisher Optima). The mass spectrometer was set to perform MS/MS on the top 7 ions using data-dependent settings and a dynamic exclusion window of 30 s. Ions were monitored over the range of 400–2000 m/z.

The mass spectra collected in this study were searched using SEQUEST (Bioworks version 3.3, Thermo Inc., San Jose CA) using an in silico tryptic peptide database assembled from the JGI modeled gene set (Thaps 3.0) (Thaps3_bd_unmapped_GeneModels_FilteredModels1_aa.fasta.gz and Thaps3_chromosomes_geneModels_FilteredModels2_aa.fasta.gz) with the addition of NCBI T. pseudonana chloroplast (EF067921) and mitochondria (DQ186202) genomes and a reversed ‘decoy’ version of these databases for false discovery rate analysis. SEQUEST parameters were set at 30% ions required per peptide, DCN of 0.1, Xcorr vs CS 1.9, 2.4, 2.9, and 1e-3 protein probability, with a false positive rate of 0.96% [35]. Changes in each protein's relative abundances across treatments were calculated using normalized spectral counts within Scaffold (Proteome Software V3.0; protein identification probabilities of 99.0% and requiring at least two tryptic peptide identifications per protein for 1264 protein identifications) (Table S4). Database search results were further processed using the PeptideProphet statistical model [69] within Scaffold 3.0 (Proteome Software Inc., Portland OR). Relative protein abundance was determined using Scaffold 3.0 for normalized spectral counting software operating on a 64bit Ubuntu Linux workstation. Spectral counts are normalized across samples in each experiment in Scaffold, including technical replicates, to allow comparison of relative protein abundance. Comparison of technical triplicate injections on the LC-MS demonstrates robust relative quantitation of the abundant proteins in the T. pseudonana proteome, with repeat injections showing similar spectral count values for each protein (Figure S2). Proteins discussed as ‘differentially expressed’ were determined by the Fisher exact test with p<0.05 [70] (Table S4).