Diverse molecular signatures for ribosomally ‘active’ Perkinsea in marine sediments

Processing of sequence data

Using 454 pyro-sequencing, we investigated the diversity of Perkinsea in a selection of European marine samples, sequencing the V4 region of the SSU rDNA [27] using both rDNA and rRNA as template. A similar DNA-based approach has been used to investigate freshwater Perkinsea [14]. We obtained sequence data from samples collected in four European coastal sites (Figure 1A), including sediment and multiple size filtrates from the sub-surface and the DCM water column samples (using a plankton net for the 2000-20 μm fraction and sequential filtration for 3–20 μm and 0.8-3 μm fractions). All V4 sequence reads (~380 bp) were assigned to eukaryotic taxonomic groups using a custom-built pipeline developed by the BioMarKs consortium (see [28]). This analysis identified 271 sequences preliminarily classified as Perkinsea.

Diversity within marine Perkinsea

To investigate the diversity and environmental distribution of the Perkinsea-like sequence tags, we conducted a phylogenetic analysis focusing on the V4 region and including the 150 sequence clusters identified (Figures 2 and 3). The regions flanking the variable V4 region are relatively conserved, while V4 stems and loops are variable [27, 30]. The phylogeny was derived from a masked alignment of 330 characters and included a mixture of sites with fast and slow patterns of variation. As our analysis was limited to the V4 region the deep and intermediate nodes of the phylogeny are poorly resolved so that the tree is only helpful for demonstrating the diversity of Perkinsea-like sequences and not the internal topology of the Perkinsea group, consistent with the aim of this study.

To identify a conservative picture of Perkinsea diversity we classified the sequence diversity into ‘cluster-groups’ on the basis of two restrictive criteria: 1) moderate topology support (>0.6/60%/60%) and 2) possession of two sequences from separate samples. Using this approach we identified 38 phylogenetic clusters labeled as cluster 1-38 on Figures 2 and 3 in addition to the morphologically characterised Perkinsus and Parvilucifera groups. 30 of these clusters represent previously undescribed marine diversity-groups. Additionally, 42 unique sequence clusters (28% - labeled with circles on the right column in Figures 2 and 3) were not grouped into ‘cluster-groups’ using our classification criteria.

The 30 new marine cluster groups show no more than 97% sequence identity between each group, suggesting they constitute taxonomically distinct groups. In contrast, the four described species that belong to Perkinsus spp. have highly similar SSU rDNA sequences (>98%, [31]). If the pattern of SSU rDNA variation in Perkinsus species is consistent across the wider diversity detected here, then these 30 cluster groups are likely to represent a diversity of forms with distinct biological and/or ecological traits, i.e. given the biology of the known Perkinsea species the diversity of sequences detected here putatively represent parasites infecting a range of host organisms.

Marine freshwater transitions

A growing body of literature has addressed the frequency in which protist groups have spread between marine and freshwater environments (e.g. [32, 33]) with varying perspectives on the number and relative ‘ease’ of these transitions dependant on the group studied and the criteria used for identifying these transitions [14, 32]. As described by Bråte et al. in 2010, marine-freshwater transitions are likely to have occurred during the diversification of the Perkinsea [14]. Our phylogenetic analyses identified the distribution of sequences recovered from both marine and freshwater environments demonstrating eight putative transitions on the phylogeny (five into freshwater and 3 into marine environments - Figures 2 and 3). However, we note that only three of these transitions are resolved in our V4 phylogenies with bootstrap support in excess of 50%. As such, additional sequencing from a range of environments combined with robust multi-gene phylogenetic analyses is required to characterise the frequency of freshwater-marine environmental transitions within the Perkinsea.

The majority of marine Perkinsea diversity is recovered from sediments

244 (92%) of the V4 sequences classified as Perkinsea were recovered from sediment. Moreover, 27% of the total sequencing effort were sampled from RNA derived cDNA libraries (Figure 1B,C and Table 1) suggesting that a significant proportion of the Perkinsea sequences were recovered from ribosomally active and intact cells. A large proportion of published environmental sequences are derived from DNA, a method that potentially detects dead organisms or extracellular DNA [34]. This is an issue arising from eDNA sequence surveys of sediment/soil environments. In contrast, extracellular RNA is thought to be less stable so that the use of rRNA can be useful for identifying ribosomally active microbes, inferring intact cells, but not distinguishing between active, senescent, dormant or encysted cells [35]. Therefore, we cannot exclude the possibility that the Perkinsea detected here are in a dormancy period or ‘dying’ whilst still maintaining transcription of a detectable RNA profile. However, these analyses identify a diverse range of ribosomally active Perkinsea in marine sediments, while in contrast recovering very little evidence of Perkinsea in the water column.

