Genetic variability in microbial eukaryotes reshapes marine biodiversity assessment in the age of amplicon sequencing

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 1 / 29 OPEN ACCESS Citation: Trubovitz S, Sandin MM, Caron DA (2025) Genetic variability in microbial eukaryotes reshapes marine biodiversity assessment in the age of amplicon sequencing. PLoS One 20(6): e0326053. https://doi. org/10.1371/journal.pone.0326053 Editor: Claudia Isabella Pogoreutz, Université de Perpignan: Universite de Perpignan Via Domitia, FRANCE Received: December 2, 2024 Accepted: May 22, 2025 Published: June 20, 2025 Peer Review History: PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal. pone.0326053 Copyright: © 2025 Trubovitz et al . This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, RESEARCH ARTICLE Genetic variability in microbial eukaryotes reshapes marine biodiversity assessment in the age of amplicon sequencing Sarah Trubovitz 1,2*, Miguel M. Sandin3, David A. Caron1 1 Department of Biological Sciences, University of Southern California, Los Angeles, California, United States of America, 2 Ocean Sciences Department, University of California-Santa Cruz, Santa Cruz, California, United States of America, 3 Institut de Biologia Evolutiva (CSIC-Universitat Pompeu Fabra), Barcelona, Spain * strubovi@ucsc.edu Abstract

Rapidly improving DNA sequencing technology has revolutionized our ability to efficiently survey the biodiversity of microbial life. We are now equipped to inves- tigate protistan richness and community dynamics on scales that would not have been imaginable with traditional observational methods. However, for most taxa the relationship between DNA sequences and morphologically-defined species is poorly understood, and morphology has remained the cornerstone of taxonomy for centuries. To better utilize the wealth of sequence data being collected, we must understand how it relates to entities such as individuals, populations, and species. Here we use a combined microscopy and sequencing approach to unveil the striking intragenomic and intraspecies genetic variation in one group of ecologically-important marine protists, the polycystine Radiolaria. Long-read 18S rRNA gene amplicon data from 173 isolated and morphologically-identified radiolarians showed that the vast majority (90%) yielded multiple sequence variants per specimen. Furthermore, every morphospecies analyzed displayed a range of different genetic signatures. Intraspe- cies genetic variability was expressed as specimens having different assemblages of ASVs, different dominant ASVs, or having no ASVs in common with other specimens of the same morphospecies. By integrating morphological and molecular information, we begin to parse the genetic richness of Radiolaria in ocean environments, as well as illuminate relationships between taxa, and their poorly-known life stages. Our

findings emphasize the need to account for protists’ taxon-specific sequence variabil- ity, particularly their intragenomic and intraspecies genetic variation, in interpreting metabarcoding diversity survey data.

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 2 / 29 Introduction DNA sequencing has provided a wealth of information that has been transforming our understanding of marine microbial eukaryote species richness, biogeography, community composition, and evolutionary relationships [1,2]. In the field of ecology, high-throughput sequencing has enabled the detection of microorganisms that are oth- erwise too rare, fragile, fastidious, or morphologically indistinct to be readily observed with classical approaches of microscopy and culture. However, at this new frontier of DNA sequencing comes new challenges regarding data interpretation and integration with previous knowledge [3]. Metabarcoding is a highly-efficient technique for measur- ing biodiversity in environmental samples, which involves amplifying and sequencing a specific region within a chosen marker gene that is common to all targeted taxa (typically the 18S rRNA gene for microbial eukaryotes, a.k.a. “protists”). In this way, metabarcoding reveals the genetic richness of a sample, but ultimately the taxonomic assignment and interpretation of those gene sequences is limited by the quality of curated sequence databases (such as the Protist Ribosomal Reference Database, PR2 [4]). Due to the laborious nature of linking morphological and molecular data, these databases so far encompass only a small fraction of protistan richness, result- ing in many metabarcoding sequences with dubious taxonomic identities and little consensus as to what magnitude of sequence variation should delineate meaningful

biological units. Furthermore, while morphology-based taxonomy remains the gold standard in protistology, it is known to be imperfect, especially for rare, poorly- studied, and morphologically-indistinct taxa [3]. This poses the question: how can we best utilize genetic sequence information that is often only loosely tethered (if at all) to the morphological concepts that have historically been considered the common denominator in biological sciences? Discrepancies between morphology-based and sequence-based taxonomic assess- ments for some protistan groups can be confounding, and for one ubiquitous group of marine planktonic protists, the polycystine Radiolaria [5], these differences are partic- ularly stark. More than a century of taxonomic inventorying and description of polycys- tine radiolarian specimens in marine sediments, plankton nets, and sediment traps has yielded approximately 445 formal morphological species [6] and perhaps up to 600– 800, counting those that await formal description: most of these in orders Nassellaria and Spumellaria [5]. By contrast, the Tara Oceans expedition revealed an astounding ~6,600 polycystine genetic operational taxonomic units (OTUs) during a sampling expedition that targeted the sunlit ocean alone, and the vast majority of those OTUs were assigned to order Collodaria (another important polycystine group) rather than Nassellaria or Spumellaria [7]. This underscores how protistan morphological species descriptions cannot necessarily be expected to align with other modes of taxonomic

