Conserved genetic markers reveal widespread diatom sexual reproduction in the global ocean

Domain: earth_science
Published: 2025-11-14
Editorial outcome: major_revision

Abstract
Sexual reproduction is an early universal characteristic of the eukaryotic life cycle, yet it is rarely observed in natural populations of micro-eukaryotes. Sex is particularly relevant for diatoms, a key group of marine and freshwater phytoplankton, where sexual reproduction counters a progressive cell size reduction due to cellular division. Here, we leveraged controlled sex transcriptome experiments of four diatom species to develop a robust method for in situ monitoring of sexual reproduction events. The resulting panel of conserved marker genes was validated for specificity and sensitivity using metatranscriptomic profiling of a natural estuarine community undergoing massive sexual reproduction of multiple species in response to increased salinity. Analysis of metatranscriptomic data linked with Metagenome-Assembled Genomes from the Tara Oceans expedition revealed widespread and coordinated expression of these markers across nine diatom genera, complemented by observations of sexual stages in automated imaging resources. Our results reveal that diatom sexual reproduction is more widespread in the global ocean than previously thought, encompassing both dominant bloom-forming species and rare taxa. Our panel of markers to detect sexual reproduction in natural environments paves the road to better understand the interplay between endogenous and environmental controls of this pivotal process, essential for the diatoms’ evolutionary success.

Introduction
Sexual reproduction is an early universal feature of the life cycle of eukaryotic organisms and is well-documented in animals and plants1. However, little is known about the environmental controls and phenology of sexual events in micro-eukaryotes2–4, including diatoms, an extraordinarily diverse and ubiquitous group of unicellular algae for which sexual reproduction has been extensively studied using laboratory experiments5,6. During the vegetative phase of their diplontic life cycle, mitotic divisions are accompanied by a gradual decrease in the average cell size of diatom populations due to constraints imposed by their rigid silica exoskeleton5. When they reach a critical sexual size threshold, cells become capable of producing haploid gametes that, upon fusion, form an expanding auxospore that restores the maximum cell size. Diatom reproductive strategies broadly fall into two categories: oogamous centric diatoms, producing eggs and flagellated sperm cells, are typically homothallic, while (an)isogamous pennate diatoms, producing non-flagellated gametes, are predominantly heterothallic (Fig. 1a).

As a result of their remarkable cell size reduction-restitution life cycle, diatom populations are expected to regularly undergo sexual reproduction to prevent critical cell miniaturization. Yet, direct field observations of this phenomenon have been sporadic7–13 because of the short duration—from hours to days—that sexual cells persist in the environment6, coupled with the challenge of identifying them through conventional light microscopy14. However, analysis of cell size distribution data15–18, population genetic studies19–23 and metatranscriptomic datasets24 have provided indirect evidence for the occurrence of sexual events.

In this study, we use an integrative approach to develop a set of sex marker genes which can detect diatom sexual reproduction within natural populations. Initially, we combine time-resolved transcriptomic data from controlled RNA-seq experiments to identify conserved and specific marker genes associated with diatom sexual reproduction. Subsequently, we assess their applicability and sensitivity through two metatranscriptomic studies. First, we investigate the potential for sexual reproduction in a microcosm experiment where we expose an estuarine diatom bloom community to sharp salinity increases. Next, we determine the occurrence of sexual events in the global ocean using metatranscriptomic data collected during the Tara Oceans expedition, linking them to individual Metagenome-Assembled Genomes (MAGs). The prevalence of sexual events is supported by automated imaging and further validated through metabarcoding and environmental data analysis. Our results suggest that the regulation of diatom sexual reproduction involves a complex interplay among local abundance, endogenous controls and environmental triggers, highlighting the need for further exploration of this key process at the base of diatom adaptation to rapidly changing aquatic environments.

Results
Conserved transcriptional activity during sexual reproduction
To identify conserved marker genes associated with sexual reproduction, we analyzed RNA-seq data from four diatom species representing different lifestyles and reproductive strategies: the centric diatom Skeletonema marinoi, the planktonic pennate diatom Pseudo-nitzschia multistriata, and the benthic pennate species Seminavis robusta and Cylindrotheca closterium (Fig. 1a, Tables S1, S2). To facilitate inter-species comparisons across experimental time series with varying sampling frequencies and durations of sexual reproduction, each time point was classified based on the predominant sexual stage: sex pheromone signaling, gametangia, gametes/zygotes, and auxospores. Multidimensional scaling of gene expression showed profound changes throughout the mating process (Fig. S1a, b). A total number of genes ranging from 1746 (S. marinoi) to 8557 (S. robusta) was differentially expressed (DE) during sexual reproduction compared to vegetative conditions (Fig. 1b). Upregulated genes included well-known diatom sex genes such as SIGs (Sex Induced Genes) and B-type cyclins25,26. S. marinoi exhibited a smaller number of downregulated genes, likely due to the absence of a pheromone-induced growth arrest25, which could lead to a general population-wide downregulation of vegetative genes, even in cells not actively engaging in sex.

Using gene families that define homology between species, we next assessed the conservation of gene expression levels across different diatom species. On average, 221±144 (st. dev.) gene families encoded DE genes at the same sexual stage across two different species (Fig. 1c). To assess the degree and direction (up- or downregulation) of this conserved response to sex, we evaluated the Pearson correlation of fold changes among shared gene families (Fig. 1c). Within the same species, correlations between adjacent time points were typically very high (ranging from 0.5 to 1) and always significantly different from 0 (Fig. S1c). Similarly, corresponding time points across species showed positive correlations (ranging from 0 to 0.5), even when comparing pennate and centric diatoms, which diverged circa 150 million years ago and employ completely different mating strategies27. Hence, a core transcriptional program is activated during sexual reproduction of both pennate and centric diatoms. Many genes involved in DNA replication and meiosis, as well as conserved signaling modules such as members of the Deltex E3 ubiquitin-protein ligase family, were upregulated in all four species during sexual reproduction (Figs. S2a, S3, S4). Conversely, shared downregulated gene families pointed towards an inhibition of primary metabolism and photosynthesis (Fig. S2b), consistent with previous reports28.

Identification of eight sex-specific transcriptomic markers
The identification of a shared transcriptional response across the four species prompted us to systematically search for conserved marker genes associated with sexual reproduction (Fig. 2). First, we selected diatom-specific gene families that were expressed in at least three out of four species and were consistently upregulated during sexual reproduction. This approach allowed us to identify markers specific to pennate diatoms, as well as those shared between centric and pennate diatoms. To ensure reliability, we prioritized markers that were sex-specific, i.e., showed negligible expression under non-sexual conditions. Therefore, gene families were ranked by their average fold change, a measure of the difference in expression between sexual and control conditions. The top-10 ranked families (M1-M10) exhibited more than 10-fold average upregulation during sexual reproduction, with the top-5 exceeding 30-fold upregulation (Fig. 3a). To verify marker specificity, we examined the S. robusta expression atlas, a dataset of 119 RNA-seq samples covering diverse conditions such as toxin exposure, nutrient depletion, environmental stress, diel cycles and bacterial interactions29,30. The top-4 ranked marker families (M1-M4) showed negligible expression across all 42 non-sexual conditions (Fig. 3b, Fig. S5), ensuring their suitability as robust markers. On the other hand, the other families showed similar expression levels under some non-sexual and sexual conditions (Fig. 3b). Hence, to minimize false positives when screening natural samples, only M1-M4 were selected as markers. Besides RNA-sequencing, we validated these markers using reverse transcription quantitative PCR (RT-qPCR) in C. closterium sexual crosses (Note S1, Fig. S6, Table S3).

In terms of their function, the M1 family encodes tubby-like proteins, suggesting a role in gene regulation, while M2 contains transmembrane lectin-like proteins. The M3 genes encode SIG7, a presumed distant homolog of the eukaryotic meiosis gene HOP225,31. Finally, M4 consists of proteins with unknown functions containing 6-10 transmembrane domains. While the conserved expression of M3 genes in gametangia aligns with their role during meiosis, M1, M2 and M4 were expressed at different sexual stages across different species, underscoring the need for further work to elucidate their specific functional roles (Fig. 3c).

Among the four markers identified, only M3 is present in both centric and pennate diatoms (Fig. 3d). M1, M2 and M4 are restricted to raphid pennate diatoms and are absent in centric and araphid pennate diatoms, an ancestral subgroup of pennate diatoms from which the raphid clade diverged approximately 120 million years ago27 (Fig. 3d). To identify additional markers for centric diatoms, we focused on genes associated with flagellum formation and function, as flagella are characteristic of sperm cells in centric diatoms. We gathered previously described proteins from the diatom flagellar body and the mastigonemes (flagellar hairs) and confirmed their expression specificity during gametogenesis in S. marinoi25,32 (Fig. 2). Four genes were finally selected as flagella markers: SIG1 (“Sex Induced Gene 1”), DNAH9/11/17 and DNAH5/8 (“Dynein Axonemal Heavy Chain”) and DRC4 (“Dynein Regulatory Complex”) (Table S6).

Collectively, we identified a set of eight sex-specific marker families: one conserved across all major diatom clades, three specific to raphid pennates, and four flagella markers for centric diatoms (Fig. 3d). These markers are well-suited for phylogenetic analysis due to several key attributes. First, each marker family forms a discrete clade of mostly single-copy orthologs (Fig. 3d); secondly, the sequence variability within individual genes is sufficient to distinguish them at the species level. Next, we used profile hidden Markov models (HMMs), which encapsulate the sequence variation in each marker family, to assess their (co-)expression patterns in environmental metatranscriptomic datasets (Fig. 2).

Broad sexual potential during an estuary diatom bloom
To assess the applicability and taxonomic resolution of the sex marker genes, we experimentally induced sexual reproduction during an intense seasonal diatom bloom in the freshwater tidal part of the Scheldt estuary in Belgium, one of the last North-West European estuaries with an intact salinity gradient33. At the time of sampling, the bloom had almost reached its peak density (chlorophyll a ~50 µg/L, Fig. S7) and was dominated by thalassiosiroid centric diatoms34. Inspired by previous laboratory studies showing that salinity can trigger sexual reproduction in certain diatom species35,36, whole-plankton samples (salinity ~0.4 ppt) were manipulated in a controlled microcosm set-up (Fig. 4a). An abrupt shift to a salinity of 10 ppt was introduced to simulate the effect of salt water incursions typical of periods of low river discharge in the estuary. Microscopic observations confirmed that this salinity shock triggered sexual reproduction in the dominant genus Cyclotella (Fig. 4b). In contrast to control conditions, where no sexual stages were detected during the course of the experiment, we observed spermatogenesis and fertilization of egg cells after 24 h of salinity treatment, followed by the development of auxospores after 48 h (Fig. S8–13). At both time points, 18S rDNA metabarcoding confirmed the dominance of Cyclotella, accounting for about 50% of the sequencing library (Fig. 4c). Phylogenetic classification of amplicon sequence variants (ASVs) revealed a single pre-dominant ASV closely related to Cyclotella scaldensis, which constituted 46% of all reads (Fig. 4d, e, Fig. S14). More than 20 other C. scaldensis-like ASVs were identified, but at much lower abundances (<2%). In addition, other Thalassiosirales species from the genera Conticribra, Skeletonema, Cyclotella (Stephanocyclus) and Thalassiosira were present at low abundances (<1%) (Fig. 4d, e, Fig. S14).

