Prior to cell culture, pyrosequencing of the 18S rRNA (3′ end), internal transcribed spacer 1 (ITS1), 5.8S rRNA, ITS2, and large subunit (LSU) rRNA genes was performed to identify the Alexandrium strain with the following paired primers: (i) ITS_F (GAGGAAGGAGAAGTCGTAACAAGG) and R732 (CCTTGGTCCGTGTTTCAAGAC), (ii) F341 (AGCGCA CAAGTACCATGAGG) and R1344 (TTCGGCAGGTGAG TTGTTACACAC), (iii) F918 (AGGTTCGTATCGATACTGACG TGC) and R2242 (CAGAGCACTGGGCAGAAATCAC), (iv) F1976 (GAAAAGGATTGGCTCTGAGG) and R3016 (AAACTAACCTTCTCACGACGGTC), and (v) F2843 (TTAC CACAGGGATAACTGGCTTG) and ETS_R (GCAGGGACA TAATCAGCACAC). Further, BLASTN analysis was performed using the pyrosequencing results (Supplementary Table 1).

A. catenella (Group I) collected from the South Sea in South Korea, was isolated using the micropipetting method, and maintained in sterilized Guillard’s (F/2) medium (Sigma, St. Louis, MO, United States) with Red sea salt of 31 PSU (Red Sea, Houston, TX, United States), and subcultured with a 1/10 dilution approximately 40 times under a light–dark cycle (12 h each) with 50 μmol⋅m^–2⋅s^–1 photons at 20°C ± 1°C. Pseudoruegeria sp. M32A2M was isolated from surface water of the South Sea in South Korea, where HABs of Alexandrium are frequently observed (Cho et al., 2020). Pseudoruegeria sp. was primarily cultured on marine agar 2216 plates (BD, Franklin Lakes, NJ, United States) at 25°C for 2 days. The cells scraped from agar plates were resuspended using gentle pipetting in F/2 media with Red sea salt and incubated for 2 days at 20°C ± 1°C for adaptation in liquid media prior to co-culture with A. catenella. Thereafter, the Pseudoruegeria sp. culture medium was filtered through a 20 μm nylon net by gravity to remove some aggregated cells. Then, the OD_600nm of filtrated cells was measured before inoculation into A. catenella cell medium. For this study, a total of 10 L of A. catenella (2,000 cells⋅mL^–1) and 1 L of Pseudoruegeria sp. (OD 0.3 at 600 nm) were cultured in F/2 medium.

For the co-culture of the two microorganisms, 100 mL of Pseudoruegeria sp. culture (OD 0.3 at 600 nm) estimated to be 6 × 10¹⁰ cells in total, were each loaded into nine flasks containing 1 L of A. catenella culture (2,000 cells⋅mL^–1) of 2 × 10⁶ cells in total. Three sets of flasks with three biological replicates were exposed to light for 2, 6, and 12 h with 50 μmol⋅m^–2⋅s^–1 photons at 20°C ± 1°C, respectively. As control groups, 100 mL of fresh medium was added into 1 L of pure-cultured A. catenella (2,000 cells⋅mL^–1), and 1 L of fresh medium was added into 100 mL of Pseudoruegeria sp. medium. Then, diluted cells were collected again using centrifugation prior to sequencing preparation.

Chlorophyll-a (Chl-a) and carotenoid were extracted from A. catenella with 90% ethanol overnight at 4°C. These compounds were measured spectrophotometrically as follows: Chl-a (mg⋅L^–1) = (12.7 × A₆₆₅) - (2.69 × A₆₄₅), and carotenoid (mg⋅L^–1) = (1,000 × A₄₇₀ - 2.05 × C_Chl–a)/245 (where A₆₆₅, A₆₄₅, and A₄₇₀ indicate the absorbance at wavelengths of 665, 645, and 470 nm, respectively, and C_Chl–a represents the Chl-a content) (Li et al., 2015a).