delineation, including concepts based on genetic similarity [8,9,10]. There is thus an urgent need to probe this mismatch between morphological and molecular richness so we can improve interpretations of the sequence datasets in genetic diversity surveys. Beyond epitomizing the disconnect between protistan morphological and molec- ular diversity, polycystine radiolarians are an ideal study system because their and reproduction in any medium, provided the original author and source are credited. Data availability statement: Sequence data underlying the analyses in this article are available from Zenodo (DOI: 10.5281/zeno- do.15550382). The raw, unprocessed sequence data are also available for download from the NCBI Sequence Read Archive (SRA) with sample accession numbers SAMN48798827 – SAMN48798999 (under BioProject PRJNA1269742). Funding: This project was support- ed by the National Science Foundation Postdoctoral Research Fellowship in Biology Award #2109767 (to ST) and the Simons Collaboration on Ocean Processes and Ecology p49802 (to DAC). MMS was supported by a postdoctoral fellowship from the Beatriu de Pinós programme of the Government of Catalonia's Secretariat for Universities and Research of the Ministry of Economy and Knowledge (grant #2021BP00068). Competing interests: The authors have declared that no competing interests exist.

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 3 / 29 relatively large size (~20–250 μm for solitary taxa and typically mm–cm for colonial taxa) permits simultaneous morpho- logical and molecular analyses, and their diverse array of siliceous skeletal structures facilitates high-resolution morpho- logical groupings. Additionally, radiolarians play integral roles in marine food webs and biogeochemical cycles, which necessitates a better understanding of their biodiversity. They are abundant in ocean water columns from the surface to abyss, from the tropics to the poles [11], occupying many ecological niches as mixotrophs, herbivores, predators, and detritivores [5,12]. They have also been shown to contribute significantly to the biological ocean carbon pump [13,14,15,16,17,18,19] and the silicon cycle [19,20]. Despite their importance, only a small fraction of polycystine radi- olarian morphospecies have been linked to their genetic sequence information [21,22,23,24,25,26,27]. In addition, very few studies have explored radiolarians’ intragenomic variation [23,28] (which in this manuscript refers to gene sequence variation within the genome of an individual specimen), or potential cryptic species [29,30] (which refers to genetically distinct subgroups within a morphologically-defined species). Here we build on emerging genetic information for polycys- tine Radiolaria by combining morphological and molecular analysis of 173 isolated specimens collected from the North

Pacific Subtropical Gyre. We employ long-read 18S rDNA amplicon data to investigate the genetic variability within individ- ual specimens and morphospecies, explore phylogenetic relationships among taxa, and discuss the implications of these findings for deciphering protistan diversity datasets in the age of amplicon sequencing. Results and discussion The dataset Our dataset included 18S rDNA long-read amplicon sequence variants (ASVs) from 173 polycystine radiolarian specimens belonging to 30 morphologically classified taxa. DNA from each individual was extracted, amplified, and underwent circular-consensus sequencing with Pacific Biosciences technology while morphology-based tax- onomy was assigned via microscopy. After bioinformatic processing, a total of 1,224 unique ASVs met criteria for analysis: a minimum of three reads in at least one specimen and taxonomically assigned to the same family as the morphologically-identified specimen they came from (see “Materials and methods”). The specimens in our dataset represented all three polycystine radiolarian orders: Collodaria (95 specimens), Nassellaria (36 specimens), and Spumellaria (42 specimens). Out of these, the Collodaria displayed the highest genetic richness (836 ASVs) despite including the lowest number of morphotaxa (9). Conversely, Nassellaria specimens had the lowest genetic richness (187 ASVs) but highest number of morphotaxa (11). Spumellaria specimens were comprised of 201 total ASVs from 10 morphotaxa. Taxonomic names of all specimens are listed in “Materials and methods,” photographic documenta-

