An improved primer set and amplification protocol with increased specificity and sensitivity targeting the Symbiodinium ITS2 region

Sample collection, DNA isolation, ITS2 primer pairs, and experimental setup

To establish a robust ITS2 NGS-based amplification protocol for Symbiodinium, we assessed the effectiveness of three commonly used primer pairs (Table 1) on samples from different environments (Table 2). Further, we tested the robustness of the most effective primer pair, SYM_VAR, by documenting the effects of PCR-to-PCR variation, annealing temperature, cycle number, and polymerase system used with regards to: specificity, sensitivity, relative sub-genera amplification, and which sequences were returned. Chimeric sequence formation was additionally assessed as a function of PCR cycles.

Sample collection: In total, 27 samples were used in our analyses from four different environments: coral tissue (CO), coral surrounding water (CSW), coral reef surface water (SRF-CO), and open ocean (SRF-OO). Samples were collected either as part of the Tara Oceans or Tara Pacific expeditions. Field permits were obtained as part of the Tara Oceans and Tara Pacific expeditions and the approving bodies were the respective governments and/or ministries for the environment of the countries where the samples were collected. The authors of this work are part of the Tara Oceans and Tara Pacific expeditions, and thus, participated in sample collection. Tara Expeditions are global scientific voyages aimed at probing morphological and molecular diversity, evolution, and ecology of marine plankton (viruses, bacteria, archaea, protists, and planktonic metazoans) in the photic layer of the world oceans to research how they are impacted by Earth’s climate changes. Tara meteorological, oceanographic, biochemical, and plankton morphological data are archived in the comprehensive Sample Registry of the PANGAEA database (https://www.pangaea.de). The Sample Registry links physical samples to metadata about sampling and analysis methodology performed on each sample. Whilst the Tara Pacific registry is still under construction at the time of writing of this manuscript, the Tara Oceans registry can be accessed at https://doi.pangaea.de/10.1594/PANGAEA.842237.

Six CO samples from three genera, 2 Porites lobata (samples: CO-0000150 and CO-0000151), 2 Pocillopora meandrina (samples: CO-0000208 and CO-0000209), and 2 Millepora platyphylla (samples: CO-0000302 and CO-0000303), were collected using SCUBA off the coast of Panama in July 2016. Corals were stored in DNA/RNA Shield (Zymo Research, Irvine, CA, USA). Two CSW and one SRF-CO samples were also collected off the coast of Panama in July 2016. CSW samples were collected (15–17 L) from water surrounding two P. meandrina coral colonies using a boat-mounted vacuum pump. SRF-CO samples were collected (18 L) directly from the surface of the reef. Water was pre-filtered through a 20 μm metallic sieve, and then vacuum filtered on a 3 μm polycarbonate filter of 142 mm diameter for 15 mins. After filtration, filters were immediately flash-frozen in liquid nitrogen.

In addition, six SRF-OO samples were collected off the coast of Panama and were processed as described above for the CSW and SRF-CO samples. Twelve additional samples were collected as part of the Tara Oceans expedition from a range of ocean basins including the Pacific, North Atlantic, and Indian Ocean. These samples were filtered to a size fraction of 5–20 μm as described in section 6a of Pesant et al. (2015).

DNA isolation: DNA was extracted using the ZR-Duet DNA/RNA MiniPrep Plus (Zymo Research, Irvine, CA, USA) kit for all coral samples. Further, DNA was extracted from all water samples as described in Alberti et al. (2017) for the 5–20 μm and 3–20 μm fractions. In brief, the protocol was based on simultaneous extraction of DNA and RNA by cryogenic grinding of cryopreserved membrane filters, followed by nucleic acid extraction with NucleoSpin RNA kits (Macherey-Nagel, Düren, Germany) combined with DNA Elution buffer kit (Macherey-Nagel, Düren, Germany).

ITS2 primer pairs: Three pairs of published ITS2 primer pairs specifically designed to amplify Symbiodinium DNA were tested. These were: SYM_VAR_5.8S2/SYM_VAR_REV (Hume et al., 2013, 2015), ITSintfor2/ITS-reverse (Coleman, Suarez & Goff, 1994; LaJeunesse, 2002), and ITS-DINO/ITS2Rev2 (Pochon et al., 2001; Stat et al., 2009). Primer details and amplification protocols are given in Table 1. For ease of comprehension, each primer pair will be referred to as the SYM_VAR, ITSintfor2, or ITS-DINO primer pair.

Experimental setup: Sequencing for this study was conducted in two consecutive efforts. In the first part, aimed at directly comparing the three primer pairs, each primer pair was used to amplify each of the 27 samples. In the second part, a subset of samples was amplified using the SYM_VAR primer pair with different PCR protocols. We tested the effect of number of cycles, annealing temperature, and polymerase system on selectivity, sensitivity, and stability to optimise the reaction. Details of the samples as well as the associated number of cycles and annealing temperatures used can be found in Table 2 and Table S1. All PCRs were carried out in 25 μL reactions with the following conditions: 12.5 μL Phusion High-Fidelity PCR Master Mix 2X (ThermoFisher Scientific, Waltham, MA, USA) or equivalent amount of Advantage 2 PCR system (Takara Bio, Nojihigashi, Japan), 1 μL of forward and reverse primer at a concentration of 10 μM, 0.75 μL of DMSO, 1 μL of genomic DNA at a concentration of between 5 and 10 ng, and 8.75 μL of ddH₂O with a GeneAmp PCR system 9700 (Perkin Elmer, Waltham, MA, USA).

Library preparation, sequencing, and sequence quality control

For both sequencing efforts, all PCRs were conducted in triplicate before being pooled for sequencing. Pooled PCR products were purified using 0.8× AMPure XP beads (Beckmann Coulter Genomics, Brea, CA, USA), then aliquots of purified amplicons were run on an Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) using the DNA High Sensitivity LabChip kit to verify lengths and were quantified using a Qubit Fluorometer (ThermoFisher Scientific, Waltham, MA, USA).