The co-cultured A. catenella and Pseudoruegeria sp. mixtures were separated through 20 μm nylon net filter membranes (Merck Millipore, Burlington, MA, United States) by vacuum to avoid cross contamination of mRNA and to further construct clearer sequencing libraries. Then, the membrane was washed with 200 mL of fresh F/2 medium. The A. catenella collected on a filter membrane was resuspended with 50 mL of F/2 medium and centrifuged at 1,503 × g at 4°C for 10 min. The bacterial cells filtrated through the membrane were centrifuged at a high g-force of 3,381 × g for 10 min to collect the bacterial cells. The centrifugation at a high g-force for bacterial cells also helped disrupt the fragile cell membrane of a few residual A. catenella. A. catenella and bacterial cell pellets were resuspended with RNAlater^TM (Invitrogen, Carlsbad, CA, United States) and incubated at 4°C for 1 h. Each pellet was stored at -80°C until further experiments.

For total RNA extraction, cells were washed with extraction buffer comprising 200 mM Tris–HCl (pH 7.5), 25 mM EDTA, 250 mM NaCl, and 0.5% SDS. Then, each pellet was instantly frozen using liquid nitrogen and ground to a fine powder. For RNA extraction of A. catenella and Pseudoruegeria sp., TRIzol (Invitrogen) and phenol (Sigma) at 65°C were added to the powdered pellet, and this was immediately followed by ethanol precipitation.

The RNA-seq library was constructed using the TruSeq mRNA Sample Prep Kit (Illumina, San Diego, CA, United States) according to the manufacturer’s instructions. Briefly, the poly-A-tailed mRNA of A. catenella was purified with oligo-dT magnetic beads. High-throughput sequencing with 100 bp single-end reads technically requires at least 100 bp RNA fragments. However, it is important to mention that RNA fragmentation causes a broad length distribution. To account for this, the purified mRNAs were further fragmented for 2 min to obtain RNA fragments with a median value of about 150 bp. cDNAs synthesized from the fragmented RNAs were ligated to adaptors and then amplified using qPCR.

For bacterial mRNA preparation, rRNAs from purified total RNAs (4 μg) were first depleted using the Ribo-zero^TM rRNA Removal Kit for Epidemiology according to the manufacturer’s instructions (Illumina). After ethanol precipitation of mRNAs (∼0.4 μg), the rRNA removal was confirmed using a TapeStation system 2200 (Agilent, Santa Clara, CA, United States). Then, the bacterial mRNAs were also fragmented to lengths of a median value of about 150 bp using a TruSeq mRNA Sample Prep Kit (Illumina), followed by ethanol precipitation.

The total library length for A. catenella and Pseudoruegeria sp. was approximately 250–300 bp by adaptor ligation and overlap extension qPCR for sequencing. The resulting library was sequenced by a HiSeq 2500 instrument with 100 bp single-end reads (Illumina).

After RNA-seq raw reads of A. catenella and Pseudoruegeria sp. were obtained from Illumina sequencing, the adapter sequences and low-quality sequences were trimmed using IlluQC software included in the NGSQCToolkit (Patel and Jain, 2012) and CLC Genomics Workbench 6.5.1 (CLCbio, Denmark) with the default setting. Then, de novo transcriptome assembly from RNA-seq raw reads of pure-cultured A. catenella was performed using Trinity v 2.0.6, with default settings for strand-specific RNA reads (Grabherr et al., 2011).

The RNA-seq reads from A. catenella were then aligned to the Trinity-assembled transcriptome by Bowtie2 to remove the misassembled sequences and their abundance was calculated using RNA-seq by Expectation Maximization (RSEM) (Li and Dewey, 2011; Langmead and Salzberg, 2012). RSEM results were transformed to fragments per kilobase of transcript per million (FPKM) values, and mis-assembled transcripts (FPKM cut-off = 1) were filtered out (Trapnell et al., 2010). Further, redundant isoform unigenes were removed by clustering with the BLAST-like alignment tool (BLAT) using tileSize = 8 and stepSize = 5, with other parameters set to default (Kent, 2002; Zhang et al., 2014b). Because this process results in a non-redundant set of the longest representative transcripts (i.e., unigenes), the analysis for various RNA isoforms by the alternative splicing or alternative polyadenylation was excluded.