tion is provided in S1 Fig, and ASV occurrence charts in S1 Data. Individual radiolarians contained multiple ASVs After stringent taxonomy and abundance-based vetting of ASV results (see “Materials and methods”), we found that 156 out of 173 radiolarian specimens (90%) exhibited multiple ASVs, with a median number of 12 ASVs per specimen (Fig 1a). This indicates that intragenomic variation was pervasive among the specimens we analyzed; however, the number of ASVs per specimen varied by taxonomic order. Collodaria specimens showed the greatest genetic richness, with a median of 32 ASVs per specimen and a relatively large range in values (1–141 ASVs per specimen) (Fig 1a). Nassellaria and Spumellaria had fewer ASVs per specimen, with median values of 3 and 5.5, respectively. Nassellarian intragenomic richness ranged from 1–57 ASVs and Spumellaria ranged from 1–40 ASVs per specimen (Fig 1a). The only family that did not yield multiple ASVs per specimen was Lophophaenidae (Nassellaria); Lophophaena hispida and Peromelissa thorac- ites each contained a single ASV (S1 Data). When the hypervariable V4 region was extracted from each long-read ASV (see “Materials and methods”), intrage- nomic richness of specimens in each taxonomic order was effectively reduced by half. There was a median number of 6 V4-region sequence variants per specimen, with Collodaria showing a range of 1–57 variants (median = 16), Nassellaria

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 4 / 29 Fig 1. Number of ASVs per specimen among three orders of polycystine Radiolaria. Each dot represents one specimen; panels show the distribu- tion of ASVs per specimen in Collodaria (pink), Nassellaria (brown), Spumellaria (green), and all taxa combined (gray). These include: the total number of ASVs that passed all vetting procedures described in “Materials and methods” (a); the subset of ASVs which constituted at least 0.1% of the number of reads per specimen (b); ASVs that made up at least 1% of reads per specimen (c); and ASVs comprising at least 10% read abundance per specimen (d). Panels (e–h) display corresponding intragenomic richness information for only the extracted V4 regions of the full-length ASVs, respectively. Box- plots denote the median, interquartile range, maximum, and minimum values (outliers excluded). All data underlying this figure is provided in S1 Data. https://doi.org/10.1371/journal.pone.0326053.g001

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 5 / 29 ranged from 1–21 (median = 2), and Spumellaria ranged from 1–37 sequence variants per specimen (median = 4) (Fig 1e). Therefore, while short-length V4 ASVs masked a significant portion of the genetic richness captured by full-length ASVs, intragenomic variation was still evident in each polycystine order albeit to a lesser degree. This result suggests that metabarcoding surveys utilizing relatively short-reads of hypervariable regions such as the V4 or V9 (e.g., [7]) are not able to detect the same level of genetic richness that can be revealed with long-read sequence data. The biological significance of protistan ASVs and their relationship to concepts such as specimens, populations, or species, remains a complex open question [3]. Previous work has suggested that ASVs almost certainly capture some intraspecies variation [3,10], but our results expand those findings by demonstrating that radiolarian ASVs can consti- tute a great deal of intragenomic variation as well. The fact that all three radiolarian orders demonstrated intragenomic variation indicates this phenomenon is not limited to species with a colonial life habit nor the presence of multiple nuclei, traits which frequently occur in Collodaria but have not been observed in any Nassellaria or Spumellaria [5,31]. Rather, it is probable that individual nuclei contain intragenomic variation for all groups included in this study (except evidently the

Lophophaenidae, but more research is needed). Furthermore, because the average specimen in our dataset yielded 12 unique ASVs, intragenomic variability is likely a significant source of the total genetic richness being measured in metabar- coding surveys of protistan diversity. The intragenomic variability we observed ranged widely among taxa, even those belonging to the same order. Therefore, while multiple ASVs could indicate multiple individuals (or multiple species) in some taxonomic groups, in others it appears to reflect the intragenomic variation within a single individual. This suggests that accurately interpreting the significance of ASV richness will require taxon-specific knowledge of genomic variability. Next, we examined ASV read abundances within specimens to assess the prominence of intragenomic variants. Thresholds were applied to determine the number of ASVs contributing at least 0.1%, 1%, or 10% of the total reads within each specimen (Fig 1b–d, respectively). These more abundant subsets of ASVs within each specimen presumably reflect their higher frequencies within the genome, potentially having a greater impact on the organism’s function and fitness compared to rarer ASVs [32]. Collodaria specimens yielded a median of 21 ASVs (range: 1–104) that comprised at least 0.1% relative read abundance (Fig 1b). At a minimum of 1% relative read abundance (Fig 1c), this decreased to a median of 5 ASVs per collodarian (range: 1–31), and at 10% (Fig 1d) collodarians exhibited between 1 and 4 ASVs with a median

value of 1 per specimen. Nassellaria and Spumellaria had fewer ASVs per specimen than Collodaria at every relative abundance threshold except for 10%, where Nassellaria exhibited the same median and range as Collodaria, and Spu- mellaria exceeded both groups by 1 ASV in terms of median and range (2 and 5 ASVs, respectively) (Fig 1d). The partitioning of reads among ASVs was investigated in greater detail by examining the relative read abundance of the top ten ranking ASVs per specimen in each order (Fig 2). Nassellaria specimens, on average, had a single, highly- dominant ASV, while a relatively low number of subdominant ASVs within their genomes comprised a minor percentage of reads. Among Nassellaria, the top ranking ASV made up a median of 92% of total reads within each specimen, versus 77% in Collodaria and only 60% in Spumellaria. Compared to the other orders, Spumellaria showed a relatively equitable distribution of reads among their dominant and top few subdominant ASVs. For example, the 2nd rank ASV in Spumel- laria comprised 23% of total reads per specimen on average, whereas 2nd rank ASVs in Collodaria and Nassellaria only made up 5% and 4% of reads per specimen, respectively. Among the extracted V4 region ASVs, the different partitioning patterns between Spumellaria and the other two orders were even more pronounced (S3 Fig). Collodaria and Nassellaria specimens had top ranking V4-ASVs make up 88% and 98% of their total reads (median values), respectively, whereas

