Oomycete metabarcoding reveals the presence of Lagenidium spp. in phytotelmata

Oomycete cultures, cox1 gene sequencing, and genus-specific primer design

The Lagenidium giganteum strain ARSEF 373 was accessed from the USDA Agricultural Research Service Collection of Entomopathogenic Fungal Cultures (ARSEF, Ithaca, NY) and was grown in a defined Peptone-Yeast-Glucose (PYG) media supplemented with 2mM CaCl₂, 2mM MgCl₂ and 1ml/L soybean oil (Kerwin & Petersen, 1997). Axenic cultures were processed for genomic DNA extraction using the Qiagen DNeasy minikit, as previously described (Olivera et al., 2016; Quiroz Velasquez et al., 2014). The genomic DNA preparations were used as templates in Polymerase Chain Reactions (PCR) in combination with the oomycete-specific cox1 primers OomCoxI-Levup (5′-TCAWCWMGATGGCTTTTTTCAAC-3′) and OomCoxI-Levlo (5′-CYTCHGGRTGWCCRAAAAACCAAA-3′). These primers were designed to overlap the standard cox1 DNA barcode used in other groups and recommended by the Consortium for the Barcode of Life (CBOL) initiative (Robideau et al., 2011b). PCR conditions corresponded to the following pattern repeated for 30 cycles: 95 °C for 30 s, 50 °C for 30 s, and 72 °C for 1 min. The resulting products were purified using the QIAquick PCR purification Kit (Qiagen, USA) and sequenced commercially using traditional Sanger technology (Macrogen USA). The generated sequences were aligned with homologous oomycete sequences obtained from the Barcode of Life Data System (BOLD) database of cox1 genes (Ratnasingham & Hebert, 2007). Alignments were performed using ClustalX with default parameters (Larkin et al., 2007). The cox1 gene alignment was used to visually identify regions suitable for genus- or species-specific primer design. Alignments corresponding to selected locations were used as inputs for the construction of sequence logos using WebLogo, version 3 (Crooks et al., 2004).

Phytotelmata sampling and plant identification

Phytotelmata were sampled from ornamental plants on the Nova Southeastern University (NSU) main campus in Fort Lauderdale, FL, USA. Four plants were selected based on two criteria, including a visual, tentative taxonomic characterization of plants as bromeliads, and the observable presence of a large volume of water within the plants axils. The precise location of each plant was recorded using the Global Position System (GPS). Phytotelmata samples consisted of a 100 mL volume of water collected using sterile serological pipettes, and transferred in sterile 50 mL conical tubes. The water samples were inspected visually for the presence of macroscopic debris and invertebrates. In addition, leaf tissues (2 to 3 cm²) were also sampled for each plant, in an effort to associate phytotelmata samples with plant taxonomic classification. The leaf samples were grounded in liquid nitrogen and processed for DNA extraction using the Qiagen DNeasy Plant Mini kit (according to the manufacturer’s instructions). The plant genomic DNA preparations were used to PCR-amplify plant barcodes using primers designed for previously characterized loci, including the trnH-psbA spacer region (Kress & Erickson, 2007; Kress et al., 2005) and the internal transcribed spacer (ITS) region of nuclear rDNA (Cheng et al., 2016) traditionally used for a wide variety of land plants, as well as the trnC-petN spacer marker used more specifically for bromeliad barcoding (Versieux et al., 2012).

Phytotelmata microbiomes DNA extractions and cox1 barcode amplification

