Transcriptome Profiling of a Toxic Dinoflagellate Reveals a Gene-Rich Protist and a Potential Impact on Gene Expression Due to Bacterial Presence

Gene Content and Gene Families

Previous MPSS analyses using well-annotated genomes have shown a strong correlation between the number of transcribed signatures and the total number of nuclear genes. In addition, these studies have demonstrated that as more libraries and conditions are examined, the number of unique signatures more closely represents the total number of predicted gene models in a genome. For example, in Arabidopsis, the number of annotated genes is 27,165 [22] and the number of unique MPSS signatures associated with protein coding regions gathered from 17 libraries is at least 29,569 [23]. Based on this correlation, we postulate that there are about 40,000 transcribed genes in Alexandrium, making it the most complex protist transcriptome yet described. It should be remembered, however, that although this number is relatively large compared to other free-living protists (e.g., 27,000 genes in the ciliate Tetrahymena thermophila [24] and 12,000 genes in the diatom Phaeodactylum tricornutum [25]), it does not account for the massive amount of nuclear DNA (ca. 150 Gb, estimated using pulse-field gel electrophoresis) in haploid Alexandrium cells. Clearly, gene number and genome size are uncoupled in these taxa. It is worth pointing out that in a recent genome size versus gene content regression study, dinoflagellates were predicted to contain 40,086 genes in the smallest genome and 92,013 genes in the largest [26].

This unusually high number of transcribed genes in Alexandrium is unlikely to represent unique functional categories; rather many may comprise large gene families that arose by extensive gene duplication events. To address this hypothesis in a conservative fashion, we first identified 4,341 expressed sequence tags (ESTs) from this strain that match perfectly and uniquely the identified set of reliable and significant MPSS signatures. Then we used KEGG Orthology (KO) [27] to functionally cluster these ESTs into families, resulting in the assignment of 1,020 KO entries to 2,169 ESTs (Table 2). The largest gene family comprises 31 members that encode peptidylprolyl isomerase (EC 5.2.1.8; cyclophilin). Subsequently, we counted the number of pairwise mismatches between signatures that correspond to ESTs clustered into the same families and ESTs belonging to different families. By comparing the numbers of pairwise mismatches between signatures from the two groups, we found that five mismatches can distinguish significantly between the two categories with p-value < 1E-10. Thus, using five mismatches as the maximum number of pairwise mismatches between signatures to obtain a rough estimate of the genome-wide distribution of gene families, we found 56 families with more than 100 members and the largest family contains 139 members (Figure 1). The largest family with members of known function contains 81 members and encodes pyruvate kinase (EC 2.7.1.40). The second largest family of known function encodes ribosomal protein L27a and contains 74 members. However, using KO-predicted families, we found cases where signatures within the same families shared low to zero identity. These cases are interpreted as duplicated genes with a relatively ancient common ancestor and the accumulation of mutations in the 3′ UTR has erased the phylogenetic signal in the signature sequences.

Examining the Alexandrium expression data drew our attention to several examples of different genes that belong to the same family and exhibit similar transcriptional profiles. For example, three S-adenosylmethionine synthetase (SAMS) genes were down-regulated in the bacterized culture. Similarly, three serine hydroxymethyltransferase (SHMT) genes were also down-regulated under this treatment. Genes encoding light-harvesting chlorophyll binding proteins followed the same pattern. In addition, four members of the ubiquitin family were up-regulated under nutrient limitation. To examine this association between gene family members and gene expression, we identified six signatures with a single mismatch between each pair with each of the six signatures having perfect matches to ESTs that encode the alpha subunit of the eukaryote translation elongation factor (EF-1α). The multiple sequence alignment (Figure 2A) of the signatures and their matching ESTs shows co-segregation of the mismatches among the signatures along with mismatches among the ESTs, suggesting these mismatches are not due to sequencing errors. Next, we found that the expression values of these six signatures are strongly correlated (Figure 2B) with a general pattern of up-regulation under nitrogen limitation and down-regulation in the presence of bacteria; i.e., when both are compared to the nutrient-replete culture. Therefore, the expression profiles among members of this (and perhaps many other) gene family are strongly correlated. This suggests that gene family expansion in Alexandrium may be a general mechanism used to enhance transcript abundance. Searching for similar patterns of co-regulation among family members, we found several families of different sizes (2, 4, and 8 family members; Table 3 ) that follow the same trend. In summary, our data indicate that dinoflagellate genomes contain large gene families with evidence for expression correlation among studied family members.