From the unigene reference, the coding regions of the assembled transcript sequences were predicted using TransDecoder¹ with the assistance of UniProt, BLASTP, Pfam, and HMMscan (HMMER v.3.2.1) (Haas et al., 2013; El-Gebali et al., 2018). For function annotation from coding sequences, BLASTx alignment was performed against the UniProt database with an E-value threshold of 10^–6 (Camacho et al., 2009; UniProt: a hub for protein information, 2014). In addition, Gene Ontology (GO) functional annotation was performed based on the annotated UniProt ID (Gene Ontology Consortium, 2014). Together, Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog ID and EuKaryotic Orthologous Groups (KOG) categorization were conducted (Koonin et al., 2004; Kanehisa et al., 2015).

Raw RNA reads from Pseudoruegeria sp. were mapped against the Pseudoruegeria sp. M32A2M reference genome in DDBJ/ENA/GenBank (GenBank assembly accession, GCA_010374725.1) using the CLC Genomics Workbench (version 6.5.1) (similarity score = 0.9 and length fraction = 0.9). GFF files were visualized using SignalMap software v.2.0.0.5 (Roche, Pleasanton, CA). Further, KEGG Ortholog ID and Clusters of Orthologous Groups (COG) categorization was performed based on the genome sequence of Pseudoruegeria sp. (Tatusov et al., 2000; Kanehisa et al., 2015).

For quantitative analysis of RNA expression, the expression levels of the unigenes of A. catenella and transcripts of Pseudoruegeria sp. were normalized by DESeq2 and calculated as reads per million (Love et al., 2014). Differentially expressed genes (DEGs) were classified when the fold change in the normalized RNA expression value of the co-cultured cells against that of the pure-cultured cells was more than twofold (|log₂Fold change (FC)| ≥ 1) with a low p-value (DESeq P < 0.05). Transcripts with low expression levels < 50 reads per million in the pure-cultured cells were excluded for DEG analysis to avoid over-interpretation. GO enrichment analysis with DEGs was conducted using the BiNGO Cytoscape application (Maere et al., 2005).

First, we identified the Alexandrium strain to be used in this study through the pyrosequencing of the LSU rRNA gene and the 18S rRNA fragment (Supplementary Table 1). We found that the Alexandrium strain was highly matched with the top four A. catenella strains screened in South Korea (DQ785887.1, DQ785886.1, AY347308.2, and DQ785885.1) in terms of total score from BLASTN analysis (Supplementary Table 1). In a previous study, we screened the marine bacterium Pseudoruegeria sp. M32A2M from regions where A. catenella was frequently found, and we characterized its genomic information through 16S rRNA sequencing and genome sequencing (Cho et al., 2020).

To investigate the algicidal activity of Pseudoruegeria sp., the Chl-a and carotenoid concentrations in A. catenella were measured according to the time of co-culture with Pseudoruegeria sp. (Figure 1A). The Chl-a and carotenoid levels started to decrease from 2 h of co-culture and decreased significantly at 12 h of co-culture to 38.7 and 28.5% of the corresponding levels in pure-cultured cells, respectively. In pure culture, A. catenella had an integrated membrane with four chains; however, the cells began to remarkably lose membrane integrity following 24 h of co-culture (Figure 1B). A. catenella has highly active flagella, as indicated by the phylum name Dinoflagellata. A. catenella showed highly active motility in pure culture (Supplementary Video 1), but after 24 h of co-culture, its motility was greatly decreased (Supplementary Video 2).

To understand how A. catenella and Pseudoruegeria sp. communicate with each other during co-culture and how their cellular functions have been altered, we investigated the changes in the cellular transcriptomes in both organisms with the passage of co-culture time. As a control, pure-cultured A. catenella and Pseudoruegeria sp. were used for subsequent comparison with co-cultured cells (Figure 2A). For co-culture, the two microorganisms were incubated together for 2, 6, and 12 h with light exposure. Mixed cells with a large difference in size were isolated using a net membrane filter (Figure 2A), because A. catenella and Pseudoruegeria sp. lengths are 20–48 μm and 1–5 μm, respectively (Fukuyo, 1985; Cho et al., 2020). From the two isolated organisms, mRNAs for eukaryotes and bacteria were individually purified and sequenced using high-throughput sequencing (Figure 2A and Supplementary Table 2). Consequently, the total mapped reads from A. catenella in high-throughput sequencing were 252,490,204 (Figure 2B). The genome sequencing of A. catenella, which has a large genome size of more than hundreds of giga base pairs, is challenging. Therefore, to create the reference for A. catenella, the RNA-seq raw reads from pure-cultured A. catenella were de novo assembled using Trinity V 2.0.6. The total number of assembled transcripts was 245,489 (Figure 2B). After the exclusion of isoforms by BLAT, 173,134 unigenes were obtained. Among the unigenes, 39,759 genes were functionally annotated through the functional retrieval procedures based on UniProt, GO, KEGG, and KOG (Supplementary Dataset 1 and Figure 2B).