We next performed deep metatranscriptomic sequencing (>100 million reads per replicate) to assess the co-expression of sex markers. A de novo transcriptome assembly yielded 868,920 individual transcripts (Figs. S15, S16, S17, Note S2). Among 102,859 expression-filtered transcripts belonging to four main Thalassiosirales clades (Cyclotella, Thalassiosira, Skeletonema, Conticribra), 21,308 were differentially expressed between the salinity treatment and the freshwater control conditions. For each of the sex markers, we used HMM searches to identify homologs in the de novo transcriptome. Phylogenetic analysis was then performed to retain only those transcripts that clustered within the clade of referenced diatom markers (Fig. 2). In general, marker genes were significantly upregulated in response to the sex-inducing salinity treatment, demonstrating their efficacy to detect sexual reproduction (Fig. 4f). All markers showed the highest expression 24 h after salinity treatment, coinciding with flagella formation and meiosis. Phylogenetic analysis, including reference proteomes of 99 different diatom accessions, identified four distinct species that co-expressed at least three sex markers, each belonging to a different Thalassiosirales clade (Fig. 4f, Fig. S18–S22). Cross-referencing the phylogenetic position of the sex markers with the 18S barcoding data confirmed that each sexually active species corresponded to one or multiple ASVs. Interestingly, sexual reproduction was not only detected in the very abundant Cyclotella scaldensis, but also in rare species, such as Thalassiosira allenii/angulata and a Conticribra sp., which represented only about 0.1% of the metabarcoding library.

Widespread diatom sexual reproduction in the global ocean
To detect sexual reproduction in the global ocean, we examined the expression of marker genes across the 52 diatom MAGs obtained from 143 locations sampled during the Tara Oceans expedition (Fig. 2, Fig. S23). Despite being selected for specificity, sex markers tend to show a low level of background expression in non-sexual conditions (Fig. S5, Table S8). Therefore, we only considered marker genes to be sexually expressed when they exceeded the expression level of at least 95% of laboratory RNA-seq samples in vegetative conditions. To identify populations emitting a sexual signal in the global ocean, we focused on diatom MAGs that co-expressed sex markers in addition to the meiotic gene SPO11-2 as a positive control. We tested the effect of requiring different degrees of co-expression (SPO11-2 + 1 marker, SPO11-2 + 2 markers, SPO11-2 + 3 markers, SPO11-2 + 4 markers) (Fig. S24). After matching the predicted sexual Tara stations with imaging data (see below), we decided to consider MAGs co-expressing at least two markers in addition to SPO11-2 as indication for ongoing sexual reproduction within these populations.

Using these criteria, we identified 32 and 181 cases with a co-expression signal for pennate and centric diatoms, respectively (a case referring to a MAG that shows a signal in a single sample). In order to evaluate whether background vegetative expression of sex markers may cause false positive hits, we defined a set of five control marker families. Apart from not being upregulated during sex, control markers have similar characteristics to pennate sex markers M1-M4 and SPO11-2: they belong to mostly single-copy ortholog families encoded in at least three of the reference species and exhibit similar expression levels to sex markers in non-sexual conditions (Figs. S25, S26, Table S9). Repeating the co-expression analysis outlined above, we identified 7 cases where control markers showed a co-expression signal. While this is lower than the sex markers, it shows that false positives may occur under background expression levels. Hence, while co-expression of sex markers data can be an indication that sex is occurring, it should be verified with independent evidence (e.g., imaging) whenever possible.

Surprisingly, co-expression of sex markers suggested ongoing sexual activity in 9 out of 12 known diatom genera across 54 locations in the global ocean (Fig. 5a, Note S4), corresponding to 38% of all investigated stations and 46% of the examined MAGs (Table S10). MAGs exhibited distinct geographical patterns of co-expression, which did not simply reflect the overall abundance and distribution of the respective genera (Figs. S27, S28). Interestingly, for the majority (5 out of 9) of genera, we observed that multiple congeneric MAGs displayed a sexual signal in exactly the same location, suggesting co-regulation of closely related species by similar environmental parameters. Chaetoceros, Fragilariopsis, Leptocylindrus, Odontella, and Skeletonema showed sexual reproduction mostly in cold waters, while signals of sexual events in Cylindrotheca and Pseudo-nitzschia were most often observed in temperate and tropical regions. Finally, Minidiscus and Thalassiosira showed a more widespread distribution of sexual reproduction that ranged from polar to tropical locations (Fig. 5a).

Sexually active Chaetoceros and Skeletonema species occurred across the entire Arctic Ocean. Screening of Imaging FlowCytoBot (IFCB) data acquired during the Tara Oceans Polar Circle expedition confirmed the presence of vegetative cells and 4-to 8-celled spermatocytes of Chaetoceros at Arctic stations 163 and 188 (Fig. 5b). Here, Chaetoceros MAGs co-expressed sexual markers in both the pico- and nanoplankton size fractions, consistent with the diameter of the observed sexual cells. At Arctic station 188, where, in addition to Chaetoceros, Skeletonema MAGs also showed co-expression of the marker genes, we identified a single auxospore of Fragilariopsis in IFCB data, which was supported by a single sex marker co-expressed with SPO11-2. (Fig. 5a, b).

In contrast, other genera exhibited more globally widespread patterns of sexual reproduction. The cosmopolitan Cylindrotheca, the only genus never occurring in cold waters, mainly showed a sexual signal in temperate and subtropical areas, including the Indian Ocean, the Pacific Ocean and the Mediterranean Sea. Pseudo-nitzschia MAGs displayed extensive sexual activity across 10 stations in the tropical Pacific Ocean (Fig. 5a).

The observed sexual signals were then compared with diatom abundance data obtained through metabarcoding of the V4 and V9 18S rDNA regions (Fig. 6a). Chaetoceros exhibited a significantly higher abundance for both 18S regions in stations where it displayed sexual activity (two-sided Mann Whitney U-test; 18S-V4: U(76, 661) = 28,772, p < 0.001, effect size = 0.262, 95% Confidence Interval = [0.116, 0.397]. 18S-V9: U(87, 779) = 40,029, p < 0.001, effect size = 0.294, 95% Confidence Interval = [0.164, 0.415].). Furthermore, chlorophyll a concentration, a proxy for bulk phytoplankton biomass, was high in the Arctic and Southern Ocean stations where Chaetoceros was undergoing sexual reproduction, both at the surface and Deep Chlorophyll Maximum (DCM) (Fig. 6b). While we cannot fully exclude that metatranscriptomic detection of sexual reproduction might be influenced by species abundance, we observed several stations with high Chaetoceros abundance that did not show any sex-specific gene expression and also observed IFCB images of vegetative Chaetoceros cells in samples without a sex signal, suggesting that the presence of sexual reproduction markers is not solely explained by coverage or abundance thresholds. Similarly, Pseudo-nitzschia, Fragilariopsis and Thalassiosira were typically abundant in locations where sexual reproduction was observed, although not significantly more than the sites without sexual activity. In contrast, the co-expression of sex marker genes by locally rare species seems to suggest that other genera were reproducing sexually at very low abundances: the median relative abundance of Odontella, Minidiscus, Skeletonema and Leptocylindrus was consistently below 1% in stations where these genera underwent sexual reproduction (Fig. 6a). Finally, sexual reproduction events were detected at various times of the year (Fig. 6b), although no clear pattern emerged when considering the prevailing environmental conditions at these stations (Fig. S29).

Methods
Retrieval of diatom protein sequences, homology prediction and functional annotation
Proteomes, which represent the complete set of proteins for a species, were retrieved from the PLAZA Diatoms v1.0 platform for the 10 diatom species listed there30. Cylindrotheca closterium protein sequences were sourced from Audoor et al.26 and S. marinoi nucleotide sequences of primary alleles were obtained from Audoor et al.26 and translated to protein by TRAPID 2.061. De novo transcriptome sequences for the basal centric diatom Leptocylindrus danicus were retrieved from the decontaminated Marine Microbial Eukaryote Transcriptome Sequencing Project (MMETSP)62–64. Subsequently, open reading frame prediction was carried out by a Transdecoder v5.0.2 run (github.com/TransDecoder), and homology inference was performed through a Diamond v2.0.14 search against the merged SwissProt v2022-02-09 and PLAZA Diatoms 1.0 proteome in –ultra-sensitive mode65.

We used the homologous gene families defined in PLAZA Diatoms v1.0 to identify homologs across species30. As C. closterium and S. marinoi are missing from the PLAZA Diatoms platform, their genes were assigned to PLAZA gene families through TRAPID 2.0 with PLAZA Diatoms as a reference database and default settings61. The resulting gene family prediction for 12 genomic diatom species and 16 other eukaryotic species was used to determine the taxonomic distribution of gene families. Finally, InterPro domain information was retrieved from PLAZA Diatoms for S. robusta and P. multistriata30, and from TRAPID 2.0 for C. closterium and S. marinoi61.

Differential expression analysis of sexual reproduction in four diatom species
We considered transcriptomic libraries profiling sexual reproduction in four different species: the centric diatom S. marinoi and the pennate diatoms P. multistriata, S. robusta and C. closterium (Table S1)25,26,28,30. To allow inter-species comparisons between time series with varying duration and sampling frequency, we assigned each sample to a mating stage based on the defining sexual cell stage present in each culture: sex pheromone signaling (S), gametangia (P, pre-gametic), gametes/zygotes (GZ) and auxospores (A) (Fig. 1a). The quality of all RNA-seq libraries was verified using FastQC v0.11.2 and summarized with MultiQC v1.766,67. Illumina adapters were detected in P. multistriata and mating type minus S. robusta libraries, which were removed with Trimmomatic v0.36 (ILLUMINACLIP:2:30:10:2)68. FastQC was repeated to confirm successful adapter removal before proceeding to read mapping. Transcript sequences for each of the four species were derived from their respective reference genome assemblies (Table S2) and indexed with Salmon v1.3.0 with a k-mer size of 31, while retaining identical transcripts (–keepDuplicates)69. To improve mapping accuracy, the full set of genomic contigs of each respective species was added to the indexed transcriptome as decoys. Subsequently, expression quantification was performed for all samples with Salmon v1.3.0 in selective alignment mode and correcting for biases using the flags –gcBias and –seqBias69. Salmon read counts were imported into R using the tximport package70. For P. multistriata, S. robusta and C. closterium, tximport was used to summarize transcript counts on the gene level, while S. marinoi reads were already mapped to gene-level transcripts. After filtering to keep only genes showing an expression level of more than one count per million (CPM) in at least three samples and trimmed mean of M value (TMM) normalization, negative binomial generalized linear models were estimated for each gene using EdgeR 3.30.371. Differential expression (DE) tests were performed to contrast the response between sexual and control samples for each species in each cell stage. For P. multistriata, gene expression of crossed cultures at each stage was compared to the average response of control cultures over both mating types, while for S. robusta and C. closterium, expression in sexual cultures was compared with their time-matched vegetative control. Meanwhile, for each S. marinoi cell stage, salt-treated sexual samples with a cell size below the sexual size threshold were compared with corresponding salt-treated cultures of large cell size (Table S1). For each contrast, a likelihood ratio test (LRT) was performed relative to an absolute log2 fold change threshold of one with the glmTreat function in EdgeR. We adopted a stage-wise testing procedure as introduced in the StageR package for R to determine an “overall” set of genes at the species level, which are DE in at least one sexual cell stage, while controlling the overall false discovery rate (FDR) at 5%72,73. To this end, p values for individual contrasts were aggregated on the gene level using the Sidak method from the Aggregation package74. FDR-adjusted p values of genes that passed the screening stage were further adjusted using the stepwise Holm procedure.