Spumellaria specimens only had a median of 61% of total reads in their top-ranking V4-ASV. Based on these findings, we speculate that intragenomic variants in Spumellaria may tend to be more functionally important compared to those in other orders, even though Collodaria yielded the highest numbers of intragenomic variants overall. One possible explanation for the intragenomic variation observed in radiolarians is that they could be polyploid. Past imaging studies have suggested that some polycystines exhibit polyploidy (e.g., [33]); however, we are not aware of any definitive evidence that polyploidy is confirmed or common in polycystine radiolarians, and our dataset is not sufficient to

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 6 / 29 test this hypothesis. Polyploidy has been better documented in other members of Rhizaria, including Acantharia, Phaeo- daria, and Foraminifera [34,35], as well as other protistan groups including ciliates [36] and diatoms [37]. Another non-exclusive explanation for the intragenomic variation we observed could be the presence of high copy numbers of the 18S rRNA gene. Some polycystines (Collodaria) have been shown to have astoundingly high gene copy numbers per cell (estimated at 5,770 ± 1,960 SSU rDNA copies for solitary specimens and 37,474 ± 17,799 for cells iso- lated from colonies [38]), similar to groups such as dinoflagellates and ciliates [39,40,41], although it has been suggested that Nassellaria and Spumellaria have fewer copies [28]. Copy number has been shown to positively correlate with cell length across marine protist clades [38] as well as with intragenomic variation in ciliates [42]. Thus, copy number may have played a role in the high intragenomic richness we observed in Collodaria, and may also partly explain the relatively lower intragenomic richness observed in other polycystine groups that tend to have smaller cell sizes. Regardless of the underlying cause, differences in intragenomic richness and ASV read abundance structure across taxa have major implications for diversity interpretations in metabarcoding studies. High-throughput environmental DNA sequencing yields an immense quantity of data that is broadly representative of diversity in a sample. However, both

the richness and relative abundance of ASVs per specimen was not uniform across major taxonomic groups, so their metabarcoding results should not be interpreted uniformly to avoid misleading conclusions regarding diversity. Global scale genetic surveys of marine protistan species richness have revealed such great diversity of radiolarian molecular taxa [7] that they exceed the number of taxa in morphology-based surveys of similar scale by more than an order of Fig 2. Distribution of relative read abundances among the top 10 ranking ASVs per specimen. Dots represent relative ASV abundances (per- cent of total reads yielded by each specimen) at each designated rank (1st–10th) for all 173 specimens. Boxplots denote the median, interquartile range, maximum, and minimum values (outliers excluded) of ASV relative abundances for each order (Collodaria = pink; Nassellaria = brown; Spumel- laria = green). Data underlying this figure is in S1 Data. https://doi.org/10.1371/journal.pone.0326053.g002

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 7 / 29 magnitude [43]. The high degree of variability we observed in radiolarians’ intragenomic richness and ASV abundance structure may underlie some of the apparent ‘species’ richness and relative abundances observed in metabarcoding diver- sity surveys. The prevalence of this situation may extend to other protistan taxa as well. Previous work has documented intrage- nomic variation in Nassellaria and Spumellaria [23,28], and in polycystine radiolarians’ sister clades, the Acantharia [23] and planktonic Foraminifera [44,45]. There is also growing evidence that distally-related protistan clades exhibit intrage- nomic variation (e.g., [42,46,47]), which taken together may partly explain the staggering number of sequence variants observed in marine metabarcoding diversity surveys (i.e., [7]). Intragenomic variation has not been explored for the vast majority of taxa nor is it typically accounted for in richness interpretations of metabarcoding data [48]. It is thus imperative that we gain a better understanding of the magnitude and structure of intragenomic variation across groups, so that we can more accurately apply metabarcoding data to address important questions in ecology, including the unified under- standing and detection of biological units, like individuals, populations, and species. Specimens of the same morphospecies showed varying ASV assemblages Morphospecies of which we were able to collect multiple individual specimens were analyzed for intraspecies genetic

sequence variation we observed in C. huxleyi specimens. The life cycle of C. huxleyi is poorly understood, but each specimen we collected was colonial, with 10s–100s of skeletonized cells matching the formal species description (sensu Strelkov and Reshetnyak [49]). It was common, however, for colonies to contain cells of varying sizes, sometimes in the process of dividing (also observed by Strelkov and Reshetnyak; Fig 19 in [49]), so we cannot rule out the possibility that differing stages of the colonies or the cells within them could be related to the genetic variation we observed. Alternatively, this intraspecific variation could be due to random neutral mutation that was not related to life stage. Although there was a lack of distinct genetic subgroups in this population (i.e., all specimens shared one, usually dominant, ASV) the fact that some specimens had different dominant ASVs and most ASVs were shared by less than half of specimens, could suggest ongoing genetic drift potentially leading to speciation, including the possible emergence of cryptic species. Either way, the multitude of ASVs observed among the 51 C. huxleyi specimens raises questions as to how to interpret that sequence diversity with respect to the morphologically-defined taxon. The other morphospecies in this study also demonstrated intraspecies genetic variation, but they did so with differing patterns of ASV presence and abundance (Table 1; Figs 4–6). Out of the 14 morphospecies examined, only one (Ptero-

canium praetextum) yielded the same dominant ASV in all specimens, although only 3 specimens could be examined (Fig 5d). There were also two subdominant ASVs found in this species; one was shared by all specimens and the other shared

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 8 / 29 by two out of three specimens. P. praetextum thus showed a relatively low level of intraspecific genetic variability, though specimens did not have identical ASV compositions. By contrast, other morphospecies showed higher degrees of genetic variation than P. praetextum. In Procyttarium primordialis? (Fig 4e), some specimens yielded different dominant ASVs (dark blue shades in heatmap), although any ASV that was dominant in one specimen was at least present in the others (light blue shades in heatmap). Even greater genetic variation was observed in Collosphaera huxleyi (S2 Fig), Disolenia zanguebarica (Fig 4b), and Euchitonia elegans-furcata (Fig 6c). For these species, the dominant ASV in some specimens was entirely absent in others (white color in heatmaps). Most morphospecies in this study included some specimens that were not linked by any shared ASVs at all, appear- ing to constitute separate genetic subgroups. This was the case for Acrosphaera murrayanna (Fig 4a), Siphonosphaera socialis (Fig 4c), Thalassicolla nucleata (Fig 4e; but see caveats regarding Thalassicolla nucleata in subsequent section), Pterocorys zancleus (Fig 5a), Acanthodesmia vinculata (Fig 5b), Eucyrtidium hexagonatum (Fig 5c), Dictyocoryne pro- funda (Fig 6a), Dictyocoryne truncatum (Fig 6a), and Didymocyrtis tetrathalamus (Fig 6b). These morphospecies each contained at least one representative specimen, or group of specimens, that had an entirely different set of ASVs than