The transcriptomes for pure-cultured and co-cultured (2, 6, and 12 h) Pseudoruegeria sp. were analyzed for the reference contigs (4,979 genes) (Cho et al., 2020; Supplementary Dataset 2). To verify experimental reproducibility between three biological replicates and to observe transcriptome changes in different culture conditions, we calculated Pearson’s correlation coefficient for the unigene reference genes using the mapped read counts normalized by DESeq2. As a result, both A. catenella and Pseudoruegeria sp. showed high Pearson’s correlation coefficient of over 0.99 between three biological replicates (Figure 2C). When observed between incubation times, A. catenella showed slight transcriptome variations for 2 h and 6 h but significant changes after 12 h (Figure 2C). On the other hand, the transcriptome of Pseudoruegeria sp. started to change significantly from 2 h of co-culture. Similarly, in the principal component analysis, A. catenella and Pseudoruegeria sp. showed distinct transcriptome changes according to the co-culture time with high reproducibility (Figure 2D). Taken together, the cellular transcriptomes of A. catenella and Pseudoruegeria sp. affected each other and showed drastic changes at 12 h of co-culture.

To investigate the reciprocal responding genes by co-culture, we selected the DEGs that satisfied the transcriptional changes over twofold with low p-values (DESeq P < 0.05) after co-culture (see volcano plots in Supplementary Figure 1). Consequently, for A. catenella, among 39,759 genes, 1,716 (4.3%) and 2,602 (6.5%) genes were found to be upregulated and downregulated, respectively (Figure 3A). In particular, 120 genes were commonly downregulated at 2, 6, and 12 h. For Pseudoruegeria sp., more than 20% of 4,246 genes were differentially regulated, i.e., 904 (21.3%) and 946 (22.3%) were upregulated and downregulated, respectively (Figure 3B).

The transcriptomes of Pseudoruegeria sp. and A. catenella also remarkably changed at 12 h of co-culture. To examine the patterns of transcriptome changes by co-culture time, all DEGs of both microorganisms were classified into 12 groups according to their fold change and expression pattern (Supplementary Datasets 1, 2) and were visualized using heatmaps (Figures 3C,D). For example, the 12th group of A. catenella, which showed a highly increased expression pattern only at 12 h of co-culture, consisted of transcripts for three cell wall-associated hydrolases (^∗), which may cause membrane destruction; a blue light receptor (^∗∗); and uncharacterized proteins (^∗∗∗). For Pseudoruegeria sp., the tripartite tricarboxylated transporter clusters (tctC-tctB-tctA) and DUF2282 domain containing predicted integral membrane proteins were assigned to the 12th group, which showed rapid increase at 12 h of co-culture, suggesting that the uptake of short carbohydrates can be active at 12 h of co-culture (Figure 3D).

To systematically investigate the cellular functions of genes changed by the co-culture of the two organisms, the DEGs in A. catenella were classified through GO enrichment analysis (Supplementary Table 3). Interestingly, the genes involved in the photosynthetic electron transport chain and PSII started to be downregulated from 2 h (Cluster 1/1) (Figure 4A). At 12 h, downregulation expanded to genes related to carbohydrate metabolism (Cluster 1/13); phosphorus/phosphate metabolic processes (Cluster 2/13); and the detection of external stimuli, including light- and abiotic factor-related genes (Cluster 3/13), as well as photosynthesis. In addition, the genes related to the movement of microtubules, which are associated with cell motility and division, as a part of the cytoskeleton (Vale, 2003), were jointly downregulated (Supplementary Table 3). These results are consistent with the drastic decrease in Chl-a and carotenoids, which are important for photosynthesis, and reduced motility and membrane disintegration of A. catenella. For upregulated clusters, the transcription of macromolecule biosynthetic process/metabolic pathways, ribosome assembly, and translation-associated genes increased from 2 h (Supplementary Table 3). At 12 h, the transcription levels of RNA synthesis/processing/modification and ribonucleotide/nucleoside-related processes markedly increased.