One hundred ng aliquots of the amplicons were directly end-repaired, A-tailed, and ligated to Illumina adapters (Illumina, San Diego, CA, USA) on a Biomek FX Laboratory Automation Workstation (Beckman Coulter, Brea, CA, USA). Then, library amplification was performed using the KAPA Hifi HotStart NGS kit (KAPA Biosystems, Wilmington,USA). After library profile analysis by LabChip GX (Perkin Elmer, Waltham, MA, USA) and qPCR quantification using PicoGreen (ThermoFisher Scientific, Waltham, MA, USA) in 96-well plates, each library was sequenced using 250 bp read length chemistry in a paired-end flow cell on a MiSeq instrument (Illumina, San Diego, CA, USA).

During sequencing, Illumina’s real time analysis was run with default settings to remove clusters with the least reliable data. Output BCL files were converted to demultiplexed fastq files using Illumina’s bcl2fastq package and in-house filtering was applied to reads that passed the Illumina quality filters (Alberti et al., 2017). Read pairs that mapped onto run quality control (QC) sequences (Enterobacteria phage PhiX174 genome) were removed using the bowtie package.

For each of the sequenced samples analysed in this study the same sequence QC was performed. Briefly, mothur_1.39.5 (Schloss et al., 2009) was used to create contigs from paired forward and reverse demultiplexed fastq.gz files using the make.contigs command. The screen.seqs command was then used (maxambig = 0, maxhomop = 5) to discard sequences putatively generated from sequencing errors. The unique.seqs command was used to create a non-redundant collection of sequences represented by a .name and .fasta file. The remainder of the QC was performed using both the .name and .fasta files produced. The pcr.seqs command (pdiffs = 2, rdiffs = 2) was used to trim the primer sequence regions from the returned sequences. In addition, this command was used to discard sequences in which the specific primer pairs could not be found—indicative of poor sequencing quality—allowing a deviation of no more than two nucleotides differences in either of the forward or reverse primer sequences. Next, split.abund was run (cutoff = 2) to discard sequences that were not found more than two times in each sequenced sample—again, to reduce incorporating sequences with sequencing errors. Sequences were once again made non-redundant using the unique.seqs command. Finally, the summary.seqs command was used to calculate the size distribution of the remaining sequences. Sequences below the 2.5% and above the 97.5% size percentiles were removed using the screen.seqs minlength and maxlength parameters. A summary of the number of sequences retained at three stages of the QC process (pre-QC, post-QC, and Symbiodinium only sequences) are given in Table S1 and Fig. S1.