Bacterial Presence and Gene Expression

Although complex and multi-species bacterial assemblages have been shown to be associated with dinoflagellates both in extra- and intra-cellular environments [28], [29], taxa appear to be limited to the Cytophaga-Flavobacterium-Bacteroides (CFB) group and the α- and γ- classes of Proteobacteria. In this study, we did not attempt to identify the prokaryotes present in the bacterized Alexandrium culture. Previous studies have shown however that members of the genera Roseobacter (α-Proteobacteria) and Alteromonas (γ-Proteobacteria) are the dominant bacterial groups associated with Alexandrium sp. [30]. Here, we focused on the effect of the presence of bacteria in the culture on gene expression in the dinoflagellate. To identify transcriptionally regulated genes, we used the nutrient-replete culture as the control condition and identified signatures that were significantly up or down regulated in the other conditions using Fisher's test. At p-value < 1E-10, we found 1,124 signatures that were differentially expressed among the three (xenic, N-limited, P-limited) treatments compared to the control (Figure 3). By relaxing the p-value to a relatively more permissive threshold of 0.05 to detect even slight changes in expression among treatments, we identified ca. 11,000 differentially expressed signatures, indicating that about 29,000 signatures are consistently expressed with non-significant differences under the culture conditions used here. In dramatic contrast to a recent study which showed that about 6% of the expressed genes in rice are uniformly expressed, housekeeping genes [31], our results suggest that about 73% of the Alexandrium transcriptome comprises a “core” component and 27% comprises the regulated component, under differing cellular or environmental conditions. Of the 1,124 signatures, 307 (27%) were differentially expressed in the xenic culture, of which 119 and 188 were up- and down-regulated, respectively. Of these differentially regulated transcripts, two sets of genes stand out because they are collectively involved in the regulation of two important cellular processes, the methionine-homocysteine cycle and photosynthesis.

Cultures

Alexandrium tamarense strain CCMP1598 was used for all treatments except for the xenic culture, for which the bacterized clone, CCMP1493, was used. Strain CCMP1493 was first isolated from a germinated cyst from Daya Bay, east of Hong Kong (latitude +22°17′60.00″ and longitude +114°17′60.00″) and identified by Enrique Balech in 1991. It was deposited in the Center for Culture of Marine Phytoplankton by Donald M. Anderson in 1992. The growth rates were f/2 – 0.4 divisions per day, f/40 N, and f/40 P – 0.095 divisions per day, and 0.37 divisions per day in the xenic treatment. In each of these treatments, the in vivo fluorescence was monitored daily, and the culture was harvested when the division rates were consistent for several days.

Expressed Sequence Tag and 454 Transcript Sequencing

Total RNA was extracted from cultures of CCMP 1598 grown under replete (f/2), nitrogen-limited (f/40 N), and phosphorus-limited (f/40 P) conditions as described above, using the Nucleospin RNA II purification kit (Clontech Laboratories, Mountain View, CA, USA) according to the manufacturer's protocol. A start and a normalized directionally cloned (3′ NotI-5′EcoR1) cDNA library was constructed from the pooled RNA as previously described [55]. The complete set of existing EST clones derived from a previous study of Alexandrium tamarense CCMP1598 [21] was then used in a DNA hybridization protocol with the normalized library [56] to generate a subtracted cDNA library for single-pass 3′ EST sequencing. We generated a total of 11,171 ESTs using Sanger sequencing of the subtracted library which were processed as previously described [57]. The clustering, which relied on the 3′ UTR regions to facilitate accuracy, was done using UIcluster v3.0.5 [58]. This procedure resulted in a total non-redundant “unigene” set of 6,723 unique clusters. These data were combined with the existing unigenes described by Hackett et al. [21] and clustered using CAP3 [59] with a 95% cutoff identity between overlapping reads to avoid over-assembly that could mask biologically significant differences among closely related gene families. This second round of clustering resulted in a Sanger-based database of 12,329 unigenes from Alexandrium. We also generated EST data from Alexandrium using ‘454’ pyrosequencing. For this procedure, equimolar amounts of total RNA from each condition were pooled and cDNA was synthesized from 1 µg of total RNA using the Clontech Super SMART PCR cDNA synthesis kit following the manufacturer's instructions with the following modifications. Second-strand cDNA synthesis was done with a single round of primer extension using a 5′ trans-spliced leader primer conjugated to Clontech's primer IIA sequence to select for full-length dinoflagellate transcripts. All dinoflagellate transcripts contain an identical 5′ trans-spliced leader sequence on mature mRNAs [60]. The product of this single round of primer extension was purified using the Qiagen PCR purification kit to remove the spliced-leader primers and the cDNA was amplified by PCR using the Clontech primer IIA according to the Clontech cDNA synthesis protocol. A single microtitre plate of 454 Titanium sequencing was done at the Arizona Genome Institute (Tucson, AZ, USA) using 5 µg of amplified cDNA. These data were assembled using gsAssembler (Roche NimbleGen, Inc., Madison, WI, USA) into contigs representing Alexandrium cDNAs. The 12,329 unigenes generated using Sanger sequencing were co-assembled with the 454-derived contigs under Seqman (DNASTAR, Madison, WI, USA) using the default settings into a total of 35,431 dinoflagellate unigenes. The combined Sanger and 454 EST data were annotated using a best-hit approach against Pfam (version 23.0) [61] with blastx.