Data-driven discovery of conserved markers for sexual reproduction in diatoms
Homologous gene families allowed the integration of differential expression results between species to select sex-specific diatom genes families that can serve as markers. We first inspected the general association of gene expression across stages and species by calculating the Pearson correlation between the average log2 fold changes of genes differentially expressed in each gene family. Next, a five-step procedure was used to identify conserved marker genes for sexual reproduction. First, we restricted our search to diatom-specific gene families that were expressed in at least three out of four sexual diatom species. The resulting families were further filtered by keeping the ones for which all expressed genes were significantly upregulated during sex in all species where they occurred. The remaining families were ranked by the average log2 fold changes observed in the cell stages where they were significantly differentially expressed. In a final step, the sex-specific expression of the top 10 ranked gene families was verified in the 119 vegetative RNA-seq samples from the S. robusta expression atlas, spanning 42 different experimental treatment groups29,30,75,76 (Table S6). Reads from the expression atlas were mapped, summarized at the gene level and normalized using TMM normalization, following the previously described pipeline. After inspecting gene expression levels in the expression atlas, the four top-ranked families (M1-4) with the highest mean fold change ranks were selected as specific markers of sexual reproduction in diatoms (Table S7).

The protein sequences of the M1-4 families were used as an input for profile HMM construction. Multiple sequence alignments for each marker gene were carried out using the MUSCLE algorithm77 implemented in the software MEGAX, version 11.0.878, and then transformed into profile HMMs using the hmmbuild command of the HMMER package, version 3.1b279.

Identification of meiotic genes and flagella markers
Based on literature, five diatom meiosis-specific genes with no known function outside meiosis were selected80: SPO11-2, MER3, MND1, MSH4 and MSH5. The DNA repair gene RAD51-A was added as an additional meiotic gene due to its strong upregulation in P. multistriata28 (Table S4). For each meiotic gene, we identified the corresponding homologous gene family on PLAZA Diatoms v1.0 and pinpointed the relevant clade containing diatom sequences and an Ectocarpus siliculosus outgroup (Table S4)30. Potential homologs in L. danicus, S. marinoi and C. closterium were identified using a relaxed BLASTp search (E < 1e-5) followed by construction of phylogenetic trees with IQ-tree (explained below) to verify correct clustering within the diatom clade.

Amino acid sequences of diatom flagellar proteins were extracted from two studies25,32. Flagella genes previously identified in L. danicus64 were not considered because of their lack of expression conservation in S. marinoi25. After identification of the homologous gene family of the best BLAST hits in PLAZA Diatoms v1.0, we inspected their phylogenetic trees to select the diatom orthologs that belonged to the same subclade as the original flagella proteins, as well as outgroup sequences for phylogenetic analysis in related species E. siliculosus or Aureococcus anophagefferens (Table S5). Flagella proteins were named after the Homo sapiens homologs. For each flagella protein, homologs in S. marinoi and L. danicus were identified using BLASTp (E < 1e-5); the clustering of the top hit(s) in the correct diatom clade was verified by constructing phylogenetic protein trees with the PLAZA diatom homologs and outgroups. As flagella markers, we selected only those genes that were expressed during sex and significantly upregulated in the differential expression analysis of S. marinoi sex transcriptomes.

For each of the flagella and meiotic genes, profile HMMs were constructed using all identified homologs from diatoms. After multiple sequence alignment with MAFFT v7.453 in –auto mode81, HMMs were constructed using the hmmbuild command from HMMER v3.1b2.

Identification of control marker gene families
To assess the possibility that sex marker co-expression in the Tara dataset may be caused by background expression in vegetative conditions, we defined a set of five control markers: C1-C4 (similar to pennate sex markers M1-M4) and PC (“positive control”, similar to SPO11-2) (Fig. S25). Control marker families were selected to have similar characteristics to the sex markers, including an expression level similar to that of the sex markers during vegetative, non-sexual conditions. Specifically, control gene families (1) must be encoded in the genome of at least three of the four reference species: S. robusta, C. closterium, P. multistriata and S. marinoi; (2) must consist of single-copy orthologs, allowing just a single duplication in one reference species; (3) must not show upregulation during sex in any of the four species. The average cross-species expression (in TPM) and its coefficient of variation (CV) were then compared to the average expression level and CV of each pennate sex marker (M1-M4, SPO11-2) in non-sexual conditions (Table S8). After removing families with poor phylogenetic resolution (e.g., HOM02SEM051941, the lowest expressed family in the dataset that encodes a suspected retroviral transposon), five control marker families (C1–C4 and PC) were selected (Table S9). Since only a few among the 2734 possible control families showed an expression level as low as the sex markers (Fig. S26), the final set of control markers had a higher expression level than the corresponding sex markers (Tables S8, S9).

Microcosm experimental design, RNA extraction and sequencing
Phytoplankton abundance and bloom dynamics in the Scheldt estuary were monitored throughout the year 2022 by taking regular Chlorophyll a measurements. Samples were filtered over 25-mm diameter Whatman GF/F glass fiber filters, followed by High Performance Liquid Chromatography (HPLC) using the method developed by Van Heukelem and Thomas82. On 20/05/2022, during the early stages of a bloom of thalassiosiroid diatoms, a natural phytoplankton sample was collected from the freshwater zone of the tidal Scheldt estuary near the city of Ghent (51°00′15.6″N 3°48′19.3″E, map created with ggMap83). The sample was filtered through a plankton net with a mesh size of 10 µm to concentrate the plankton which was then distributed over large 150-mL culture flasks. Culture flasks were supplemented with nutrients making up WC medium with half-strength N, P and Si: 42.5 mg/l of NaNO3, 4.4 mg/l of K2HPO4 and 14.2 mg/l of Na2SiO3.9H2O and kept at 21°C in a 12/12 day/night rhythm at 5 µmol photons m−2 s−1. Half of the cultures were brought to a salinity of 10 ppt. To monitor the progression of sexual reproduction and identify species present in the sample, microscopic pictures were taken by a ToupTek TP105100A digital camera attached to an Axiovert 40C inverted microscope fitted with a LD ACHROPLAN 20x objective and a 0.5x C-Mount (Carl Zeiss GmbH, Oberkochen, Germany).

After 24 h and 48 h, two salt-treated and one control samples were harvested by filtration on a Versapor filter with a pore size of 3 µm (Pall Corporation, NY, USA) and flash frozen in liquid nitrogen. Total RNA extraction was carried out using a RNeasy Plant Mini Kit (Qiagen) as detailed by Bilcke et al.29. RNA purity and quantity was assessed by NanoDrop and a bioAnalyzer Eukaryote Total RNA pico chip (Agilent) was run to verify RNA integrity. Unstranded library preparation was performed using the NEBNext Ultra RNA Library Prep Kit for Illumina including a polyA selection step to enrich for messenger RNA. Finally, 150 bp paired-end sequencing of reads was performed on the NovaSeq 6000 S4 platform at Novogene Ltd. (UK).

Microcosm 18S metabarcoding
To assess the community composition throughout the microcosm experiment, we collected a natural sample at the sampling site and two replicates for each treatment for metabarcoding analysis. DNA samples were collected by vacuum filtration water until saturation, using a 0.22 µm MF-MilliporeTM filter. Immediately afterwards, the filter was flash frozen and stored at −80°C, until further analysis. To assess the natural community of the sampling site, one sample was collected in the field at the moment of sampling.

DNA extraction was carried out with the DNeasy Powerlyzer Microbial Kit from Qiagen (Hilden, Germany). The polymerase chain reaction (PCR) amplification was targeting the ribosomal RNA gene. The 18S V4 region was amplified using the TAReuk454FWD1 (5′-CCAGCASCYGCGGTAATTCC-3′) and the TAReukREV3 (5′-ACTTTCGTTCTTGATYRA-3′) primers84. PCR and library preparation were performed following D’Hondt et al.85. Paired-end (2×300 bp) sequencing was performed with the Illumina MiSeq technology (Illumina, San Diego, US) by Genewiz (Leipzig, Germany).