Sequence data analysis—primer pair comparisons

To assess which of the three primer pairs performed the best, each primer pair was assessed for specificity (i.e. preferential amplification of Symbiodinium DNA) and sensitivity (i.e. ability to amplify Symbiodinium from Symbiodinium-rare environments). Additionally, potential within-Symbiodinium taxonomic biases were assessed by comparing the proportions of sequences for each primer pair that belonged to each of the nine major lineages of Symbiodinium, clades A–I (Pochon & Gates, 2010). Our assumption was that if no or only small clade biases exist, the proportion of sequences from each clade should be similar between primers. Consequently, a large deviation in clade proportions by any one of the primer pairs was regarded as indicative of taxonomic bias.

The specific PCR conditions used for each primer pair are given in Table 1 with PCR reagents as detailed above. The sequences returned from each PCR were annotated using blastn and the NCBI ‘nt’ database according to their closest match to one of the following categories: Symbiodinium, dinoflagellate, stony coral, Hydrozoan, uncultured eukaryote, other, or ‘no match’. To do this, the .fasta file produced from the initial QC from each sample was run against NCBI’s ‘nt’ database with the max_target_seqs argument set to 1 and an output format string of ‘6 qseqid sseqid staxids sscinames sblastnames evalue’. For sequences to be categorised as Symbiodinium, an e-value >50 was required (representing approximately a 100% coverage match at 80–85% nt identity) in addition to the closest match being of Symbiodinium origin. Additionally, closest match subject sequences were screened for two sequences in particular (JN406302 and JN406301), which are mis-annotated as Symbiodinium (highly divergent from any other Symbiodinium sequences) in the ‘nt’ database. Thus, query sequences matching these sequences were annotated as ‘Unclutured eukaryote.’ Notably, before controlling for this, 58% of all ‘Symbiodinium’ sequences amplified by the ITSintfor2 and ITS-DINO primer pairs (0% for the SymVar primer pair) in the SRF-OO samples matched these sequences as their closest match (Fig. S2).

blastn was also used to associate Symbiodinium sequences to one of the nine clades. To do this, a .fasta file was created for each sample that contained all sequences that had previously been categorised as Symbiodinium. This file was used as an input for blastn with the max_target_seqs argument set to 1 and an output format string of ‘6 qseqid sseqid evalue pident gcovs’. A custom BLAST database was used that contained a single representative sequence for each of the nine clades (Data S1).

Sequence data analysis—optimising the SYM_VAR PCR protocol

Having determined SYM_VAR as the most effective primer pair for the amplification of the Symbiodinium ITS2 region, we undertook a second sequencing effort to optimise the PCR protocol. Specifically, we investigated: (1) how consistent the PCR was (i.e. PCR-to-PCR effect) by conducting PCRs with the exact same reaction conditions but sequenced on different sequencing runs, (2) whether increasing the annealing temperature could increase sensitivity in non-coral sample environments, and (3) whether a decrease in the number of PCR cycles would still allow for a satisfactory amplification of Symbiodinium whilst decreasing chimeric sequence return.

PCR conditions were the same as for the SYM_VAR_5.8S2/SYM_VAR_REV primer pair (Table 1), but with a 30 s denaturing step, following the manufacturer’s recommendation (Phusion High Fidelity PCR; ThermoFisher Scientific, Waltham, MA, USA). Taxonomic identities and clade breakdown (Symbiodinium sequences only) of sequences returned from each sample (annealing temperature and cycle number comparisons only) were determined in the same way as for the initial primer pair comparisons described above.