other specimens. These cases in which no ASV was shared across all representatives of a morphospecies could be indic- ative of cryptic species. Notably, the morphological subspecies of Didymocyrtis tetrathalamus (D. t. tetrathalamus and D. t. coronatus; Fig 6b) were not found to be genetically distinct, as they shared some ASVs. Therefore, while Didymocyrtis tetrathalamus contained multiple genetic subgroups, the subgroups did not appear to directly correspond to the morpho- logical distinctions defining the two subspecies. Instead, their morphological differences may be attributable to a gene other than 18S or potentially an epigenetic mechanism. Genetic subgroups have previously been reported in another spumellarian morphospecies, Spongotrochus glacialis, but in that case genetic differences corresponded with habitat par- titioning and a minor difference in skeletal morphology (spine length) [30]. These examples demonstrate that intraspecies genetic variation was the norm rather than the exception among Col- lodaria, Nassellaria, and Spumellaria. Moreover, a single dominant ASV was rarely shared among all representatives of a given morphospecies (only 1 out of 14 species) and most morphospecies (9 out of 14) included some specimens with no Table 1. Categories of intraspecies genetic variation observed in this study. Morphospecies (# specimens) The same ASV was dom- inant in all specimens Some or all specimens had different dominant ASVs, but at least 1 ASV was shared by all There was no ASV shared by every specimen of the morphospecies

Collodaria Acrosphaera murrayanna (6) X Collosphaera huxleyi (51) X Disolenia zanguebarica (2) X Procyttarium primordialis? (3) X Siphonosphera socialis (6) X Thalassicolla nucleata (10) X Nassellaria Acanthodesmia vinculata (2) X Eucyrtidium hexagonatum (2) X Pterocanium praetextum (3) X Pterocorys zancleus (19) X Spumellaria Dictyocoryne profunda (6) X Dictyocoryne truncatum (3) X Didymocyrtis tetrathalamus (10) X Euchitonia elegans-furcata (5) X https://doi.org/10.1371/journal.pone.0326053.t001

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 9 / 29 shared ASVs at all, suggesting separate genetic subgroups that may reflect cryptic species (Table 1). Polycystine mor- phospecies therefore cannot be easily delimited based on their 18S rDNA sequences alone, because each one contains uncharted genetic richness and an unknown number of genetic subgroups with different ASV compositions. One way to begin to untangle the degree of genetic differentiation in polycystine morphospecies is to compare the genetic distance within species to the genetic distance between species. Characterizing these genetic distance thresholds could provide an approach for estimating species richness in environmental metabarcoding studies, when morphological species counts are not possible. We quantified radiolarian intragenomic, intraspecies, and interspe- cies genetic variation by computing mean p-distances (calculated in MEGA 12 [50]). P-distance is the proportion of sites with non-matching bases across all positions of two aligned sequences. Mean p-distance provides a metric of the average degree of difference between all sequence pairs comprising a pre-defined group (e.g., between all ASVs within a specimen or a morphospecies) [50]. Both full-length ASV data and the extracted V4 hypervariable region sequence variants were used to compute intragenomic, intraspecies, and interspecies p-distances (Fig 7a and 7b, respectively). Fig 3. Intraspecies genetic variation in Collosphaera huxleyi. Each bubble denotes one ASV; the size of the bubble reflects the total read abun-

dance of that ASV across all C. huxleyi specimens collected in this study (natural log). The position of bubbles along the y-axis represents the percent of C. huxleyi specimens that contained a given ASV (with ≥3 reads). Only one ASV was shared by 100% of specimens (top left), and this ASV also had the greatest number of total reads. ASV bubbles are colored with semi-transparent blue fill; the less transparent the blue color appears, the more ASVs were shared by a similar percent of C. huxleyi specimens. The solid blue region at the bottom of the plot reflects the high number of ASVs present in only one specimen (330 ASVs) and their wide-ranging read abundances within those specimens (bubble sizes). Embedded photos show a typical live colony specimen (top) and skeletal remains of three cells from a colony (bottom; scale bar = 50μm). The data underlying this plot can be found in S1 Data, and visualized as a heatmap in S2 Fig. https://doi.org/10.1371/journal.pone.0326053.g003

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 10 / 29 Fig 4. Heatmaps illustrating intraspecies and intragenomic ASV patterns in Collodaria. ASV presence and relative read abundance are shown for each morphospecies or group of morphotaxa examined in this study: Acrosphaera murrayanna (a); Disolenia zanguebarica (b); Siphonosphaera socialis (c); Sphaerozoidae colonies (d); and solitary Collodaria (e). Morphological classifications of each specimen are shown on the y-axis of each heatmap

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 11 / 29 Results showed that there is indeed greater genetic difference between morphospecies than within them, but the typical p-distances separating intragenomic and intraspecific from interspecific genetic diversity varied widely between taxonomic groups (Fig 7). Among Collodaria, intragenomic variation was expressed by p-distances ranging from 0–0.176 (median: 0.009); Spumellaria had a comparable median intragenomic distance of 0.009 but lower range in values (0–0.037), while Nassellaria specimens showed the lowest median and range (median: 0.001; range: 0–0.032). Intraspecies genetic varia- tion was higher than intragenomic distances for all orders, with median p-distance values of 0.037 for Collodaria, 0.004 for Nassellaria, and 0.013 for Spumellaria. However, orders displayed very different magnitudes of intragenomic and intraspe- cies genetic variation, with collodarian intragenomic distances on-par with the intraspecies distances in Spumellaria and even surpassing intraspecies distances in Nassellaria. Similarly, interspecies genetic variation was higher than intraspe- cies variation for all orders, but p-distance values were substantially different between them (Fig 7). The median p- distance separating Collodaria morphospecies in this study was 0.203 (range: 0.083–0.222), 0.113 (range: 0.026–0.164) for Nassellaria morphospecies, and 0.027 (range: 0.011–0.141) for Spumellaria morphospecies. When these analyses

were performed using the extracted V4-region ASVs (Fig 7b), the same patterns held but p-distances were greater at each hierarchical level. This is presumably because the relatively short (~440 bp) V4 hypervariable region is expected to contain a higher proportion of site differences between ASVs at each hierarchical group level compared the full-length (~1450 bp) ASVs, which are comprised of both conserved and hypervariable regions. These results indicate that p-distances have promise for delimiting polycystine species-level taxa in metabarcoding studies, but appropriate p-distance thresholds vary considerably based on the length and characteristics of the target amplicon regions, and the order-level classification of the ASVs in question. Additional single-specimen sequencing efforts involving more Radiolaria from different taxonomic groups will improve our ability to define p-distance thresholds for estimating species richness with metabarcoding data. These findings highlight the challenge of detecting specific morphospecies of interest in DNA metabarcoding sur- veys using the existing structure of reference sequence databases (i.e., PR2 [4]). Reference sequences derived from isolated specimens of known protistan species are vital for interpreting genetic diversity data. But these individual refer- ence sequences will, in many cases, fail to capture the range of prevalent sequence variants produced by species in the environment. This can lead to overestimation of species richness and evenness in a sample when interpreting molecu-

lar data relative to morphospecies data. If we wish to detect a given morphospecies in a DNA survey, it would be more appropriate to organize reference sequence databases by groups of multiple sequence variants that co-occur within the same specimen, and within specimens of the same morphospecies, thereby explicitly accounting for intraspecies genetic variation. Conversely, an approach that is less rigidly bound by morphospecies concepts may ultimately provide us with a more natural characterization of protistan richness. Work on Foraminifera has shown that cryptic species often occupy narrower niches than the morphospecies to which they belong [51]. This means we risk overestimating species’ biogeo- graphic ranges and niche breadths if we fail to consider their genetic variation and look only through the lens of morphol- ogy. Ecological models predicting future geographic shifts in plankton communities could thus be improved by knowledge of intraspecific genetic variation. Such models utilize species’ presumed habitat preferences to predict geographic shifts in plankton communities (e.g., [52]), but if genetically distinct subgroups of some morphologically-described species have different habitat preferences, then our assumptions about their future geographic shifts may be inaccurate. Therefore, interpreting taxa as assemblages of multiple shared ASVs, likely with more than one possible assemblage per morpho- type, could allow the use of metabarcoding data to shed light on their ecological preferences and roles that might other-

wise remain obscured behind classical morphospecies concepts or ungrouped, discrete ASVs. with a random numeric identifier. The x-axis of each heatmap denotes the ASVs (represented by tick marks) present in each specimen of a morphos- pecies or group of morphotaxa. Relative read abundances of ASVs per specimen are indicated by the shade of blue (light = low abundance, dark = high abundance); ASVs not present in a particular specimen are indicated by white background. Representative micrographs of morphotaxa are shown on the right of each heatmap. All data underlying this figure can be found in S1 Data. Photos of every specimen can be found in S1 Fig. https://doi.org/10.1371/journal.pone.0326053.g004

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 12 / 29 Fig 5. Heatmaps illustrating intraspecies and intraspecimen ASV patterns in Nassellaria. ASV presence and relative read abundance are shown for each group of morphospecies and genera: Pterocorys (a), Acanthodesmidae (b); Eucyrtidium hexagonatum (c); Pterocanium praetextum (d); and other Nassellaria (e). Morphological classifications of each specimen are shown on the y-axis of each heatmap with a random numeric identifier. The juvenile specimen could not be identified based on morphology; instead, it was classified as family Acanthodesmidae based on its dominant ASV match to a reference sequence (99.9% confidence). The x-axis of each heatmap denotes the ASVs (represented by tick marks) present in each specimen of a morphospecies or group of morphotaxa. Relative read abundances of ASVs per specimen are indicated by the shade of blue (light = low abundance, dark = high abundance); ASVs not present in a particular specimen are indicated as white background. Representative micrographs of morphotaxa are shown on the right and below each heatmap. All data underlying this figure can be found in S1 Data. Photos of every specimen can be found in S1 Fig. https://doi.org/10.1371/journal.pone.0326053.g005

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 13 / 29 Fig 6. Heatmaps illustrating intraspecies and intraspecimen ASV patterns in Spumellaria. ASV presence and relative read abundance are shown for each group of morphospecies or morphotaxa: Dictyocoryne (a); Didymocyrtis (b); Euchitonia elegans-furcata (c); and other Spumellaria (d). Morpho- logical classifications of each specimen are shown on the y-axis of each heatmap with a random numeric identifier. Juvenile specimens were classified

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 14 / 29 Intraspecies genetic variation has clear implications for the study of modern and future protistan biodiversity and biogeography, but it can also reshape our understanding of their fossil record, including historical responses to climate change. In paleontological investigations genetic information is typically not available, and morphology is the primary data used to characterize protists’ evolutionary and ecological responses to past environmental shifts (e.g., [53,54,55]). If morphologically-defined species contain differing degrees of genetic variation, possibly including cryptic species, then the impact of some morphospecies extinction or origination events in the fossil record could have had broader ecologi- cal significance than previously believed. As intraspecies variation becomes better understood for each morpho- taxonomic group, we may be able to better estimate the impacts of macroevolutionary events throughout the group’s history. Insights gained from morpho-molecular relationships Long-read metabarcoding provides a unique opportunity to examine both biodiversity and phylogenetic relationships with the same dataset [56]. Most metabarcoding studies utilize short reads of ~300–500 base pairs, making them heav- ily reliant on matches to full-length reference sequences for interpretation of evolutionary relationships between taxa. Because the length of our reads spanned most of the 18S rRNA gene (~1450 base pairs), we were able to construct a

phylogeny from our data that directly informs on evolutionary relationships between ASVs and is decoupled from refer- ence sequence similarity (S4 Fig). In this way, we investigated the connection between morphological taxonomy, total ASV composition, and the phylogenetic relationships of specimens’ dominant ASVs, similar to the approach used by others [24,25,26]. Our results variously revealed that genetic information can in some cases: (1) support morphospecies concepts, (2) indicate taxonomic concepts in need of revision, or (3) potentially inform on radiolarians’ poorly-understood life stages. (1) Supporting morphospecies concepts. One example of morphospecies supported by genetic information was found in the Spumellaria genus, Dictyocoryne. Although D. profunda and D. truncatum have similar morphological characters (Fig 6a), our results showed that they are indeed genetically distinct. All of the dominant ASVs in D. truncatum specimens branched separately from any of the dominant ASVs in D. profunda specimens (S4 Fig). Furthermore, none of the dominant ASVs found in the D. truncatum specimens were present in any D. profunda specimens, and vice versa (Fig 6a). This finding suggests that ASV data can be used to successfully distinguish them in a metabarcoding survey. The specimens assigned to Dictyo- coryne sp. had only partially-developed skeletons (either due to breakage or incomplete growth), making it difficult to determine their species identity from morphology (S1 Fig). Based on their patterns of shared ASVs and phylogenetic

clustering, however, we were able to infer that Dictyocoryne sp. specimens #1–4 belong to D. profunda rather than D. truncatum (Fig 6a; S4 Fig). This exemplifies a case in which molecular data can be used to assist in species classifi- cation, even when morphological characters cannot be observed. There was one ASV present in very low abundance in Dictyocoryne-truncatum-03 (5 reads) that was also present as a subdominant ASV in three D. profunda specimens. However, due to the extremely low number of reads in Dictyocoryne-truncatum-03, we interpret this ASV as environ- mental contamination that was inadvertently sequenced (see “Materials and methods”), rather than an authentically shared ASV between the two Dictyocoryne species. as Spumellaria and Pyloniidae based on the identity of their matching reference sequences (>99% confidence). The x-axis of each heatmap denotes the ASVs (represented by tick marks) present in each specimen of a morphospecies or group of morphotaxa. Relative read abundances of ASVs per specimen are indicated by the shade of blue (light = low abundance, dark = high abundance); ASVs not present in a particular specimen are indicated as white background. Representative micrographs of morphotaxa are shown on the right of each heatmap. All data underlying this figure can be found in S1 Data. Photos of every specimen can be found in S1 Fig. https://doi.org/10.1371/journal.pone.0326053.g006

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 15 / 29 (2) Highlighting taxonomic concepts in need of revision. Results from specimens in the Sphaerozoidae (Collodaria) clearly indicated that a taxonomic revision is needed at the genus level. Each of the Sphaerozoidae genera, Collozoum (no spicules), Rhaphidozoum (simple spicules), and Sphaero- zoum (branching spicules) were found to be polyphyletic (Fig 8). Their dominant ASVs were scattered throughout the family Sphaerozoidae, indicating that spicule presence and shape may have evolved independently in multiple clades and should not be considered diagnostic characters at the genus level. These results align with the findings of Biard et al. Fig 7. Mean p-distances of ASVs within individual specimens and morphospecies, and between morphospecies of the same taxonomic order. Panel (a) shows the average intragenomic, intraspecies, and interspecies p-distances for each taxonomic order, based on full-length ASV compositions. Results of the same analyses are shown for the extracted V4-region sequences in panel (b). Only the ASVs derived from specimens morphologically- classified to the species level were used in this analyses. All data underlying this figure is given in S2 Data. https://doi.org/10.1371/journal.pone.0326053.g007

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 16 / 29 Fig 8. Phylogenetic context of ASVs generated during this study. PR2 reference sequences and clades are shown in black (includes a total of 258 polycystine sequences and 18 Acantharea outgroup sequences). The dominant ASVs from specimens in this study are shown in green text with a numeric identifier and the morphology-based taxonomic name. Families or clades in which any of our specimens’ dominant ASVs clustered are

PLOS One | https://doi.org/10.1371/journal.pone.0326053 June 20, 2025 17 / 29 [24], who also showed that Collozoum and Sphaerozoum were polyphyletic, but we speculate that Rhaphidozoum may be polyphyletic as well. ASVs from our Rhaphidozoum specimens did not group in the same clade as the Rhaphidozoum acuferum reference sequences (see Fig 8); however, since only two specimens were collected during this study and few reference sequences were available, more data from this genus is needed to confirm its polyphyly. Our results also suggest that morphospecies concepts of solitary Collodaria may be invalid and instead synonymous with various colonial taxa. Phylogenetic analysis indicated that all solitary Collodaria specimens in this study (morpholog- ically assigned to Thalassicolla nucleata and Procyttarium primordialis?) either shared their dominant ASV with a colonial specimen, or their closest phylogenetic relative was a colonial specimen (Fig 8). One ASV (coll-002) was dominant in 9 specimens, morphologically-classified as Sphaerozoum sp. (colonial), Procyttarium primordialis? (solitary), or Thalassi- colla nucleata (solitary). Because this ASV comprised between 54% and 87% of reads in these specimens, it is not likely the result of environmental contamination or technical error and probably indicates an authentic shared ASV. Additionally, these specimens were picked from three different net tow samples on North Pacific Subtropical Gyre cruise HOT 339, discounting the possibility that their shared dominant ASV is attributable to physical contact during sample collection.

ASV coll-007 was also shared by both solitary and colonial taxa, and ASV coll-012 was dominant in several Thalassi- colla nucleata (solitary) specimens while also being most closely phylogenetically related to the dominant ASV in colonial Rhaphidozoum specimens (Fig 8). We additionally found that the dominant ASV in two Sphonosphaera socialis colonies (coll-003) was taxonomically assigned to the Thalassicolla melacapsa (solitary) reference sequence in PR2 with 99.5% confidence, suggesting these morphospecies may be conspecific. Furthermore, the ASVs mentioned above did not form a monophyletic clade, which is consistent with the findings of Biard et al. [24] that Thalassicollidae is polyphyletic (Fig 8). Together, these results support the hypothesis that solitary Collodaria are conspecific with colonial taxa [24,57,58], rather than a distinct taxonomic group, as they have previously been considered (e.g., [21,22]). Solitary Collodaria should thus no longer be considered valid species in biodiversity surveys, but instead a morphotype that occurs across multiple clades. (3) Informing on radiolarian life stages and development. Our results indicate that some juvenile polycystine specimens lacking skeletal development may be identifiable to the fam- ily level or lower based on their ASV composition and phylogenetic clustering. These specimens were devoid of diagnostic taxonomic characteristics under microscopy, so they may only be distinguishable and possible to classify with sequencing

methods. One putative juvenile specimen was genetically classified to the family Acanthodesmidae based on genetic sim- ilarity, but due to the poor representation of this family in reference sequence databases classification at lower taxonomic levels was not possible (S5 Fig; Fig 5e; S1 Data). Another putative juvenile was assigned to the family Pyloniidae, due to its genetic match to Larcopyle butschlii reference sequences (S5 Fig; Fig 6d). There were two other putative juvenile specimens for which all ASVs were assigned to Spumellaria, but no genetic matches were found below the order level (S1 Data). These juvenile specimens shared their dominant and many subdominant ASVs (19 ASVs shared in total) with a mature spumellarian specimen that was morphologically classified as lithellid sp. 3 (sensu Trubovitz et al. [55]) (Fig 6d; S5 Fig). Therefore, all three specimens likely belong to the same open-nomenclature taxon which, based on its phylogenetic clustering, may comprise a novel Spumellaria clade (S4 Fig). highlighted in green and collapsed for clarity. The number of dominant ASVs clustering in each group is shown with yellow bubbles and the number of specimens per clade, family, and dominant ASV is denoted with blue bubbles. Colored boxes headlined with ASV identifiers depict photographs of all specimens that shared the same dominant ASV, and correspond to the colored highlighting in the expanded Collodaria portion of the phylogenetic tree.