Massively Parallel Signature Sequencing (MPSS)

The same sources of mRNA used to construct the cDNA libraries were also used to generate the MPSS libraries to ensure comparability between the EST sequences and MPSS data. Additionally, mRNA was extracted from cultures of CCMP1493 grown under replete (f/2) for the xenic condition. The cDNAs were captured according to Illumina's protocols as described in Erdner and Anderson [19] and Brenner et al. [20]. Briefly, the cDNA was digested with DpnII and then amplified using PCR. Each cDNA was tagged by a 32-base synthetic oligonucleotide. The tagged cDNAs were then hybridized to their complementary 32-base tags that were covalently attached to microbeads. Each bead has only a single type of tag, but it is present in excess and generally, about 100,000 copies of a cDNA can be bound to a single bead. The result was a library of microbeads, in which each bead contained about 100,000 identical copies of a cDNA fragment that was derived from a particular mRNA. Libraries of approximately 2 million microbeads were loaded into flow cells for sequencing, which was performed simultaneously on each microbead in a cell. The result was a 21 bp signature sequence for every bead (hence every mRNA species) in the sample. The sequences from one or more flow cells (each containing a portion of the sample) were combined to form a set of 350,000 signatures. All of the signature sequences in a data set were then identified and compared to all other signature sequences, and identical sequences were grouped and counted.

Expression Data Analysis

Using blastn, MPSS signatures were matched to the assembled unigenes. To increase the sensitivity of the blastn search, the option of filtering low complexity regions was disabled, the word size was set to two nucleotides, and the expected e-value was relaxed to 1E+3. The search results were validated such that no more than three mismatches between a signature and a matching EST were allowed and a perfect match within the four nucleotides of the DpnII site (GATC) was necessary. For each signature, the matching ESTs were ordered by the identity scores rather than original e-value-based by blast. A matching EST with the maximum identity score and minimum e-value was designated as the most likely signature-matching EST. The annotations of the matching ESTs were directly transferred to the signatures. The unigene set generated by this study, including all previous EST data available from Alexandrium can be accessed from the public project web site: http://dbdata.rutgers.edu/alexbase. This web site also provides the MPSS expression data and matching ESTs. The combined Sanger/454 sequencing derived unigene set and the MPSS tag expression data over the different culture conditions are also available as supplementary files Figure S1 and S2, respectively.

Differential Expression Analyses

Signature frequencies were transformed to transcript per million (TPM) normalized values, where a signature-normalized value equals the signature frequency divided by the sum of the frequencies of all signatures in a library. Signatures with ambiguous nucleotides (i.e., other than A, C, T, and G) or repeats of sizes (string of the same nucleotide) >7 nucleotides were excluded. Additionally, signatures with frequencies less than 4 TPM under all conditions were also discarded. Pairwise Fisher's exact tests [62] were performed to determine the statistical significance of the differential expression patterns between the different treatments. Considering two libraries and of signatures with frequencies for signature is and respectively, then the matrix (i.e., the contingency table) for the Fisher's test was prepared as following,

Finally, the probability was adjusted using the BH method [63] to control the false discovery rate. To determine the differentially expressed signatures, we used p-value < 1E-10 as a consistent threshold between all pairwise comparisons.