As a final optimisation of the SYM_VAR amplification protocol, we tested whether the polymerase enzyme system used could affect the outcome of the PCR amplification. For 11 samples (Table 2), we contrasted amplification with Phusion High-Fidelity PCR Master Mix (ThermoFisher Scientific, Waltham, MA, USA) and the Advantage 2 PCR Enzyme System (Takara Bio, Tokyo, Japan). PCR amplification with the Advantage 2 PCR Enzyme System were conducted as shown in Table 1, but with 30 cycles.

To further assess the potential effect of PCR protocol modifications, we additionally assessed whether PCR-to-PCR variation, or any of the protocol optimisation modifications would affect which sequences were amplified. To assess this, the top 20 most abundant sequences (by relative abundance) were identified and plotted by decreasing abundance for pairings of sequenced samples for which PCR conditions were identical, or only differed in one of the investigated parameters (i.e. cycle number, annealing temperature, or polymerase system). This analysis was conducted for two samples, AW-0000035 and IW-0000015 and only for Symbiodinium sequences contained in each. It should finally be noted that although we did not test for it here, variability between different PCR cyclers (e.g. in ramping rate and accuracy of temperature binning) also has the potential to cause variation in a given PCR.

An online version of the final optimised PCR protocol is available through the Doi: 10.17504/protocols.io.n8edhte.

Design of the SYM_VAR primer pair

The SYM_VAR primer pair were designed to minimise clade bias and coral-host complementarity. Originally, the reverse primer of the SYM_VAR primer pair that binds in the 28S region of the rDNA, SYM_VAR_REV, was designed to be complimentary to the forward primer SYM_VAR_FWD that binds in the 18S (Hume et al., 2013). This primer pair was designed to amplify both ITS regions as well as the 5.8S region of the rDNA. These primers were designed with the incorporation of degenerate bases based upon a multiple sequence alignment containing representatives of all Symbiodinium clades so that, where possible, the primer sequences had an equal number of nucleotide mis-matches to each of the clades. Later, to be able to amplify only the ITS2 region of the rDNA, the SYM_VAR_5.8S2 primer was designed (Hume et al., 2015). Compatibility to the SYM_VAR_REV primer and a minimisation of clade bias were maintained whilst complementarity to coral host ITS2 rDNA sequences was also taken into account by including representative sequences from P. lobata and P. lutea (ITS2 sequences from which had been amplified using alternative primers) to the multiple sequence alignment (Hume, 2013). It should be noted that this primer pair was successfully used to recover Symbiodinium DNA from Acropora downingi, Cyphastrea microphthalma, Favia pallida, Platygyra daedalea, Porites harrisoni, and P. lutea in Hume et al. (2015).

Sampling permits

The following permits were required for the collection of the samples analysed as part of this study: Panama scientific sampling: SE/AP-18-16; Panama CITES: SEX/A-72-16, no. 05567; Portugal scientific sampling: 080522/2009/Proc.o E.17.b; Maldives scientific sampling: FA-D2/33/2010/02; Chile scientific sampling: 13270/24/74.

Primer pair comparisons for amplification of Symbiodinium ITS2

In the CO samples, the SYM_VAR primer pair exclusively amplified Symbiodinium DNA. For the ITSintfor2 and ITS-DINO primer pairs Symbiodinium sequences represented 84% (SD 11%) and 79% (SD 12%) of total reads, respectively. Remaining sequences associated with Scleractinia, Hydrozoa, or resulted in no BLAST match (Fig. 1).

In the CSW water samples, the SYM_VAR primer pair maintained a high specificity for Symbiodinium (98%; SD 3%). In contrast, the alternative pairs produced 51% (SD 23%) and 42% (SD 24%) of Symbiodinium reads, respectively, with the remaining sequences characterised as dinoflagellate or falling into the uncultured eukaryote category (Fig. 1).