Phytotelmata samples were vacuum-filtered through 47 mm diameter, 0.45 µm pore size nitrocellulose filters (Millipore), as previously described (Mancera et al., 2012), and the microbial fauna retained on these filters was subjected to DNA extraction using the MoBio PowerWater DNA isolation kit (according to the manufacturer’s instructions). A similar workflow (vacuum filtration and DNA extraction) was used to process negative control water samples. These samples consisted of 100 mL of water collected at a drinking water fountain located on the NSU campus, as well as a 100 mL of seawater collected off the coast of Hollywood Beach, FL, USA. The resulting metagenomic DNA preparations obtained from phytotelmata and negative controls samples were initially PCR amplified using the oomycete-specific cox1 primers OomCoxI-Levup and OomCoxI-Levlo and the reaction parameters described above. Products of these PCR reactions were visualized on agarose gels. Subsequently, aliquots (one µl, non purified) corresponding to the products obtained using the OomCoxI-Lev primer set were used as templates for a second round of amplification. These nested PCR reactions were performed using the Lagenidium-specific primers (LagCox F & R) under stringent conditions (30 cycles of the following pattern: 95 °C for 30 s, 68 °C for 30 s, and 72 °C for 1 min). Products of these PCR reactions (expected size: 525 bp, nested within the 700-bp products obtained using the OomCoxI-Lev set) were visualized on agarose gels, cloned using the Invitrogen TOPO technology and processed for commercial Sanger sequencing (Macrogen USA). Resulting sequences were evaluated through homology searches and phylogenetic analyses as described below.

Oomycete community assessment through cox1 metabarcoding