Furthermore, the functional changes from the perspective of metabolic pathway were investigated through KEGG analysis. As the reference genome for KEGG analysis, the first fully genome-sequenced eukaryotic marine phytoplankton, Thalassiosira pseudonana, was used (Armbrust et al., 2004). Overall, A. catenella expressed a drastically changed transcriptome at 12 h (Figure 4B). However, energy metabolism-associated genes, including photosynthesis and oxidative phosphorylation, uniquely began to be downregulated continuously from 2 h, and the number of downregulated genes gradually increased with co-culture time (Figure 4B and Supplementary Figure 2A). At 12 h, energy metabolism genes associated with carbon fixation, nitrogen metabolism, and sulfur metabolism were additionally downregulated (Figure 4B and Supplementary Figure 2A). Furthermore, carbohydrate metabolism genes associated with glycolysis, the tricarboxylic acid (TCA) cycle, and the pentose phosphate pathway were intensively downregulated at 12 h (Figure 4B and Supplementary Figure 3).

The results of GO and KEGG analysis together show that the transcription of genes associated with photosynthesis and oxidative phosphorylation in A. catenella is primarily deactivated by co-culture with Pseudoruegeria sp., which seems to significantly affect carbon fixation- and carbohydrate-associated metabolism upon long-term exposure. In detail, we investigated the transcriptomic changes in photosynthesis-related oxidative phosphorylation among metabolic pathways at the gene level (Figure 4C). Among photosynthesis pathways, PSII is the first light-dependent oxygenic protein complex (Supplementary Figure 3). It consists of D1 protein (PsbA), which binds chlorophyll P680, and D2 proteins (PsbD) as reaction-center proteins. It is also surrounded by a CP43/CP47 (PsbB/PsbC) core complex of antenna and light-harvesting complexes to employ chlorophyll and pigments to absorb light. It also forms a complex with Cyt b559 subunits (PsbE/PsbF), oxygen-evolving complex (PsbO and PsbP), and other complex proteins. The transcription levels of psbA, psbB, psbC, psbD, and psbE in the PSII complex were very high, ranking as 1st, 2nd, 13th, 22nd, and 30th, respectively, in the overall transcript expression profile (Supplementary Dataset 1); however, they were all significantly downregulated upon short-term exposure to Pseudoruegeria sp. (Figure 4C). The core petB (cyt b₆) and petD (cyt b₆f-IV) subunits of the cytochrome b₆f complex that transfers electrons from PSII to PSI (Kurisu et al., 2003) were gradually downregulated with an increase in co-culture time (Figure 4C and Supplementary Figure 3). Consecutively, similar to PSII, PSI reduced the transition from NADP to NADPH through light-driven reactions of photosynthesis. PsaA (A1) and PsaB (A2), which are photosynthetic reaction center proteins among the PSI complex, were also downregulated upon co-culture (Figure 4C). This is consistent with the downregulation pattern of PSII center proteins (D1 and D2 proteins). The downregulation of the center proteins in the two photosystems can inhibit the assembly of the complex, including the antenna proteins and pigments. In addition, the transcription of atpA and atpB of F1, which is responsible for hydrolyzing ATP in the ATP synthase complex, was highly downregulated (Figure 4C). ATP from ATP synthase and NADPH from PSI are necessary for carbon fixation. As a result, the downregulation of photosynthesis-associated core genes consequently led to the downregulation of carbon fixation metabolic processes, including rubisco, glycolysis, and pentose phosphate pathways (Supplementary Figure 3).

The DEGs in Pseudoruegeria sp. were functionally classified through GO enrichment analysis (Supplementary Table 4). In particular, translation-associated genes (Cluster 1/3) and electron transport chain genes (Cluster 2/3) were continuously upregulated from 2 h onward (Figure 5A). Unlike the observations in A. catenella, ATP biosynthetic process genes (Cluster 1/7) and energy generation by respiration and TCA genes (Cluster 2/7) were highly upregulated at 12 h. However, genes for biosynthetic processes of vitamins/porphyrins were downregulated at 2 and 6 h of co-culture, and genes for transmembrane transport and some response networks to external stimuli were also downregulated at 12 h of co-culture (Supplementary Table 4).