In the SRF-CO sample, Symbiodinium amplifications were considerably lower compared to the CO and CSW samples. SYM_VAR amplified 43% of Symbiodinium reads compared to 1% in both alternative primer pairs. The remaining sequences in the SYM_VAR pair were primarily of dinoflagellate origin. In the other two primer pairs, although a large proportion were also dinoflagellate in origin, sequences primarily fell into the uncultured eukaryote category (Fig. 1).

In the SRF-OO samples, the ITSintfor2 and ITS-DINO primer pairs were more similar in their amplifications compared to the SYM_VAR pair. In contrast to the other primer pairs, SYM_VAR amplified less ‘uncultured eukaryotic,’ ‘other,’ and ‘no match’ sequences, and more dinoflagellate and Symbiodinium sequences. Specifically, amplification of Symbiodinium was achieved in 12, 7, and 8 of the SRF-OO samples using the SYM_VAR, ITSintfor2, and ITS-DINO primer pairs, respectively. All of the five libraries that achieved a Symbiodinium amplification above 5% of total reads were amplified with the SYM_VAR primer pair. In total, the SYM_VAR primer pair achieved a Symbiodinium amplification of 7% (SD 15%) compared to <0.5% (SDs < 0.5%) for the alternative primer pairs (Fig. 1).

Across all samples in which all three primer pairs amplified >5% of Symbiodinium reads, the proportions of sub-genus taxa amplified (Symbiodinium clades) were comparable (Fig. 1).

Establishment of an optimal PCR protocol using the SYM_VAR primer pair

First, we assessed whether the same Symbiodinium ITS2 sequences were returned from PCRs performed on the same samples, with identical reaction conditions, but sequenced on different sequencing runs. For the two samples tested, the 20 most abundant sequences were present in both PCRs with relative abundances of each sequence being comparable (Fig. 2).

Next, we increased the annealing temperature from 56 to 59 °C to see whether a gain in Symbiodinium selectivity, and therefore sensitivity, could be achieved without yielding sub-genus taxonomic bias. In the CSW and SRF-CO samples, increasing annealing temperature increased the proportion of returned Symbiodinium sequences from 95% to 98% and from 28% to 47%, respectively. However, in both samples a small increase in the proportion of clade D versus clade C was observed, from 45 to 47% and from 57% to 66%, respectively. In the two SRF-OO samples tested, no Symbiodinium amplification occurred at either 56 or 59 °C (Fig. 3). The 20 most abundant sequences were found in both libraries for the CSW sample tested. However, two of the 20 sequences were missing from one of the libraries in the SRF-CO sample tested (the 10th and 16th most abundant sequences; Fig. 2).

Additionally, we tested whether Symbiodinium could be sufficiently amplified using a decreased number of PCR cycles to reduce chimeric sequence formation or taxonomic bias. To accomplish this, we compared the above-used protocol to one with only 25 cycles. In the CSW and SRF-CO samples there was a greater proportion of Symbiodinium amplified at 35 compared to 25 cycles (∼5% greater in the CSW samples and 15% in the SRF-CO samples). However, in both of the SRF-OO samples no Symbiodinium amplification was seen at either cycle conditions. Only non-Symbiodinium, dinoflagellate taxa sequences were amplified. Of the nine samples considered in the chimeric analysis (the three CSW and SRF-CO samples considered above, as well as six CO samples) greater number of chimeric sequences were putatively identified at 35 cycles versus 25 cycles in three of them with an increase of 1% in each case (Table 3). The proportions of sub-genus taxa amplified (Symbiodinium clades) were comparable (Fig. 3) and all but one (the 20th) of the most abundant sequences were amplified in common (Fig. 2).

Last, we also tested to see whether there would be any effect of using different polymerase systems. In the CO and SRF-OO samples tested, no difference in proportion of returned Symbiodinium sequences was observed between polymerase systems used (100% and 0% amplification for all libraries, respectively). However, in the CSW and SRF-CO samples we observed a relative decrease of Symbiodinium reads in the libraries that used the Advantage 2 system (Fig. 4). Despite, this relative decrease, the polymerase system used had minimal effect on which Symbiodinium sequences were returned (top 20 most abundant; Fig. 2).