The phytotelmata cox1 libraries were prepared for single molecule real time (SMRT) sequencing using recommended protocols available from Pacific Biosciences (PacBio multiplexed SMRTbell libraries). The workflow included a two-step PCR amplification as previously published (Pootakham et al., 2017). First, fusion primers were custom designed by combining the OomCoxI-Levup and OomCoxI-Levlo primer sequences described above with the PacBio universal sequence. These primers were HPLC purified and further modified by the addition of a 5′ block (5′-NH4, C6) to ensure that carry-over amplicons from the first round of PCR were not ligated in the final libraries (Integrated DNA Technologies). The first PCR reaction used these primers to amplify cox1 fragments from all four phytotelmata metagenomic DNA preparations. Resulting products were gel-extracted and served as templates for the second PCR reactions. The second reaction used the PacBio Barcoded Universal Primers (BUP) so that unique combinations of (symmetrical) forward and reverse barcoded primers were associated with each phytotelmata samples. Products of the second amplification were purified (DCC, Zymo Research), and sent to the University of Florida Interdisciplinary Core for Biotechnology Research (ICBR) where amplicons were pooled in equimolar concentrations and further processed for library construction and SMRT sequencing. The PacBio raw reads were demultiplexed and assessed for quality at the ICBR. Quality control processing included eliminating poor quality sequences, sequences outside the expected amplification size (ca. 810 bp) and sequences that failed to include both flanking, symmetrical barcodes. High quality reads served as inputs for homology searches to assign taxonomic identification down to the genus level, using BLAST2GO (Conesa et al., 2005). Sequences homologous to Lagenidium spp. were further processed for thorough phylogenetic analyses. These sequences were trimmed to eliminate flanking 5′ and 3′ regions, and evaluated for redundancy (100% homology) and OTU clustering using the ElimDupes tool (http://www.hiv.lanl.gov/). Selected sequences were included in the alignment described below.

Phylogenetic analyses

The cox1 gene sequences generated from axenic cultures and environmental samples were aligned with homologous oomycete sequences using ClustalX (Larkin et al., 2007). Most orthologous sequences were downloaded from the BOLD database (Ratnasingham & Hebert, 2007) as described above. However, the alignment was also complemented with orthologous Lagenidium spp. sequences available from GenBank, including the cox1 sequenced fragments recently generated from Lagenidium spp. isolates collected on mammalian tissues (Spies et al., 2016). The complete cox1 alignment consisted of a 620-character dataset that contained 62 taxa. The position of the shorter, Sanger-based environmental sequences was inspected visually and confirmed based on the location of the Lagenidium-specific primers. The jModeltest program (Darriba et al., 2012) was used to identify the most appropriate maximum likelihood (ML) base substitution model for this dataset. The best-fit model consistently identified by all analyses was the Generalized Time Reversible model with a gamma distribution for variable sites, and an inferred proportion of invariants sites (GTR+G+I). ML analyses that incorporated the model and parameters calculated by jModeltest were performed using PhyML3.0 (Guindon et al., 2010). ML bootstrap analyses were conducted using the same model and parameters in 1,000 replicates. The phylogenetic tree corresponding to the ML analyses was edited using FigTree v. 1.4.4.

Lagenidium giganteum cox1 gene sequence analysis

The cox1 fragment generated from the Lagenidium giganteum strain ARSEF373 was 683 bp long, and its sequence was deposited in the GenBank/EMBL/DDBJ databases under the accession number MN099105. Homology searches (not shown) demonstrated that the generated sequence was 100% identical to cox1 sequences reported from two other strains of L. giganteum (strains ATCC 52675, and CBS 58084, with cox1 sequences publicly accessible under the accession numbers KF923742 and HQ708210, respectively). Both strains ARSEF 373 and ATCC 52675 were originally isolated from mosquito larvae, according to culture collection records. Further comparisons (not shown) indicated that sequences from these mosquito-originating strains appeared divergent from the cox1 fragments sequences generated from multiple strains of L. giganteum f. caninum that have been reported as mammal pathogens, yet also retained the ability to infect mosquito in laboratory settings (Vilela et al., 2015). These results highlight the potential of molecular barcodes such as cox1 to distinguish between the known Lagenidium strains.

Unsurprisingly, the entomopathogenic L. giganteum cox1 sequences were also different from sequences characterizing more phylogenetically-distant oomycetes, including Lagenidium, Pythium and Phytophthora spp., as well as other Peronosporales. These differences provided a basis to develop Lagenidium giganteum-specific primers, and the location ultimately selected for primer design is illustrated in Fig. 1. The specificity of the designed primers relied especially on the reverse primer, that is located on a region that is immediately (40 bp) upstream the OomCoxI-Levlo primer (Fig. 1). This region was characterized by the presence of a 5′-ATCA-3′ motif that was showed to be prevalent in Lagenidium: alignments demonstrated that it was present on all the publicly available cox1 sequences (41 sequences total) obtained from L. giganteum (both mosquito and mammal strains) as well as L. humanum (Fig. 1). In contrast, the motif was not found in L. deciduum sequences (3 sequences), and was found only sporadically in Pythium and Phytophthora sequences (most notably in Py. helicandrum, Py. carolinianum, and some strains of P. ramorum, P. cactorum and P. infestans). As a result, the reverse Lagenidium-specific primer was designed to incorporate the reverse complement sequence 5′-TGAT-3′ at its 3′ end, and overlapped additional polymorphic sequences between Lagenidium and other Peronosporales. The primer sequences were finalized at 5′-ACTGGATCTCCTCCTCCTGAT-3′ for the reverse primer (LagCox R), and 5′-TAACGTGGTTGTAACTGCAC-3′ for the matching forward primer (LagCox F).

Environmental detection of Lagenidium spp. in phytotelmata using Sanger sequencing

A total of four plants were selected for analysis (Fig. 2), and information about these plants taxonomic identification through barcoding is available in File S1. The oomycete- and Lagenidium-specific cox1 primers were used in combination with metagenomic DNA preparations representative of the four plant phytotelmata (Fig. 2E). As illustrated in Fig. 2, the first round of amplification, using oomycete- specific cox1 primers, consistently produced detectable amplicons of the expected size (ca. 700 bp) for all plant-based water sources, but not the control water sources, suggesting the presence of oomycetes in the four sampled phytotelmata. Similarly, the nested PCR amplifications, using Lagenidium-specific primers (Fig. 1) and stringent PCR conditions, also produced fragments of the expected, 525 bp- size, leading to the production of twelve high-quality sequences (three per plants) publicly available in GenBank under the accession numbers MN099114–MN099125. Homology searches demonstrated that all twelve of these newly-obtained, environmental sequences were more similar to Lagenidium spp. cox1 sequences than other any oomycete barcodes (not shown). However, alignments also revealed that none of the environmental sequences were 100% identical to the previously published Lagenidium spp. cox1 sequences obtained from known strains maintained in axenic cultures (based on the 484 bp fragment length), suggesting a yet-unsampled diversity within the Lagenidium genus. Using a traditional 97% distance level to build Operational Taxonomic Unit (OTUs), the twelve Sanger-based sequences clustered in two distinct OTUs. The first OTU consisted of the Lagenidium humanum cox 1 barcode (accession number KC741445) clustered with the three sequences obtained from P3 (these three sequences were identical) and two identical sequences from the P1 phytotelma. All other environmental sequences (three identical sequences from the P4 phytotelma, as well as one unique sequence from P1, and three unique sequences from P2) clustered in a second OTU that included all known cox1 sequences from L. giganteum, including the L. giganteum f. caninum cox1 barcodes. These preliminary findings strongly suggested that all environmental sequences corresponded to Lagenidium spp. cox1 genes, and that the mosquito pathogen Lagenidium giganteum is present in phytotelmata. Relative abundance calculations indicated that L. giganteum (7 reads out of 12, or 58%) appears more frequently than L. humanum-like isolates (5 reads out of 12, or 42%). The sampled sequences, albeit limited in number, also validated the newly designed primers as specific for the genus Lagenidium. All sequences were incorporated in the phylogenetic analyses described below, in an effort to more precisely determine their taxonomic nature.

Assessment of Lagenidium spp. presence in phytotelmata microbiome using cox1 PacBio sequencing

A total of 23,857 reads were retained, demultiplexed and processed for bioinformatics analyses (File S1. Analyzed PacBio sequence datasets (available in the NCBI Sequence Read Archive data under accession numbers SRX6359420 –SRX6359423 as part of Bioproject PRJNA550619) included 7,852, 6,576, 5,151 and 4,278 reads for phytotelmata P1 to P4, respectively. Homology searches indicated that only a minority of these filtered reads (227 reads, or 0.9%) could not be assigned a taxonomic classification at the phylum/genus levels. Most sequences were classified into two major eukaryotic phyla, corresponding to animals and protists (Fig. 3). Animal sequences appeared to exclusively belong to insects and related taxa (Fig. 3), consistent with the hypothesis that phytotelmata are actively used environments for a specialized fauna of invertebrates. Protist sequences were further divided into oomycete and non-oomycete subgroups, and, as anticipated, oomycete sequences represented the majority of protist sequences in three sampled communities (Fig. 3). Oomycetes were found especially prevalent in phytotelmata P3 and P4, where they accounted for 79 and 90% of the sequences, respectively. Oomycetes represented 49% of the sequences in the P1 phytotelma, where the sequence distribution was characterized by a large proportion (40%) of invertebrate sequences (Fig. 3). These invertebrate sequences virtually all corresponded to a single OTU closely related to an unidentified Arachnida cox1 barcode (data not shown). In contrast to the P1, P3 and P4 samples, the P2 filtered reads contained a majority of non-oomycete sequences (Fig. 3), with an overrepresentation (82%) of OTUs homologous to the freshwater diatom genus Sellaphora (not shown). Oomycete sequences in P2 represented only 12% of the total sequences generated for this phytotelma (Fig. 3). These results pointed to the promises of using SMRT-based, long read cox1 sequences to assess the oomycete communities of selected environments but also suggested that the primer sequences, or the amplification conditions, used for these analyses may need to be refined in order to limit the production of amplicons from organisms that are phylogenetically close to oomycetes, such as diatoms. Overall, oomycete barcodes were detected in all phytotelmata, and sequence classifications at the genus level revealed a total of 10 oomycete genera, including Achlya, Aphanomyces, Halophytophthora, Haptoglossa, Lagenidium, Phytophthora, Phytopythium, Pythiogeton, Pythium and Saprolegnia. As illustrated in Fig. 3, Pythium, followed by Lagenidium, represented the most prevalent genera in the oomycete communities of all phytotelmata. Sequences identified as Pythium spp. corresponded to 90%, 93%, 39% and 97% of all oomycetes reads for plants P1–P4, respectively (Fig. 3) In agreement with the Sanger-based analyses, sequences homologous to Lagenidium spp. cox1 barcodes were detected in all samples. These sequences accounted for 7.2%, 1.7%, 59.8% and 0.3% of all oomycete reads, for phytotelmata P1 to P4, respectively, indicating that Lagenidium was present at low frequencies when compared to Pythium, except in the case of the P3 sample (Fig. 3). Also in agreement with the Sanger-based analyses, none of the reads identified as Lagenidium spp. were identical to the previously published L. humanum cox1 sequence fragment. However, a small number of reads (<1%) were shown to be 100% homologous to the mosquito pathogen L. giganteum cox 1 gene sequence (accession numbers HQ708210 and KF923742): 3 reads (out of 279) in the P1 sample and 1 read (out of 2,345) in the P3 dataset. OTU clustering at 100% distance level recognized identical reads within and between samples, and revealed that a single sequence was consistently the most predominant Lagenidium barcode across all four phytotelmata: this predominant sequence was represented by 103 reads out of 279 (37%) for P1, 3 reads out of 14 (21%) for P2, 1,215 reads out of 2,435 (50%) for P3 and 3 reads out of 13 (23%) for P4. Using a lower distance level for OTU clustering (97%), virtually all PacBio reads (>99%) clustered with these predominant sequences (not shown), and were associated with the L. humanum barcode. Finally, further sequence alignments compared reads obtained through Sanger vs. PacBio technologies. These comparative analyses showed that the overrepresented PacBio reads for P1-P4 were 100% identical to the sequences obtained using Sanger-based technologies for the P3 sample, highlighting the concordance between the two Lagenidium spp. barcode detections.

Phylogenetic analyses

The generation of novel Lagenidium-like cox1 sequences using both traditional and Next-Generation sequencing technologies prompted comprehensive phylogenetic analyses that incorporated these environmental barcodes within a robust alignment of sequences obtained from axenic cultures. The phylogram inferred from Maximum Likelihood analyses (ML) is presented in Fig. 4. The tree was rooted with representatives of the saprolegnian oomycete clade (Fig. 4), and focused on the peronosporalean clade, which includes the well-established Phytophthora and Pythium genera, as well as the more basal Albugo spp. (McCarthy & Fitzpatrick, 2017). The tree topology was very consistent with previously published oomycete phylogenies (Beakes, Glockling & Sekimoto, 2012; Lara & Belbahri, 2011; Spies et al., 2016), and depicted several Lagenidium species within a monophyletic clade and as sister taxon to a cluster containing a strongly supported monophyletic grouping of Phytophthora spp. and a paraphyletic assemblage of Pythium lineages (Fig. 4). The branch leading to Albugo spp. remained basal to this Phytophthora-Pythium-Lagenidium cluster. Although all Pythium species appeared monophyletic, deeper nodes, indicative of relationships between various Pythium spp., were characterized by weak statistical support. Similarly, poor bootstrap support prevented the confirmation of a recently proposed Lagenidium sensu stricto classification that regrouped L. giganteum, L. humanum and L. deciduum, and was inferred from a six-gene phylogeny reconstructions that included cox1 gene sequences (Spies et al., 2016). However, the present analysis confirmed the strongly supported, monophyletic association between L. giganteum and L. humanum (Fig. 4). All of the environmental sequences obtained from phytotelmata clustered within this Lagenidium clade, strongly validating the metagenomic approach, and the preliminary taxonomic identifications inferred from homology analyses. The environmental barcodes, independently from the amplification strategy and sequencing technology used to obtain them, segregated into two different groups: some sequences, including the most represented sequences generated using NGS technologies, appeared as sister taxa to L. humanum (99% bootstrap support), whereas another group of environmental sequences were strongly associated with the L. giganteum isolated from mosquito larvae (94% bootstrap support). Interestingly, no sequences appeared close to the L. giganteum f. caninum clade, or close to the more distant L. deciduum (Fig. 4), suggesting that, although the metabarcoding approach used in this study revealed a previously sub-sampled diversity within the genus Lagenidium, the sampling strategy may have biased the detection of Lagenidium spp. towards species that inhabit very specific ecological niches. The phylogenetic analyses clearly indicated that oomycetes such as L. giganteum and (possibly) L. humanum are present in phytotelmata, and that the metabarcoding approach described in this study provides a basis for the detection and isolation of novel Lagenidium strains independently of host-dependent baiting or occasional observations of infections.