Sampling

Four European coastal stations were sampled (Figure 1A) as part of the work of the BioMarKs consortium (http://biomarks.scrol.fr): offshore Oslo (Norway, GPS position 59°15′N, 10°42′E), Naples (Italy, GPS position 40°48.5′N, 14°15′E), Blanes near Barcelona (Spain, GPS position 41°40′N, 2°48′E) and Roscoff (France, GPS position 48°46′N, 3°57′W). Each station was sampled over three depths (sediment, DCM and sub-surface) using the same sampling protocol (as described in [37] - Figure 1A). Environmental conditions and sampling area are described in [28]. Briefly, 30 to 50 litres of seawater were collected at the sub-surface and DCM either using a plankton net (for the fraction between 20–2,000 μm) or using Niskin bottles (for sampling of fractions less than 20 μm) coupled to a CTD sensor. Water samples were then size-fractioned using different pore size polycarbonate filters of 142 nm diameter (20 μm, between 3–20 μm and finally between 0.8-3 μm). Each filter was flash frozen and stored at −80°C for further analysis. Sediment samples were taken from a sediment core. Small aliquots of the surface sediment material (~1 cm³) were frozen and stored at −80°C for molecular analysis.

DNA/RNA extraction and 454 tag sequencing

For water column samples, DNA and RNA were extracted simultaneously using the NucleoSpin RNA L kit (Macherey-Nagel, Düren, Germany). For sediment samples, DNA and RNA were isolated using the PowerMax Soil DNA Isolation kit and the PowerSoil total RNA Isolation kit (MoBio, USA). DNA and RNA quality were confirmed using gel electrophoresis (1.5% agarose gels) and quantified using a NanoDrop ND-1000 Spectrophotometer. To avoid contamination by DNA in the RNA extractions, DNAse from the TurboDNA kit (Ambion, Carlsbad, CA, USA) was used to remove traces of DNA. Extracted RNA (100 ng) was reverse transcribed into cDNA using random primers and the Superscript III RT kit (Invitrogen, Carlsbad, CA, USA) following the protocol outlined by the manufacturer.

Universal eukaryotic primers TAReuk454FWD1 (5′-CCAGCASCYGCGGTAATTCC-3′) and TAReukREV3 (5′-ACTTTCGTTCTTGATYRA-3′) were used to sample the V4 region (~380 bp) of the SSU rDNA [27] using polymerase chain reaction (PCR) amplification. The primers were adapted for 454 sequencing with an A-adapter-tag forward and a B-adapter reverse as outlined in the 454 sequencing instructions. PCRs were performed in 25 μl mixtures of 1X Master Mix fusion High Fidelity DNA polymerase (Finnzymes, Thermo Scientific, Espoo, Finland), 0.35 μM of each primer, 3% dimethyl sulfoxide and 5 ng of template DNA or cDNA. PCR reactions consisted of an initial denaturation step at 98°C for 30s, followed by 10 cycles of: 10s at 98°C, 30s at 53°C and 30s at 72°C and then 15 cycles of 10s at 98°C, 30s at 48°C and 30s at 72°C. All PCR products were conducted in triplicate, checked using agarose gel electrophoresis (1.5% agarose gels), pooled and purified using NucleoSpinExtract II (Macherey-Nagel, Düren, Germany), eluted in 30 μl of water, and quantified using NanoDrop ND-1000 spectrophotometer. A final quantity of 200 ng of PCR product was then selected for 454 sequencing. Amplicon sequencing was carried out using a 454 GS FLX Titanium system (454 life sciences, Branford, USA) installed at Genoscope (http://www.genoscope.cns.fr), France.

Analysis of 454 reads of the V4 area SSU and phylogenetic analysis

Only reads with exact forward and reverse primer sequences and an estimated sequence error of ≤ 0.1% were retained for further analysis. Reads were assigned to taxonomic groups by co-clustering of sample sequences with those from a custom-built reference SSU rDNA database PR²[38] truncated to the V4 region. Reads were assigned to the Perkinsea when they were more similar to a reference Perkinsea sequence than to any other sequence in the PR² database, in terms of global alignment identity. This process identified 271 sequences tentatively classified as Perkinsea (each sequence has been labelled with a sequence number followed by the sequencing ID, see Additional file 5: Table S4 for details).

All existing Perkinsea SSU rDNA sequences (both environmental and from cultured organisms) plus a selection of 31 published sequences that encompass all the other major Alveolata lineages were recovered from the NCBI non-redundant nucleotide database and assembled into a reference dataset of 67 sequences (Additional file 2: Table S1 and Additional file 3: Table S2). The Perkinsea 454 V4 sequences were aligned to the previous reference dataset using Muscle [39], as implemented through the multiple alignment-editing program Seaview [39, 40]. The alignment was then improved manually, with particular attention to the V4 region. Ambiguously aligned characters were masked and excluded from the alignment prior to phylogenetic analysis. A preliminary tree was used to identify long-branch or highly novel sequences that could potentially represent chimerical sequences. Candidate chimerical sequences were investigated further by visual inspection of the alignment according to methods described by Berney and co-authors [29]. Of the 271 putative Perkinsea 454 tags identified from the Biomarks dataset, 265 marine V4 454 sequences branched with the Perkinsea and were retained for the final analyses (Table 1 and Additional file 5: Table S4). All 271 sequences are available in the European nucleotide archive (https://www.ebi.ac.uk/ena) under accession numbers PRJEB5698. We have included the six putative chimeras in the submission so these can be checked historically, these six sequences have the reference numbers 37–0005, 54–0139, 134–0005, 138–0140, 206–0147 and 268–0287.

Two datasets were then created, with different alignment masks: 1) a dataset encompassing the complete SSU rDNA sequence alignment and including a wide selection of Alveolata lineages (Alveolata SSU dataset composed of 1,437 positions and 377 sequences) and 2) a second dataset restricted to an alignment mask of the V4 region and focusing only on the Perkinsea phylotypes including sequences from Bråte et al.[14] (Perkinsea V4 dataset; 330 positions and 351 sequences). As the 454 sequences only encompassed the V4 region, for the first alignment, all missing positions in the Alveolata SSU dataset were encoded as gaps (consistent with the approach used in [14]). Prior to phylogenetic analyses we used the program Modelgenerator v0.85 [41] to determine the best model parameters for the two datasets. For the Alveolata SSU dataset a general time reversible model was selected with a discrete gamma distribution of the substitution rates (8 categories) and a proportion of invariable sites of 0.14 (GTR + Γ + I; gamma distribution shape parameter of 0.32). For the Perkinsea V4 dataset a GTR + Γ model was selected, with a gamma distribution shape parameter of 0.32 and 8 rate categories.

We then conducted Bayesian analyses using MrBayes v3.2.1 [42]. For both datasets we used the covarion parameter and a Γ rate correction with nst = 6 (equivalent to the GTR substitution model). The chains were run for 5,000,000 generations with two replicate tree searches both with 4 chains with a heat parameter of 2. Trees were sampled every 250 generations. In both analyses the MrBayes runs reached a stationary phase by 500 generation samples, and so the first 500 samples were discarded (as the burnin), and a consensus topology calculated from the remaining trees. For both analyses, the covarion model was compared to the non-covarion via Bayesian model comparison. This should be done using Bayes factors (the ratio of the respective marginal likelihoods for the two models) [43]. Unfortunately, the high dimensionality of parameter space makes the marginal likelihood term computationally intractable to evaluate directly. Therefore, the simplest, if somewhat imperfect, method of estimating the marginal likelihood is that of the modified [44] harmonic mean estimator [45] as implemented in the Trace package v1.4 [46] using 1,000 bootstrap pseudo-replicates. These analyses demonstrate that the use of covarion parameters produced an improved tree search (Additional file 6: Table S5).

For both datasets, support for the tree topology was evaluated by the bootstrap method and using the Bayesian posterior probabilities (PP) from the MrBayes runs [42]. Bootstrap support values (BV) were estimated using RAxML v7.0.3 [47], with 1,000 pseudo-replicates. For the Perkinsea V4 dataset, we also conducted a LogDet distance analysis [48] with 1,000 pseudo-replicates, as implemented in the Seaview [40] tree calculation module. This extra analysis was included to account for the possibility of compositional biases in the sequences [49]. We did not conduct LogDet analysis for the Alveolata SSU dataset because the large number of missing characters resulted in poor bootstrap results (which was not an issue for the likelihood and Bayesian analyses).