Reduced representation sequencing accurately quantifies relative abundance and reveals population-level variation in Pseudo-nitzschia spp.

Certain species within the genus Pseudo-nitzschia are able to produce the neurotoxin domoic acid (DA), which can cause illness in humans, mass-mortality of marine animals, and closure of commercial and recreational shellfisheries during toxic events. Understanding and forecasting blooms of these harmful species is a primary management goal. However, accurately predicting the onset and severity of bloom events remains difficult, in part because the underlying drivers of bloom formation have not been fully resolved. Furthermore, Pseudo-nitzschia species often co-occur, and recent work suggests that the genetic composition of a Pseudo-nitzschia bloom may be a better predictor of toxicity than prevailing environmental conditions. We developed a novel next-generation sequencing assay using restriction site-associated DNA (2b-RAD) genotyping and applied it to mock Pseudo-nitzschia communities generated by mixing cultures of different species in known abundances. On average, 94% of the variance in observed species abundance was explained by the expected abundance. In addition, the false positive rate was low (0.45% on average) and unrelated to read depth, and false negatives were never observed. Application of this method to environmental DNA samples collected during natural Pseudo-nitzschia spp. bloom events in Southern California revealed that increases in DA were associated with increases in the relative abundance of P. australis. Although the absolute correlation across time-points was weak, an independent species fingerprinting assay (Automated Ribosomal Intergenic Spacer Analysis) supported this and identified other potentially toxic species. Finally, we assessed population-level genomic variation by mining SNPs from the environmental 2bRAD dataset. Consistent shifts in allele frequencies in P. pungens and P. subpacifica were detected between high and low DA years, suggesting that different intraspecific variants may be associated with prevailing environmental conditions or the presence of DA. Taken together, this method presents a potentially cost-effective and high-throughput approach for studies aiming to evaluate both population and species dynamics in mixed samples.