We used the DADA2 pipeline to assign sequencing data to amplicon sequence variants (ASVs) for each sample86. Quality control, trimming and filtering was done using the ‘filterAndTrim’ function. TrimLeft removed the primers and trimRight removed 10 bp at the end of each read. Afterwards, forward and reverse reads were truncated after 250 and 220 nucleotides (truncLen). Sequences with EE higher than 2 and ambiguities (maxN) were removed. Pair ends were merged with a minimum overlap length of 12 bp, allowing for one mismatch using the function ‘mergePairs’. The 18S sequences were assigned to the Protist Ribosomal Reference database (PR2) version 4.14.0 (https://github.com/pr2database/pr2database/releases/tag/v4.14.0)87. The ASV table was further analyzed using the phyloseq package88. Contaminant sequences were removed following the method of Davis et al.89. Reads belonging to non-identified taxa at kingdom and infra-kingdom, as well as to the phylum Metazoa and the classes Embryophyceae and Streptophyta were deleted. ASVs with a relative abundance lower than 1e-5% were omitted.

A deeper classification of the 32 Thalassiosirales ASVs that made up more than 0.1% of the library in at least one sample was performed by constructing a phylogenetic tree that combines ASV sequences with reference Thalassiosirales 18S nucleotide sequences from NCBI GenBank. Multiple sequence alignment was carried out with MAFFT v7.453 and a phylogenetic tree was constructed with IQ-tree v2.1.390 with 1000 ultrafast bootstrap repeats. Barplots showing ASV abundance were added to this tree using the ggtree and ggtreeExtra packages for R91.

De novo assembly, taxonomic classification and differential expression analysis of Scheldt metatranscriptome
The quality of raw Illumina reads was checked with FastQC v0.11.2 and summarized with MultiQC v1.766,67. Raw reads were normalized with the bbnorm function from BBMap v38.98, with settings target=50 min=20 minkmers=1592. Reads of the six metatranscriptome samples were pooled to produce a single de novo co-assembly using rnaSPAdes 3.14.1 with a k-mer size of 4993, resulting in a set of 868,920 transcripts with an average length of 1113 nucleotides and an N50 value of 1545 bp. Detection and translation into protein of open reading frames (ORFs) was performed with Transdecoder v5.0.2 (available at github.com/TransDecoder). Homology-based ORF detection was enabled through Diamond v2.0.14 searches in –ultra-sensitive mode65 against the SwissProt v2022-02-09 and PLAZA Diatoms 1.0 proteomes.

For taxonomic classification of reads and transcripts, a custom diatom-enriched Kraken2 v2.1.2 database94 in protein mode was created combining Kraken2 RefSeq protein sequences (Archaea, Bacteria, human, fungi, plant and protozoa), all PLAZA Diatoms v1.0 proteomes30, the decontaminated MMETSP sequences of 53 diatom species, 87 Thalassiosirales-related proteomes95, 7 diatom proteomes from Phycocosm (phycocosm.jgi.doe.gov: Maypse1, Fracro1, Conwei1, Chaten1, Nitput1, Nithil2, Cylclo1) and 3 Pseudo-nitzschia proteomes96. The MMETSP transcriptomes were first translated by transDecoder as explained above for L. danicus62,63. We then used Kraken2 v2.1.294 to taxonomically classify the 868,920 rnaSPAdes metatranscripts, and in parallel to classify the raw reads for control (2 concatenated samples) and sexual (4 concatenated samples) conditions in paired-end mode.

Expression quantification was performed by mapping back all paired-end reads to the indexed rnaSPAdes de novo transcriptome (k-mer size 31) with Salmon v1.8.069 (settings: --gcBias, --seqBias, -l A). After importing read count tables with tximport, genes with a CPM > 0.1 in at least two samples were retained, followed by TMM normalization. Next, negative binomial generalized linear models were fitted and a LRT was performed in EdgeR v3.30.3 using the glmLRT function, comparing gene expression between salt-treated and control samples71. To correct for changes in the relative abundance between samples, we used the Kraken2 taxonomy to perform the differential expression pipelines separately for each of the four major Thalassiosirales clades: Cyclotella (29,910 filtered transcripts), Thalassiosira (31,895 filtered transcripts), Conticribra (17,418 filtered transcripts) and Skeletonema (23,636 filtered transcripts). While assembly-wide differential expression data was initially used to determine which taxa co-expressed sex markers (Figs. S18–22), the abundance-corrected data were used to further investigate the co-expression of sex markers within each taxon (Fig. 4). Raw unscaled abundance-corrected CPM expression values of sex marker hits are displayed in Fig. S30.

Detection and phylogenetic analysis of marker genes in Scheldt metatranscriptome
To detect diatom homologs of sex marker genes in the metatranscriptome assembly, HMM searches were performed. Profile HMMs of the data-driven sex markers (M1-4) and flagella markers (Fig. 2) were queried against the translated rnaSPAdes metatranscriptome using HMMER v3.1b2 (bitscore > 50 and E < 1e−10). Protein sequences of hits were combined with the original diatom marker proteins and selected outgroup sequences to delineate the specific clade of interest. Multiple sequence alignments were produced with MAFFT v7.45381 and gaps were trimmed with TrimAl 1.4.197 with a gap threshold of 0.5 except for DNAH flagella proteins, where a gap threshold of 0.2 was used due to the large number and diversity of HMM hits. Maximum-likelihood phylogenetic trees were constructed based on trimmed alignments with IQ-tree v2.1.390 with 1000 ultrafast bootstrap repeats and automatic ModelFinder substitution model selection. Scheldt metatranscriptome proteins clustering within the clade of the original marker genes were selected (“phylogenetic selection”, Fig. 2) for downstream assessment of their differential expression, taxonomic classification and gene expression levels.

Scheldt marker transcripts were phylogenetically situated in the Thalassiosirales using 87 proteomes of Thalassiosirales and related centric diatoms95. Sex marker homologs were selected from these 87 proteomes using HMM searches and phylogenetic selection as described above. These reference markers and the original marker protein sequences (homologs in 12 reference diatom genomes and potential non-diatom outgroups) were merged with the selected Scheldt hits, followed by multiple sequence alignment, trimming and phylogenetic tree construction as explained in the previous paragraph. Phylogenetic trees were midpoint rooted with the Phangorn package for R. Trees were visualized with annotation and metadata using the ggtree and ggtreeExtra packages for R91.

Detection of marker genes in Tara Oceans MAGs
We next identified homologs for each of the eight marker genes and SPO11-2 in the MAGs reconstructed from Tara Oceans metagenomic reads98 (Fig. 2, Table S11). In brief, protein sequences encoded on a total of 683 MAGs were retrieved from https://www.genoscope.cns.fr/tara, out of which 52 diatom MAGs were searched using the hmmsearch command from HMMER v3.1b2. In the case of the gene M1, we masked the carboxy-terminal Tubby-like domain before querying the Tara Oceans database to avoid that the widespread and phylogenetically conserved domain resulted in false positives. The output hits of sexual marker genes in MAGs were filtered by selecting only sequences belonging to diatom MAGs, showing a HMMsearch score ≥ 50 and an Evalue lower than 1e−10.

Next, a phylogenetic selection was performed to select robust diatom homologs of each marker (Fig. 2). First, the protein sequences of the selected Tara hits were merged with the query diatom markers (Table S7) and, when available, homologs from L. danicus and Chaetoceros tenuissimus99. To improve the separation of the correct marker clade, we added diatom sequences from the mitosis-specific SPO11-3 paralog to the alignment of the meiotic gene SPO11-2. Similarly, in the case of the two closely related markers that code for the dynein axonemal heavy chain (i.e., DNAH5/8 and DNAH9/11/17), we always included the reference sequences of the sister markers as an outgroup.

Multiple sequence alignments were then performed using MAFFT v.7.49081 followed by automatic removal of poorly aligned positions with trimAl v.1.4.197 and maximum-likelihood phylogenetic tree construction with IQ-TREE v2.1.390, using the same settings used for the microcosm experiment. Trees were midpoint rooted and visualized and edited in R (packages: ggtree, treeio, ape, phylotools) and FigTree version 1.4.3100. For each phylogenetic tree, all homologous Tara genes clustering within the original marker clade were selected, taking into account tree topology and bootstrap support values (Fig. 2).

Determining co-expression of sex marker genes in Tara Oceans
After identifying sex marker genes on diatom MAGs, metatranscriptomic data associated with each MAG were analyzed to identify diatom sexual reproduction events at a global scale. Tara Oceans metatranscriptomic reads corresponding to various size fractions were mapped against the MAGs using BWA-MEM2 for high-efficiency mapping, followed by processing with Samtools for sorting and duplicate removal. MAGs expressing less than 10,000 reads or 1000 genes in a sample were filtered out. MAG-level transcripts per million (TPMs) were calculated for each gene relative to the total number of reads mapped to each MAG in each sample. By calculating MAG-level TPM values, we corrected for the relative abundance of each MAG across samples, ensuring that the expression of sex marker genes was normalized relative to the overall transcriptional activity of the same MAG. This enabled consistent comparisons across stations and provided a robust basis for applying gene-specific expression thresholds determined for single species in laboratory setting, as explained below.

In parallel, expression thresholds were calculated based on the four-species bulk RNA-seq dataset (Table S8). Salmon TPM values were imported into R using tximport70. For each sex marker family, the 95% TPM percentile was calculated across all genes in that family during vegetative conditions. This amounted to 6 vegetative samples per gene for P. multistriata, 119 vegetative samples per gene for S. robusta (expression atlas), 9 vegetative samples per gene for C. closterium and 8 vegetative samples per gene for S. marinoi (4 salt-treated and 4 untreated controls).

Next, the Tara expression dataset was filtered to retain only samples in which diatom MAGs expressed the positive control gene SPO11-2. For those MAGs, the co-expression of sexual markers was defined as the case where multiple sex marker genes encoded on the same MAG had MAG-level TPM values exceeding the 95th percentile expression threshold (Table S8) and were supported by at least two reads within the same Tara sample (Table S12). Co-expression of control markers C1-C4 with their positive control PC was determined in the same manner (Table S9). The geographical distribution of sex signals was visualized using R with ggplot2101 and Scatterpie102 packages, by grouping the information at the genus level. Assembly metrics for Tara Oceans diatom MAGs were compared between sex and non-sex samples using the Mann-Whitney U test (Fig. S31, Table S13). The following metrics were analyzed: ANVIO completion, that represents the proportion of the genome that is assembled, reflecting the completeness of the MAG assembly; ANVIO redundancy, referring to the redundancy of the MAG assembly and thus to the level of duplicate or repetitive sequences within the assembled genomes; N50, i.e., a measure of assembly contiguity, defined as the length of the shortest contig at which half of the total assembly length is contained in contigs of that length or longer; num_contigs, i.e., the total number of contigs in the MAG assembly; total length, i.e., the total base pair length of the assembled genome. To account for multiple comparisons, p values were adjusted using Bonferroni correction.

Integration of metadata associated with Tara Oceans
The distribution of detected sexual reproduction events across sampling stations was analyzed for different size groups and depths. We considered the following size classes: 0.8–5 µm (picoplankton), 3–20 µm and 5–20 µm (nanoplankton), 20–180 µm (microplankton), 180–2000 µm (mesoplankton), and the following depths: surface (SRF) and DCM (Deep Chlorophyll Maximum). Furthermore, we integrated information derived from metabarcoding data (V4 and V9 regions of the 18S rDNA) to further explore the putative relation between sexual reproduction events and species abundance (Note S3). OTUs tables for V4 and V9 regions were downloaded from https://zenodo.org/records/4043757 and https://zenodo.org/records/3768510, respectively. For each genus, a Mann-Whitney U test was applied to assess differences in distribution of abundances, with the Bonferroni correction to account for multiple tests. An adjusted p < 0.05 was considered statistically significant. Results are reported as U(m, n), p, effect size, and 95% confidence intervals. Full statistical values are reported in the Results for Chaetoceros (V4 and V9). Finally, in order to identify common conditions in locations with diatom sexual reproduction events, we analyzed environmental parameters, including seasonality. Chlorophyll a concentrations were downloaded from https://doi.org/10.1594/PANGAEA.875579. To compensate for missing physicochemical values in the Tara Oceans in situ dataset, climatological data of the World Ocean Atlas 2018103 were preferred.

IFCB data104 from the Tara Oceans Polar Circle expedition were explored through the EcoTaxa105 platform. All publicly available images as of 19 December 2023 were screened for the presence of vegetative cells and sexual stages (spermatogonia, oogonia and auxospores). Afterwards, observations of vegetative and sexual stages were cross-referenced with the predicted sexuality at each station of the nine diatom genera that co-expressed sexual markers (Chaetoceros, Cylindrotheca, Fragilariopsis, Leptocylindrus, Minidiscus, Odontella, Pseudo-nitzschia, Skeletonema, Thalassiosira). Results of this microscopic screening are available on Zenodo (see “Data Availability”).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Reviewer's conclusion\nThe submitted work gives the impression of a very profound study, touching upon one of the key stages of the life cycle of\ndiatoms – the stage of sexual reproduction.\nAs for the essence of the study, it should be noted here its versatility and complexity, which was achieved by involving\nspecialists of different profiles. In this regard, it should be emphasized that for me, as a narrowly focused specialist, it is\nextremely difficult to give an exhaustive assessment of all the issues discussed, from problems of reproductive biology to\nitems of molecular genetics and bioinformatics.\nSince I am sufficiently familiar with the field of reproductive biology of diatoms, I can note that the presented work contains a\nnumber of new findings that significantly expand our understanding of the mechanisms of regulation of the sexual process\nand sexuality of diatoms.\nA set of genes were identified that are differentially expressed during sexual reproduction compared to vegetative\nconditions. What is very important, a set of expressed genes related to sexual rerproduction was established that is universal\nfor centric and pennate species of diatoms. At the same time, marker genes characteristic only of centric and only of pennate\nwere identified. These genes allow the sexual process to be registered directly in a natural population.\nUsing this approach, the authors were able to analyze the results of metagenomic analysis performed on large areas of the\nWorld Ocean. This allowed the authors to obtain evidence of the commonness, ubiquity, and high frequency of sexual\nreproduction of particular diatoms in the World Ocean. Previously, cases of sexual reproduction and auxospore formation\nwere considered rare events.\nAnother important result of the work can be considered the confirmation that the identified genetic markers are well-suited for\nphylogenetic analysis due to several key attributes. First, each marker family forms a discrete clade of mostly single-copy\northologs; secondly, the sequence variability within individual genes is sufficient to distinguish them at the species level.\nFrom the standpoint of understanding the species and the species concept, the identified marker genes involved in sexual\nreproduction are much more acceptable for establishing species boundaries than the commonly used chloroplast or\nribosomal genes. The accumulation of nucleotide substitutions in the latter is not sufficient evidence of completed\ndivergence in the genetic history of species; they merely correlate with the emerging, ongoing, or completed divergence of\nspecies.\nThe manuscript does not contain any formal errors; the text is grammatically and stylistically verified and does not require\nadditional revision. The work meets the expected standards, at least in field of reproductive biology of diatoms. The authors\nprovided enough detail in the methods for the work to be reproduced\nFrom my point of view, the paper can be published in the current form.\nNickolai Davidovich", "reviewer_2": "• What are the noteworthy results?\nThe life cycles of phytoplankton are probably their least understood characteristics. Diatoms are among the most\necologically and biogeochemically important phytoplankton, one of the most striking groups, and yet their life cycle has\nremained strikingly mysterious. The structure and process of formation of the silica frustule in most diatoms is such that, with\nevery cell division, one of the daughter cells is slightly smaller than the parent, diminishing the average size of a clone\nthrough time. Cell size in most species studied is only restored through a sexual phase. This also means the sexual phase\nimplies a loss of carbon and silicon from the phytoplankton to sedimenting silica valves as new ones are synthesized. It has\nlong been pointed out that there are diatoms that manage to avoid the size reduction, and also making sexuality no longer\nobligate, and this means that resolving the diatom life cycle also might offer insight into one of the biggest questions of\neukaryotic biology, which is why sex exists and persists in spite of obvious costs. The problem is that the sexual stages in\ndiatom life cycles are difficult to observe, so there is very little information from natural environments.\nThe present study offers a major breakthrough, that has been years in the waiting. Early differential gene expression studies\nstaring over 25 years ago, then by population genetics studies and, in the past decade, with full transcriptomic studies, have\nlooked for molecular markers of the diatom life cycle that might allow detection of diatom sex in the field. So this study really\nbuilds on 25 years of work to propose robust sets of markers based on a synthesis of careful lab studies with new models,\nand, for the first time, uses mesocosm and direct sampling of the environment to look for diatom sex across essentially the\nentire ocean.\nUsing a combination of the new molecular markers in the Tara dataset, combined with imaging cytometry, they clearly detect\nevents of mass diatom sex in two arctic ocean stations. This confirms the sparse and chance observations from traditional\nmicroscopy based methods that mass sex events do occur and that they can involve many very similar species\nsimultaneously (itself surprising). Demonstrating a possible toolset for sex detection is by itself an impressive achievement\nthat the diatom community has been waiting for a long time for, and this alone is worthy of publication at a high level\nbecause it means we might finally have what we need to effectively study how diatoms really work in the ocean.\nThe final result, which is potentially most surprising, is the fact that they report detection of sexual markes in a large portion\nof stations, suggesting it is not that rare. if in fact diatom sex occurs not just in very brief pulses, but also is occurring\ncontinuously at a low level as the authors propose, that would be a very important finding both for diatom biology and, even\nmore fundamentally, for biological oceanography.\n• Will the work be of significance to the field and related fields? How does it compare to the established literature? If the work\nis not original, please provide relevant references.\nAs discussed above, I consider the work significant. It might be of great significance. The reason it might be signficant\nbeyond phytoplankton biology and oceanography is, for diatom sex to work, especially in raphid pennate diatoms which (at\nleast in the lab) require opposite mating types to pair, the ocean must have an unseen structure. Could this structure perhaps\nbe the mucus parachutes identified in the latest Science paper by the group of Manu Prakash (Chajwa et al., Science 386,\neadl5767 (2024). DOI: 10.1126/science.adl5767)? Or something else?\nIt also could have importance for understanding the evolutionary maintenance of sex in eukaryotes.\n• Does the work support the conclusions and claims, or is additional evidence needed?\nMost of the work is very rigorous and well supported. I do have one major concern, which is about the interpretation of sex\nmarkers detected in metatranscriptomes when there are no images (e.g., by IFCB) to confirm this. The authors interpret this\nresult as supporting possibility that sex is occuring at a low level most of the time. However, their own transcriptomic data in\nthe lab (Fig. 3 and supplementary figures S4 and S5) as well as previous transcriptome studies both in diatoms as well as\nmany other eukaryotes, show that transcripts of even the most sex-specific genes usually can be detected at low levels from\nvegetative cells. That implies that, even if only vegetative cells are present, reads from such genes would eventually be\ndetected in the metatranscriptome is deep enough.\nTo try to rule this out, the authors arbitrarily considered detection only when the meiotic recombination subunit spo11-2 as\nwell as at least two other sexual phase markers were simultaneously detected. This would be expected to be a very robust\nstrategy. However, again, if the metatranscriptomes are deep enough, eventually one might expect at least some cases of\ndetection of reads from spo11-2 and at least two other genes even when there are only vegetative cells present. I could not\ntell from the present manuscript if they were able to assess this possibility.\n• Are there any flaws in the data analysis, interpretation and conclusions? Do these prohibit publication or require revision?\nDefinitely the issue raised above is quite complex when using metatranscriptome and MAG data, with great variability in\nmapping, and the text implies that the authors have put in an important effort to trying to resolve the problem. However, it is\nnot clear the size and estimated completeness of MAGs, nor do they report how much non-sexual-specific transcripts map to\nthe MAGs. Reporting such information more clearly could help.\nI think there might at least two possible strategies that might help to make their case more convincing.\nIn one approach, they might select genes that are expressed at similarly low levels in vegetative cells, but that are not\nupregulated in the sexual phases. How often are three of such genes from the same organism (MAG) detected at the same\ntime? Clearly such an approach in itself could have weaknesses. It may be that the conditions under which the genes\nselected are expressed are found frequently in the environment but not in the lab. But, if enough low-expressed vegetative\ngenes are selected, one could use a type of bootstrapping by testing all possible three-gene combinations of a much larger\nset, assuming that it is unlikely that all similarly low-expressed genes in the lab are stimulated by the same environmental\ntrigger in nature.\nAnother possibility might be to conduct simulations with data using the sexual markers they have selected and the lab\ntranscriptomic data they have. Perhaps it would be possible to estimate (by repeated simulations; again, perhaps a type of\nbootstrapping), how likely it would be to detect coexpression of the sexual phase genes in the metatranscriptomes from\nvegetative cells, assuming the species contributing to the metatranscriptoms have similar enough vegetative cell\ntranscriptomic profiles to the representatives in the lab.\nArguably both approaches suffer that they might not be strictly quantitative because of the need to extrapolate expession\ndata from four species in the lab to the thousands of species in the field. Nevertheless, it still would be very useful. Perhaps\nthe vegetative expression levels in the lab are higher, perhaps lower than other species in nature, but it would at least offer an\napproximate estimation of the likelihood of false positives (are they likely, unlikely, or very unlikely?).\n• Is the methodology sound? Does the work meet the expected standards in your field?\nI find most of the work quite rigourous. Some appears to be based on synthesizing several previously published lab results,\nand is quite rigorous. The mesocosm approach is quite strong also. The results from the arctic stations, complemented by\nIFCB, is quite convincing, although I have some minor comments below that might help.\n• Is there enough detail provided in the methods for the work to be reproduced?\nYes. But I do suggest providing more details on the MAGs and metatranscriptomic data in this paper, as mentioned above.\nIn summary, I think this work should be tracked for eventual publication. I think they can probably do some important revision\nto address my one principal concern. As I point out, they may be able to make the case qualitatively convincing by offering\nsome more details, but I do hope they can do one of the\nMore detailed comments, suggestions, and corrections follow:\nFig. 1. It’s necessary to go to Fig. S1 to understand abbreviations. Could be nice to have in the main Fig. 1 or the text to\nmake it more fluid. S: pheromone signaling, P: gametangia, pre-gametic, GZ: gametes and zygotes, A: auxospores)\nLines 84-86: “S. marinoi exhibited a smaller number of downregulated genes, which is likely the result of the lack of a sex-\ninduced growth arrest in this species.” That could be an explanation, but it could also be just that only a portion of the\npopulation is induced, so upregulation is easy to detect than downregulation. That is, if 50% of the population is induced and\n50% isn’t, than in 50% of the cells there will be genes which are completely unexpressed in the vegetative stage to maximal\nexpression, and in those cells there will be genes that go from maximal expression to essentially 0. But at the population\nlevel that would be an average change of 0 to 50% max expression for the most upregulated genes but only a change from\n100% max expression to 50% max expression for the most downregulated genes. There might also be differences between\ncentrics and pennates in how the regulation works.\nParagraph starting on line 181 refers to Fig. 3, when I think they mean Fig. 4.\nLines 191-193: “High levels of spermatogenesis and fertilization of egg cells occurred… followed by the development of\nnumerous auxospores after 48h” What are “high” levels? Is there any way to quantify this? If not, it would be perhaps better\nto use a phrase that is more clearly qualitative or semi-quantitative, such as that stages were easily detected in the treatment\nbut not the control. Would be good at least in supplementary figures to have qualitative images from treatment vs control.\nFig. 4F and associated text. F) boxplots comparing the expression levels during sexual reproduction of markers belonging to\neach of the four taxa, calculated as the mean CPM after 10 ppt salinity treatment.” It’s not so clear that if the combined reads\nof all markers is the best way. That is, the transcripts per gene per sexual cell are likely not the same for all genes (I would\nexpect a cell would might need to synthesize less meiotic recombination proteins than flagellar components, for example, to\nget the job done). Is there a way to make some normalized index that could take into account info from the lab transcriptomic\ndata? It might be a way to make it possible to evaluate more precisely the contributions of different species to the sexual\nmetatranscriptome. This isn’t strictly necessary for the purposes of this study, but it could be very relevant to see more\nprecisely how sex could be simultaneous or not among closely related co-occuring species.\nLines 281, 320-329: A challenge, as I understand from Notes S3, is that their methodology relies only on detection of marker\ngenes, as they did not find a way to do relative expression. Their strategy to rely on co-expression of spo11-2 with at least\ntwo other markers is very valuable, and it is quite convincing in the stations where IFCB images confirmed sexual stages\ncould be detected from the corresponding genera.\nI offer some technical suggestions above when discussing it. Whether or not these or similar suggestions can really be\nimplemented for the present study, I would suggest more cautious language. For instance, “In contrast, other genera were\nsexually active at very low abundances” (line 281) might be re-worded to acknowledge explicitly that they have only\ndetected an indicator suggesting sexual activity. More importantly “Our findings challenge this notion by identifying a sexual\nsignal across the overwhelming majority of Tara Oceans stations worldwide, thus suggesting that diatom sexual events\noccur frequently in the global ocean” might also be toned down."}, "author_response": "We thank the Reviewer, dr. Davidovich, for his thoughtful review and especially for offering his\nthoughts on the use of sex markers for species delimitation through phylogenetic analysis. In fact,\nin a separate study in preparation, we are trying to identify sex-specific genes that are diverging\n(i.e., evolving different sequence variants) during an ongoing and incomplete speciation event in\nSeminavis robusta. In general, we expect that cell-cell recognition proteins may be the first to\ndiverge during prezygotic isolation, but it should be noted that the diverging genes may differ from\ncase to case, especially between different mating systems (heterothallic vs homothallic diatoms).\n\nWe sincerely thank the Reviewer for their supportive comments. We especially appreciate the\nsuggestions to improve our analyses, which we have implemented as described below.\n\nHow the “immotile” raphid pennate diatoms manage to pair up in the water column is indeed a\nfascinating topic, and the observed co-expression of sex markers in several planktonic pennate\ndiatoms supports the idea that they are able to do so. Previously, a couple of mechanisms have\nbeen proposed: passive pairing due to sinking [1] or localization in thin layers to improve pairing\n[2]. However, mucus parachutes are another very promising possibility as a way to bring cells\ntogether and prevent them from sinking too quickly. We added this hypothesis to the paper!\n[1] Font-Muñoz, J. S. et al. Collective sinking promotes selective cell pairing in planktonic pennate\ndiatoms. Proceedings of the National Academy of Sciences of the United States of America 116,\n15997–16002 (2019).\n[2] Scalco, E., Stec, K., Iudicone, D., Ferrante, M. I. & Montresor, M. The dynamics of sexual\nphase in the marine diatom Pseudo-nitzschia multistriata (Bacillariophyceae). Journal of\nPhycology 50, 817–828 (2014).\n\nThis is a pertinent remark, and something that we had been considering as well. We would like to\nnote that the Supplementary Figure S4, which the Reviewer refers to, shows the expression of\n“meiotic-specific” genes (as defined in [1]), not sex markers. For these meiotic genes, there is\nindeed a relatively high background expression in vegetative cells, possibly due to high levels of\nmitotic recombination in diatoms [2]. In fact, this result was the initial reason why we did not use\nknown meiotic genes as markers for sex, but rather chose a data-driven approach to identify\nmarkers that are more specific for sex.\nIn response to this question, we gathered all available data from the lab RNA-seq experiments\n(including the S. robusta expression atlas) and plotted all expression for vegetative (control) and\nupregulated sexual (mating) conditions to get a better sense of the background expression of the\nproposed sex markers. These plots have been added as a new Supplementary Figure.\n\nWhile the data-driven approach identified more specific sex-related genes than the meiotic genes\nshown in Supplementary Figure S4, the Reviewer raises a valid point that even for these markers,\nthere is some expression in vegetative cells (for instance, the considerable expression of marker\nM3 in Cylindrotheca closterium). To avoid false positives, we relied on the co-expression of\nmultiple independent markers, but we agree that this may not be sufficient to fully control for false\npositives. It is not an easy problem to solve: due to the complexity of Tara metatranscriptome data\nand the associated MAGs, there is no obvious background model of false positives that we could\ncalculate. However, based on the Reviewer’s suggestions, we improved our bioinformatics\npipeline to assess and reduce the level of false positives (see below).\n\nFinally, in response to “the authors interpret this result as supporting the possibility that sex is\noccuring at a low level most of the time”, we would like to clarify that this surprising finding is not\nonly supported by the Tara analyses but also by the microcosm analysis, where, at the time of\nsampling, at least four thalassiosiroid species (and more conspecific genotypes) could be\nsexualized by a salt treatment, suggesting that these populations are “ready to reproduce”, and\npart of the population may continuously undergo sexualization as they are transported into more\nsaline estuarine waters. Interestingly, while our current manuscript was under review, a very\nthorough meta-analysis of sex in natural diatom populations was published, going all the way\nback to the 19th century [3], also suggesting that sex may occur more frequently than would be\nimplied by recent mass-event observations. We now reference this review paper in the discussion.\n[1] Patil, S. et al. Identification of the meiotic toolkit in diatoms and exploration of meiosis-specific\nSPO11 and RAD51 homologs in the sexual species Pseudo-nitzschia multistriata and Seminavis\nrobusta. BMC Genomics 16, 930 (2015).\n[2] Bulankova, P. et al. Mitotic recombination between homologous chromosomes drives genomic\ndiversity in diatoms. Current Biology 1–12 (2021).\n[3] Mann, D. G. & Edlund, M. B. The Ecology of Diatom Reproduction. in Diatom Ecology 59–83\n(John Wiley & Sons, Ltd, 2024).\n\nWe appreciate the reviewer’s concerns regarding the size and completeness of the MAGs. In the\nrevised manuscript, we have addressed these points by providing a new supplementary table\n(Table S13) detailing key assembly metrics for each MAG—such as total length, N50, estimated\ncompleteness, and redundancy—alongside their taxonomic assignments and whether they\nexhibit a sexual reproduction signal. In addition, these assembly metrics were compared between\nsex and non-sex MAGs using the Mann-Whitney U test, coupled with Bonferroni correction to\nadjust p-values accounting for multiple comparisons. The results, shown in Supplementary Figure\nS30, were visualized using boxplots with significant annotations. The comparison of assembly\nmetrics for Tara MAGs between sex and non-sex samples reveals that sex MAGs are more\ncomplete, as evidenced by a statistically significant higher Anvio’ completion in the sex MAGs\ngroup (p < 0.05, Mann-Whitney U-test with Bonferroni correction). This was expected: since we\nrequire these MAGs to encode multiple sex markers, MAGs with a sexual signal are likely to be\nrelatively complete. In contrast, no significant differences were observed for the other metrics,\nincluding Anvio’ redundancy, N50, number of contigs and total length.\nAs a first step to control false positives due to vegetative expression, we have implemented\nexpression thresholds. Previously, we considered a marker to be expressed if it showed any\nmapping, in order to avoid complications associated with expression quantification in Tara (the\nmain complication being how to control for differences in relative abundance of a species between\nsamples). In the revised version, we took inspiration from the Reviewer’s suggestion to determine\n“baseline expression levels”. We used the marker’s expression levels in vegetative conditions\nfrom the experimental data to set expression thresholds above which a marker gene would be\nconsidered expressed in Tara. Specifically, we selected the 95% percentile of vegetative\nexpression for each marker family as the threshold, measured in transcripts per million (TPM). In\nother words, a marker would only be considered sexually expressed if it occurs in less than 5%\nof bulk RNA-seq samples when cells are in vegetative conditions. As shown in the figure below,\nthese thresholds are typically quite low (around 1 TPM), but, for example, Marker M3 has a higher\ncutoff value due to the aforementioned background expression in Cylindrotheca closterium.\n\nIn parallel, we calculated MAG-level TPM values for all Tara stations. Since the expression\nthresholds were determined for single (laboratory) species, we calculated Tara TPMs relative to\nthe total expression of each individual MAG in a given sample, thereby correcting for the relative\nabundance of that genotype. It must be noted that this approach does not correct for differences\nin the completeness of MAGs, although this effect is likely limited because the MAGs co-\nexpressing sex markers are overall highly complete (median value > 75%, Fig. S30).\nAfter applying the quantitative expression thresholds, the number of Tara Stations with a sex signal\ndecreases from 71 to 54. Notably, we still observe similar patterns of sex in the global ocean,\nincluding “rare” species that express sex markers, and we find similar correlations with abundance\nand seasonality. Generally, the sexual signal remained when sexual cell stages were observed in\nimaging flowcytobot data (e.g., massive Chaetoceros spermatogenesis in station 188). For\nFragilariopsis, we now only find a single sex marker expressed with SPO11-2 at Station 188, even\nthough we observed - admittedly only a single - an auxospore in that sample. We now discuss this\nevent as a possible false negative.\nAll relevant Figures, supplementary figures, tables and Results sections have been updated to\nreflect the new expression threshold approach. In addition, we have included the expression\nthresholds as a Supplementary Table and uploaded the dataset containing the TPM expression of\nall genes encoded on Diatom MAGs to the public Zenodo repository.\nIn addition, as we agree that the SPO11-2 + 2 markers cutoff is somewhat arbitrary, we now\npresent the results for different numbers of co-expressed markers (SPO11-2 + 1, 2, 3 and 4\nmarkers) in a supplementary Figure. This allows for an assessment of which sexual hits are\nsupported by more evidence, hence being more likely to be true positives, like the cases of some\nArctic Ocean stations and two tropical Pacific stations, where different genera (Chaetoceros,\nLeptocylindrus and Skeletonema in the first case, and Thalassiosira, Pseudo-nitzschia and\nMinidiscus in the latter) exhibit a strong sexual signal by co-expressing a higher number of sexual\nmarkers together with SPO11-2.\n\nIn addition to the expression thresholds explained above, we believe it is useful to consider the\nbehavior of non-sexual vegetative genes in Tara as a baseline for the level of co-expression that\ncould be expected by chance. Specifically, we now defined non-sexual “control markers” that are\nexpressed at a similar level to the sex markers in vegetative conditions, and we assessed how\noften these control markers would be considered co-expressed in Tara, by following the same\nthreshold rules as the sex markers. Since we found almost no gene families with such low\nexpression level (see below) and as these families require manual curation and phylogenetic\nselection (which is very time-consuming and not straightforward to automate), we limited\nourselves to a set of 5 control families that are as similar as possible to the pennate sex markers\nin all aspects (SPO11-2 + M1-M4). We chose to apply this approach only to the pennate markers\nand not to the centric (flagella) markers, because we only have one centric reference species and\nlack an expression atlas to accurately assess expression across a large number of vegetative\nconditions.\nApart from never being upregulated during sex, the control marker families had to be as similar\nas possible to the sex marker families: (A) encoded in the genome of a similar set/clade of\ndiatoms, (B) mostly single-copy ortholog families, and (C) having similar expression levels (both\naverage and variability) to the sex markers in vegetative conditions. To select these families, we\nused a data-driven approach similar to the one used for the pennate sex markers (see\nsupplementary figure below). Note how selecting for conserved families in step 1 already removed\nmost families showing very low expression.\n\nNext, we identified families with a similar expression level and coefficient of variation to the sex\nmarkers during vegetative conditions. However, among the 2734 families, those with an\nexpression as low as the sex markers were very rare, further supporting the specificity of our sex\nmarkers. In fact, only one family exhibited a lower average expression level than marker M1 (see\nfigure below), and further inspection revealed it is a retroviral transposon with poor conservation\nand phylogenetic descent, making it unsuitable for use. As a result, we had to select control\nmarkers with a higher baseline expression than the sex markers, which may affect the results by\nsuggesting a higher level of false positives.\n\nUltimately, after selecting control markers that are as similar as possible to M1-M4 and SPO11-\n2, we ended up with five families: C1-C4 and PC (“positive control”, similar to SPO11-2):\nWe performed HMM searches against Tara MAG genes and conducted phylogenetic selection of\npennate diatom orthologs to compare with the “pennate sex markers” M1-M4. We then considered\nhow many Tara samples were co-expressing the positive control + 2 control markers and found\nthat this occurred in seven cases (a case refers to a single MAG in a single sample). While this is\nobviously lower than what happens with the sex markers (we recorded 32 cases of pennate\ndiatoms co-expressing the sex markers and 181 cases of centric diatoms), it does show that false\npositives can occur at the natural background expression levels of sex markers.\nWe have added this control example to the manuscript and clarified that, although co-expression\nmay signal that sex is occuring (e.g., in case of the IFCB imaged sites), false positives remain\npossible. We have also added a section in the discussion where we highlight the need for some\nkind of control for false positives in observational metatranscriptome datasets like Tara. As\nsuggested by the Reviewer, we also discuss how the approaches based on expression in the lab\nrely on the key assumption that expression levels in natural species are similar to those in\nlaboratory species.\n\nTaken together, based on the Reviewer’s advice, we improved the Tara bioinformatics pipeline in\nthe following ways: (1) introducing data-driven expression thresholds, (2) visualizing the effect of\ndifferent co-expression levels, and (3) introducing a simple set of control genes to assess the\npotential of false positives.\n\nThank you for the suggestion, we have clarified the abbreviations in the top left corner of the figure\nand in the caption.\n\nWe agree with the Reviewer’s reasoning. In diatom culture experiments, there is always a specific\npart of the population that does not engage in sex (often much more than 50%). Therefore, it\nshould indeed be easier to detect upregulated genes (whose fold change should approach infinity\nin the case of a perfect marker gene) than downregulated genes (whose expression can only\ndecrease to the level of non-sexual cells in the culture). This is reflected in the volcano plots (Fig.\nS1), which are asymmetric for all species: the top upregulated genes are much more significant\nand have higher fold changes than the top downregulated ones.\nThis comment made us realize that we hadn’t sufficiently clarified this point in the manuscript.\nSpecifically, the three pennate diatoms use sex pheromones to induce a growth arrest. This\nmeans that even a part of the population not actively engaging in sex can still be affected, leading\nto downregulation of cell cycle and photosynthetic genes across the whole population [1,2]. In\ncontrast, no growth arrest and no sex pheromones affecting the entire population are known for\ncentric diatoms. This likely results in a lack of downregulated genes in centrics.\nWe have rephrased this point as follows: “S. marinoi exhibited a smaller number of downregulated\ngenes, which is likely the result of the lack of a pheromone-induced growth arrest that could cause\na general population-wide downregulation of vegetative genes even in cells not actively engaging\nin sex.”\n1. Annunziata, R. et al. Trade-off between sex and growth in diatoms: Molecular mechanisms\nand demographic implications. Science Advances 8, eabj9466 (2022).\n2. Audoor, S. et al. Transcriptional chronology reveals conserved genes involved in pennate\ndiatom sexual reproduction. Mol Ecol 33, e17320 (2024).\n\nThank you for spotting this oversight, we have corrected all references to this figure in this\nparagraph.\n\nWe agree that the statement was too subjective, based on solely qualitative observations. Unlike\nin single-species laboratory cultures, we could not quantify the level of sexual reproduction in this\ncommunity because it consisted of hard-to-distinguish Cyclotella/Thalassioriales species.\nAdditionally, we could not differentiate empty frustules that had released sperm cells from dead\ncells. The estuary system also contained a significant amount of sediment and detritus, which\nmade cell counting challenging.\nIn more qualitative terms, however, we never observed any sexual stage in the control (untreated\nmicrocosm) when screening these cultures. In contrast, after the salt treatment, a majority of\nCyclotella-sized centric cells appears to have differentiated into either sperm cells (empty\nfrustule), egg cells (surrounded by sperm cells), and later auxospores.\nAs suggested by the Reviewer, we have included original, uncropped microscopic images of\nuntreated and salt-treated cultures taken during the microcosm experiment as new\nSupplementary Figures. In addition, we have revised the text as follows:\n“In contrast to control conditions, where no sexual stages were detected during the course of the\nexperiment, we observed spermatogenesis and fertilization of egg cells after 24 h of salinity\ntreatment, followed by the development of auxospores after 48 h (Figures S8-13).”\n\nThe idea of using these markers to determine the intensity of a sexual reproduction event for each\nspecies is indeed very interesting. Specifically, the Reviewer suggests that the expression level\nof individual markers in the metatranscriptome, relative to their expression in laboratory\ntranscriptomes, can be used to estimate the proportion of a population undergoing sex. For\ninstance, if a population of 100% sexual cells expresses DRC4 at 100 CPM, then an expression\nlevel of 20 CPM in the metatranscriptome would suggest that approximately 20% of cells are (1)\nsexual cells (2) of that particular species.\nHowever, calculating such an “expression index per sexual cell” for each marker gene does not\nseem feasible with the current data for three reasons:\n1. In lab experiments, only part of the population typically engages in sex. The level of sexual\nreproduction in the “reference” RNA-seq datasets we use is variable, making it difficult to\ndetermine what expression level can be expected for each gene in a population that is\n100% sexual.\n2. The expression of sex marker genes may depend on the specific stage of sexual\nreproduction. For instance, meiotic recombination genes might show a higher expression\nlevel just before gametogenesis, followed by a burst of expression in flagella genes during\ngametogenesis. Hence, the expected expression level of each marker gene would depend\non the specific timing and mixture of sexual stages in the population in both the “reference”\nlab population and the metatranscriptome population.\n3. Different species in the metatranscriptome community may have a different efficiency in\nRNA extraction or quality, making quantitative comparisons between species unreliable.\nWhen we split the boxplots by marker to compare expression levels between species (see below,\nfor simplicity only showing the highest expressed paralog), we do not observe the expected trend\nwhere certain markers always have higher/lower expression than the others. This is possibly due\nto point (2) above: the different species may be in different (sub)stages of sexual reproduction at\nthe time of sampling.\n\nTaken together, we believe that the current box plots reflect the reality quite well: there is a general\ntrend where expression level follows species abundance, but there is considerable variability in\nthe expression level between markers, which is too big to allow for any kind of quantitative\nprediction.\nWe hope that the set of markers proposed here will stimulate future research into this topic. More\nextensive time series of sexual reproduction in nature - whether in an induced setting, like in this\npaper, or by tracking the development of a bloom or sex event - would make it possible to link sex\nintensity to expression levels and connect expression timing with the succession of developmental\nstages.\n\nWe have toned down these statements, and now also acknowledge that co-expression of sex\nmarkers is only an indicator of sex. Some examples below.\n“In contrast, the co-expression of sex marker genes by locally rare species seems to suggest that\nother genera were sexually active at very low abundances: (...)”\n“Although the co-expression of sex markers is only an indirect proxy for sexual activity, the\nidentification of a sexual signal across the majority of Tara Oceans stations worldwide suggests\nthat diatom sexual events may occur more frequently in the global ocean than previously thought.”"}, {"round": 2, "reviewer_comments": {"reviewer_2": "I congratulate the authors for their work. The revised manuscript is even more rigorous and I think it represents a\nbreakthrough both for diatom biology and for pelagic microbial oceanography more broadly.\nThey reached my highest expectations in responding to my suggestions for trying to use \"control\" genes to assess potential\nfalse positives. I expect the bioinformatic work to satisfy this suggestion was not trivial. I point out that this type of approach\nwill be offer a robust strategy to investigators seeking to use metatranscriptomic and MAG data for detecting other types of\nevents that are challenging to observe (e.g., biological interactions, rare use of alternative nutrient or energy sources).\nI also appreciated that they responded to my suggestions by adding some phrases to speculate succinctly on the\nimplications of their observations for how planktonic microbes might interact.\nI have some very minor comments the authors can correct on their own as they see fit, without any need for a re-review on\nmy part:\np. 9, line 49 - space missing \"ofsperm\"\np. 16, line 248 \"Widespread diatom sexual reproduction in the global ocean To detect sexual reproduction in the\" Seems\nlike a section title is not properly formatted\np. 16, lines 253-255: \"Therefore, we only considered marker genes sexually expressed above a threshold that corresponds\nto the expression in maximum 5% of laboratory RNA-seq samples in vegetative conditions.\" The sentence needs\nclarification. What is meant by \"in maximum 5%\"? Does it mean it is only upregulated in a maximum of 5% of all vegetative\nconditions tried? Or that it was only detected in 5% of vegetative samples?", "reviewer_3": "As a non-specialist in diatom life cycle, my comments mainly concern the methodology used to study the metatranscriptome\ndatasets. In line with the comments of the first 2 reviewers, I found the work done to select diatom marker genes specific to\nsexual reproduction very good. I think that this gene set could be useful for a rapid detection of sexual events in natural\npopulations. I am also quite convinced by their conclusions that sexual reproduction occurs in many different oceanic\nregions. My main concern is the important lack of normalisation procedure (sorry if I missed this in the Mat&Met) to draw\nconclusions about the quantification/prevalence of these reproductive events in diatom populations (both microcosm and\nTara datasets). Details for each analysis are given below.\nMy main concern with the microcosm experiment is the way the differential expression was calculated. If I understand the\nmethodology correctly, differential expression levels were calculated for all genes in the community and not for each taxa\nindependently. Consequently, observed expression differences could be due to variations in cell number rather than gene\nregulation. As mentioned lines 209-210, most genes are downregulated because due to the mortality of non-diatom taxa. So,\nI think all diatom genes will appear upregulated (sex marker genes and the others) because their relative abundance has\nincreased in the community. Could the authors show that sex marker genes are more upregulated than other diatom genes ?\nI think this is necessary to conclude that the salt treatment induced the expression of sex genes.\nI have seen in the revised version that the authors have added 3 vegetative genes as a control. Is it possible to add the\nexpression of these genes in Figure 4e to show that not all diatoms genes are upregulated ? Another possibility could be to\nuse all Cyclotella transcripts as a baseline and see if Cyclotella sex gene markers represent a higher proportion of its\ntranscripts in the salt condition (and do the same for each diatom genus).\nI don’t understand what is the conclusion of Figure 4f. As Thalassiosira is less abundant than Cyclotella, all Thalassiosira\ngenes will have a lower read count. I do not see anything specific to the sexual gene markers here (except if I missed\nsomething in the normalization method).\nThalassiosira and Conticribra do not appear in the 18S barblot (Figure 4c) because of their low abundance. Is it possible to\nadd a panel showing the variation of the 18S relative abundance of the four clades highlighted in panels e and f ?\nThe methodology employ to find sex genes in Tara Oceans MAGs is good and I’m convince by the expression of diatom sex\ngenes in many oceanic regions (Figure 5). The use of imaging data is a nice complement to prove the presence of sexual\nevents.\nHowever, I’m not convinced by the conclusion that “there is more abundance of Chaetoceros in stations where there is a\nsexual activity” (lines 302-306). How can we exclude the possibility that sex genes are detected only because Chaetoceros\nis abundant enough to exceed the sequencing depth threshold? Are there oceanic stations with high abundance of\nChaetoceros without the expression of sex genes ? If so, this could be very interesting and should be presented. Otherwise,\nwe could interpret the results as sex genes are expressed at low level everywhere and only when the species is abundant\nenough we can detect it within metatranscriptome data.\nIn the discussion (lines 368-369), to say that sexual reproduction is more prevalent in some stations, I think the expression of\nsexual genes needs to be normalized to the expression of all other genes for each diatom species. Same problem for\nThalassiosira in station 123 (lines 380-382).\nIs it possible to quantify the proportion of cell ongoing sexual reproduction with imaging datasets ? I think this would be a\ngood indication of a higher prevalence in a sample.\nOther points :\nFigure 1c : Statistical tests on Pearson correlations are missing.\nLine 654 : Could you specify which metatranscriptomic samples you used for the assembly ? I assume a co-assembly of the\n6 samples ?\nLine 675-676 : Did you apply EdgeR differential expression directly to the 850,000 transcripts ? I would be curious to see the\nMA plot or violin plot (coloured by taxa).\nFigure 4e (and related SupFigures) : Could you explain how this normalisation between 0 and 1 was done ? Is it a simple\nratio per gene of the CPM values ? I think it is also important to see the CPM counts here (or at least the average per line).\nLines 721 to 728 : These phylogenies used to select marker genes within Tara MAGs should be available somewhere to\nunderstand how the MAG proteins were selected according to the phylogenetic placement (as SupFigure or in a public\nrepository).\nLine 359-360 : “more frequently in the global ocean“ is not clear for me. Do you mean “higher frequency of sexual events in\nnatural diatom populations than previously though” ? I think this has not been shown in the manuscript. Instead, you could\nwrite “sexual events occur in more oceanic regions than previously though”."}, "author_response": "We thank the Reviewer for their supportive comments and for re-examining our manuscript in\nsuch detail. We are happy that our improvements to the previous version of the manuscript were\nsatisfactory.\nWe have implemented the remaining improvements suggested by the Reviewer below.\n\nWe thank the Reviewer for critically assessing the methodology of our metatranscriptome\nanalyses. We indeed previously did not perform normalization for the microcosm differential\nexpression experiment, and have now implemented this in the revised manuscript. For the Tara\ndataset we did include a correction for taxon-abundance but this was not explained clearly\nenough in the manuscript. Below, we provide a point-by-point response to your concerns.\n\nIt is correct that we did not perform any kind of normalization that\ntakes into account the taxon-specific abundance shifts between\ntreatments in the microcosm experiment. However, we reasoned\nthat, despite some shifts in community composition (most\nstrikingly the demise of Leptophryidae parasites after salt\ntreatment), the changes in the relative abundance of\nThalassiosirales diatoms were not sufficient to fully explain the\nobserved upregulation of sex marker genes. As can be seen from\nthe 18S rDNA barcoding barplots on the right, Cyclotella\nincreased in relative abundance by about 40%, but some of the\nother Thalassiosirales clades even decreased after a salt\ntreatment, a pattern inconsistent with the concurrent upregulation\nof sex marker genes.\n\nWhile we previously already annotated raw microcosm RNA reads and transcripts with Kraken2,\nwe needed a better taxonomic resolution to answer your questions, especially for the\nThalassiosirales. Therefore, we created a new Kraken2 protein database that includes 97 new\ndiatom proteomes, among which 87 thalassiosiroids. Apart from offering a higher taxonomic\nresolution, this also led to a better capture of reads: the number of unassigned reads decreased\nby ~10%. The Methods section and Supplementary Figures 15 and 16 were updated to reflect\nthis change in the Kraken database.\nWe next applied Kraken2 to the transcripts from the microcosm de novo transcriptome. To verify\nthe accuracy of the extended Kraken database, we checked the sex marker gene hits, which we\npreviously manually annotated using phylogenetic analysis. Kraken2 correctly assigned all of\nthese sex marker hits to the order Thalassiosirales and 26/27 transcripts were correctly placed\namong the four Thalassiosira clades that we use in the paper. Hence, these Kraken2 taxonomy\npredictions are accurate and can be used to determine whether expression trends are affected\nby community composition.\nWe first compared the differential expression results\nbetween the four most common taxa (Bacteria, Diatoms,\nFungi+Animals, and plants) (see barplots). Surprisingly,\ndiatom transcripts were less likely to be upregulated in\nresponse to salt compared to the other taxa (11.5% for\ndiatoms vs 15.1% on average for other taxa). The\nmajority of the sex marker hits (21/27) was among the\n10% most upregulated Thalassiosirales transcripts in the\ndataset, further supporting that their response is not just\na side-effect of a general shift in abundance.\n\nBecause changes in abundance may\nonly affect the differential expression at\na lower taxonomic level (e.g. the\nCyclotella species that increased in\nabundance), we visualized the fold\nchanges for each of the four\nThalassiosirales clades separately. As\ncan be seen from the volcano plots on\nthe right side, Cyclotella is indeed\nslightly skewed towards upregulated\ngenes, but the majority of Cyclotella\ngenes (about 80%) are not\ndifferentially expressed.\n\nTaken together, these results indicate that the upregulation of diatom sex markers is not a\ngeneral consequence of increased diatom abundance. However, we agree that normalization by\ngenus-level abundance would be the most rigorous approach. We have therefore repeated all\nanalyses relative to libraries with all genus-level reads (see below).\n\nAs suggested by the Reviewer, we subdivided the microcosm RNA-seq dataset into the four\nThalassiosirales clades (Cyclotella, Thalassiosira, Conticribra and Skeletonema) using the new\nKraken2 taxonomy. By performing a separate differential expression analysis for each of these\nclades, we corrected for differences in library size (total number of reads captured) due to\nchanges in relative abundance for that specific taxon. The results of these statistical tests were\nnearly identical to the DE analysis over all taxa we performed earlier (26/27 sex markers\nidentical).\nWe also recalculated the counts per million expression values for each clade separately, and\nthis also did not change the outcome (see comparison below). We have updated the\nmanuscript, figures and Zenodo datasets with these new per-clade normalized DE analysis and\nfigures.\n\nWe agree that this panel was mainly showing differences in abundance of the individual species\nand so is not really informative. We have removed it in favor of the barplot of 18S barcoding\n(suggested below).\n\nWe added barplots showing the 18S relative abundance of the four main Thalassiosirales\nclades during the microcosm experiment as a new Figure 4d.\n\nWe thank the reviewer for this insightful comment. We agree that it is important to distinguish\nbetween biological regulation and detection limits driven by species abundance. First, we would\nlike to clarify that we apply a normalization strategy that mitigates the confounding effect of\nspecies abundance. Specifically, we compute MAG-level TPM values, meaning that the\nexpression of each gene is normalized to the total expression from the same MAG. This\nwithin-MAG normalization ensures that genes whose expression scales linearly with MAG\nabundance do not increase in TPM as the species become more abundant (e.g. Chaetoceros).\nTherefore, our fixed TPM expression thresholds for detecting sexual reproduction (based on\nexperimental data) are independent of whether a species is rare or abundant in a sample.\nOf course, as suggested by the Reviewer, it is possible that sexual reproduction is more\nwidespread than detected, and that sex markers are missed in some samples due to insufficient\nsequencing depth. However, we note that we have successfully detected sex marker expression\nin several diatom MAGs that are locally very rare (e.g., <1% of mapped reads), suggesting that\nsex signals can be detected even when abundance is low.\n\nAs for Chaetoceros: this genus indeed shows significantly higher abundance in stations where\nsexual reproduction signals are detected, as shown in our statistical analysis (lines 302-306).\nHowever, we have also observed several stations with high Chaetoceros abundance but no\nevidence of sex marker expression (see Supplementary Figure 28 and additional bar plots\nprovided below). This suggests that high abundance alone is not sufficient to explain the\npresence of sex marker expression.\nTo acknowledge the valid concern raised by the reviewer, we now include a sentence in the\nmain text to acknowledge that abundance-driven detection cannot be fully excluded, but it is\nunlikely to be the only explanation, as sexual reproduction signals are not detected in all\nhigh-abundance stations.\n\nWe thank the reviewer for this thoughtful comment. We agree that normalization is critical to\ndraw meaningful conclusions about gene expression levels. As correctly pointed out, the\nexpression of sex marker genes must be interpreted in the context of the overall transcriptional\nactivity of the corresponding diatom species.\nAs explained above, we calculated MAG-level TPM values, meaning that each gene’s\nexpression was normalized by the total number of reads mapped to the same MAG in each\nsample. This ensures that the expression of sex marker genes is normalized relative to the\noverall transcriptional activity of the same MAG, effectively correcting for both transcript\nabundance and the underlying MAG abundance across samples. This normalization allows for\nrobust comparisons across stations.\nTo clarify this methodological point, we revised the manuscript as follows:\nRevised text:\n“By calculating MAG-level TPM values, we correct for the relative abundance of each MAG\nacross samples, ensuring that the expression of sex marker genes is normalized relative to the\noverall transcriptional activity of the same MAG. This enables consistent comparisons across\nstations and provides a robust basis for applying gene-specific expression thresholds\ndetermined for single species in a laboratory setting.“\n\nQuantification of sexual and vegetative cell stages in the IFCB/EcoTaxa data is not possible\nbecause of limitations of this dataset, in particular the variability in data quality between\nsamples. Some stations have thousands of pictures, others only have a few or captured just air\nbubbles and no cells.\nHowever, we previously performed a qualitative survey of the IFCB dataset: for each station, we\nrecorded whether we observed vegetative cells and/or sexual cell stages for each genus, and\ncross-referenced this with the sex marker signal. This analysis also supports that sex is not\nalways occurring when a genus is abundant: we often identified vegetative cells of genera that\ndid not display a sex signal (metatranscriptome) and sexual cell stages (IFCB).\nThis dataset can be found in our Zenodo repository, or through this link:\nhttps://zenodo.org/records/14963127/files/Ecotaxa_IFCB_observation_sheet.xlsx?download=1 .\nWe have more clearly referred to this in the context of Chaetoceros abundance.\n\nThank you for spotting this. We added a matrix with adjusted p-values of statistical testing on\nPearson correlations to Fig. S1 and discuss the results in the main text.\n\nYes, that is correct, we have now clarified this in the Methods section.\n\nOriginally, we had indeed applied EdgeR to the entire dataset of >800.000 transcripts. You can\nfind volcano plots from this statistical test, separated for the four taxa that we focus on, pasted\nat the beginning of this rebuttal.\nNote that based on your feedback, we have now changed the procedure and performed\ndifferential expression analyses separately for each Thalassiosirales clade (hence taking into\naccount differences in library size due to relative abundance in each sample).\n\nIndeed, the “normalized CPM” we use for heatmaps are CPMs that are scaled to a maximum of\none per gene. This is now explained more explicitly in the captions.\nWe used gene-level scaling because the genes with the highest expression during sex made it\nimpossible to see the expression dynamics for other genes in the heatmap. E.g. Sig1 in\nConticribra reaches about 2500 CPM during sex (0.25% of the entire transcriptome) which\nmakes other genes nearly invisible when using a single non-normalized scale. What’s more, we\noften identified multiple isoforms (or duplicates) of the same marker in the same taxon, and\nwhile all of them were upregulated after the salinity shift, their baseline expression differed\nsubstantially. Because we wanted to visualize that each marker hit (transcript, isoform) in itself is\nspecific for sex, we used the normalization so that also the lower expressed ones are visible.\nWe agree that it is more fair to also show the raw (non-scaled) CPM. Because the heatmaps\ncould not clearly display this, we made separate barplots with raw CPMs for all marker hits in\nFig. 4f, and included them as a new Supplementary Figure (S30).\n\nWe completely agree that the phylogenetic selection of orthologs of sex markers should be\ntransparent. To this end, the phylogenetic trees are available on the public Zenodo repository,\nboth for the microcosm and the Tara Oceans phylogenetic analyses. In fact, these trees were\nalready available on Zenodo for the previous versions of this manuscript\n(https://zenodo.org/records/14963127, under “Phylogenetic_selection.zip”). However, this wasn’t\nclearly stated in the data availability section. We now clarified this.\n\nWe thank the reviewer for this helpful comment. We agree that our original wording may have\nimplied a quantitative increase in the frequency of sexual events, which we do not directly\nmeasure. Our intended message is that because diatom sex is traditionally considered rare in\nnature, the detection of co-expression of sex marker genes across more than 50 Tara Oceans\nstations is already noteworthy. While we do not quantify frequency, the widespread geographic\noccurrence of this signal suggests that sexual reproduction may be more common than\npreviously assumed.\nTo reflect this more clearly, we have revised the sentence as follows:\nRevised text (lines 358-360):\nAlthough the co-expression of sex markers only serves as an indirect proxy for sexual activity,\nthe identification of a sexual signal across more than 50 Tara Oceans stations worldwide is\nnotable, given the long-standing view that diatom sex is a rare event in nature\nTo ensure consistency, we also revised the corresponding sense in the abstract:\nRevised abstract sentence:\n“Our results reveal that diatom sexual reproduction is more widespread in the global ocean than\npreviously thought, encompassing both dominant bloom-forming species and rare taxa.”"}]

Decision letter