Furthermore, the DEGs in Pseudoruegeria sp. were classified using KEGG, with the reference being Ruegeria pomeroyi, which is a member of the marine Roseobacter clade and is phylogenetically close to Pseudoruegeria sp. (Park et al., 2014). Carbohydrate metabolism genes involved in glyoxylate and dicarboxylate metabolism, TCA cycle, and glycolysis were particularly highly upregulated (Figure 5B and Supplementary Figures 2B, 4). Furthermore, energy metabolism, amino acid metabolism, membrane transport, and cellular community-associated genes were together highly upregulated upon co-culture (Figure 5B). Especially, the transcriptomes of the genes associated with oxidative phosphorylation among the category of energy metabolism, such as NADH dehydrogenase (Nuo), succinate dehydrogenase (Sdh), cytochrome C reductase (Pet), cytochrome C oxidase (Cta), and F-type ATPase (Atp), significantly increased after co-culture with A. catenella (Figure 5C).

In addition, the transcription levels of membrane transport genes were observed to be highly upregulated following 2 h of co-culture (Figure 5D). In particular, the genes associated with branched-chain amino acid (BCAA) ABC transporter (livG-F-H-M-K); microcin C transporter (yejA-B-E-F); and oligopeptide transporter (oppA-B-C-D-F), which play a role in transferring dipeptides, branched-chain amino acids, oligopeptides, and proteins into the cytoplasm, were highly upregulated (Figure 5D). In addition, tripartite ATP-independent periplasmic (TRAP) transporter (dctP-Q-M) and tricarboxylate transporter (tctA-B-C) genes used for the uptake of tricarboxylates such as citrate and C4-dicarboxylates as metabolites in the TCA cycle, were also substantially upregulated (Figure 5D).

In both KEGG and GO analyses, genes associated with vitamins and the metabolism of cofactors were downregulated or upregulated according to the type of vitamin. Several genes for porphyrin and chlorophyll metabolism for cobalamin (vitamin B12) biosynthesis and the metabolism of thiamine (vitamin B1), pyridoxine (vitamin B6), and nicotinamide (vitamin B3) were downregulated (Figure 5E and Supplementary Figures 2B, 5).

Interestingly, algicidal effect by co-culture with Pseudoruegeria sp. was also observed in the supernatant after the removal of cells from co-culture, indicating that Pseudoruegeria sp. secretes algicidal compounds into the extracellular environment. To predict the metabolites secreted by Pseudoruegeria sp. that could cause this algicidal effect, we searched the secondary metabolite synthesis clusters, which had been predicted using Antismash (Cho et al., 2020). As a result, several secondary metabolite-synthesizing gene clusters encoding bacteriocin, NRPS/T1PKS, terpene, β-lactone, ectoine, L-homoserine lactone, and lasso peptide were found from the Pseudoruegeria sp. genome sequences (Figure 5E and Supplementary Table 5). Further, when the transcriptome changes of these clusters were examined, the core genes of the peptide toxin bacteriocin, which is known to be harmful to some mammals and other bacteria (Cornut et al., 2008; Cotter et al., 2013), and the lasso peptide gene cluster, which encodes a type of bacteriocin, were highly upregulated under co-culture. In addition, the clusters associated with the synthesis of β-lactone, which is known as a natural algicidal molecule, showed high RNA expression levels and significantly increased transcriptomic expression during co-culture (Figure 5E).

If the algicidal molecule secreted from bacteria is a peptide or protein-based molecule, the bacteria may also use their own secretion system. Accordingly, we monitored the transcription of bacterial secretion systems, including the sec-dependent system, e.g., secA, secDF, secYEG, secB, yajC, ftsY, and ffh, and the twin-arginine translocation system, e.g., tatA and tatBC, for the secretion of peptides/proteins such as bacteriocin. The transcription of such genes also started to increase from 2 h of co-culture and indicated the highest change at 12 h of co-culture (Figure 5F).