Evolutionary genomics of the cold-adapted diatom Fragilariopsis cylindrus

The Southern Ocean houses a diverse and productive community of organisms. Unicellular eukaryotic diatoms are the main primary producers in this environment, where photosynthesis is limited by low concentrations of dissolved iron and large seasonal fluctuations in light, temperature and the extent of sea ice. How diatoms have adapted to this extreme environment is largely unknown. Here we present insights into the genome evolution of a cold-adapted diatom from the Southern Ocean, Fragilariopsis cylindrus, based on a comparison with temperate diatoms. We find that approximately 24.7 per cent of the diploid F. cylindrus genome consists of genetic loci with alleles that are highly divergent (15.1 megabases of the total genome size of 61.1 megabases). These divergent alleles were differentially expressed across environmental conditions, including darkness, low iron, freezing, elevated temperature and increased CO2. Alleles with the largest ratio of non-synonymous to synonymous nucleotide substitutions also show the most pronounced condition-dependent expression, suggesting a correlation between diversifying selection and allelic differentiation. Divergent alleles may be involved in adaptation to environmental fluctuations in the Southern Ocean.

The Southern Ocean houses a diverse and productive community of organisms 1,2 . Unicellular eukaryotic diatoms are the main primary producers in this environment, where photosynthesis is limited by low concentrations of dissolved iron and large seasonal fluctuations in light, temperature and the extent of sea ice [3][4][5][6][7] . How diatoms have adapted to this extreme environment is largely unknown. Here we present insights into the genome evolution of a cold-adapted diatom from the Southern Ocean, Fragilariopsis cylindrus 8,9 , based on a comparison with temperate diatoms. We find that approximately 24.7 per cent of the diploid F. cylindrus genome consists of genetic loci with alleles that are highly divergent (15.1 megabases of the total genome size of 61.1 megabases). These divergent alleles were differentially expressed across environmental conditions, including darkness, low iron, freezing, elevated temperature and increased CO 2 . Alleles with the largest ratio of non-synonymous to synonymous nucleotide substitutions also show the most pronounced condition-dependent expression, suggesting a correlation between diversifying selection and allelic differentiation. Divergent alleles may be involved in adaptation to environmental fluctuations in the Southern Ocean.
The pennate diatom genus Fragilariopsis is especially successful in the Southern Ocean, with the cold-adapted species F. cylindrus ( Fig. 1a) regarded as an indicator species for polar water [8][9][10] . It is frequently found to form large populations in both the bottom layer of sea ice and the wider sea-ice zone, including open waters 9 (Fig. 1b). Sea ice is characterized by temperatures under 0 °C, high salinity and, owing to the semi-enclosed pore system within the ice, low diffusion rates of dissolved gases and exchange of inorganic nutrients 11 . However, unlike in ice-free surface waters of the Southern Ocean 12 , dissolved iron is not considered to be limiting to phytoplankton growth within sea ice 13 . Most phytoplankton in the Southern Ocean face inclusion into sea ice every winter and are released again in summer when most of the sea ice melts 14 ; certain species such as F. cylindrus have therefore evolved adaptations to cope with this drastic environmental change. Thus, comparative analyses of the genome of the psychrophile F. cylindrus with those of diatoms that evolved in temperate oceans provide an opportunity to obtain insights into how this species has adapted to conditions in Southern Ocean surface waters.
We found many loci with highly divergent alleles in the diploid F. cylindrus draft genome sequence. To resolve the divergent alleles from paralogous genes, we independently carried out Sanger and PacBio sequencing and used haplotyped Sanger-finished fosmids to validate the haplotype-resolved genome assemblies (Supplementary Data 1-3). Using complementary approaches, we found that the F. cylindrus genome assembly consists of 15.1 Mb of loci with highly divergent alleles that were assigned to different scaffolds. The remaining 46 Mb of sequence consists of alleles similar enough to be assembled onto the same scaffold (Supplementary Information 2-5). The haplotype assembly size of the genome (61.1 Mb; Extended Data Table 1) was confirmed by quantitative PCR with reverse-transcription (qRT-PCR) (57.9 Mb). The genome completeness according to the Core Eukaryotic Genes Mapping Approach 15 is 95.6% and the nuclear scaffold N50/L50 is 16/1.3 Mb, corresponding to assembly size (Extended Data Table 1).
The haplotype-resolved genome contains 21,066 predicted protein-coding genes (Extended Data Table 1) with 6,071 genes (29%) being represented by diverged alleles (Allele sets 1 and 2, Supplementary Data 1). Sequence divergence between alleles was up to 6%, but this was still significantly less (Mann-Whitney, P < 0.001) than that between paralogous genes (Extended Data Fig. 1 and Supplementary Information 4,5).
We compared the F. cylindrus genome with those of Thalassiosira pseudonana 16 and Phaeodactylum tricornutum 17 (Extended Data  Table 1), both of which live in temperate and neritic marine environments 18,19 characterized by higher water temperatures, turbidity and concentrations of dissolved iron. The haploid gene content of F. cylindrus is enriched for two conserved metal-binding protein families (structural classification of proteins (SCOP) fold families; Supplementary Information 6). When accounting for its genome size, it is enriched for copper-binding but not iron-binding proteins ( Fig. 1c), and it contains a disproportionate abundance of domains belonging to the plastocyanin/azurin-like family fold (SCOP ID 49504). Copper-containing plastocyanin may facilitate photosynthetic electron transport, reducing the need for iron 20 . There also appear to be more zinc-binding proteins in the F. cylindrus genome than in the other genomes, with 121 proteins containing zinc-binding myeloid-Nervy-DEAF-1 (MYND) domains, compared to 7 in T. pseudonana and 12 in P. tricornutum (Fig. 1d). MYND domains facilitate protein-protein interactions and are involved in regulatory processes 21 ; most of those in F. cylindrus appear to be lineage-specific ( Supplementary Information 6). Evolutionary genetic analysis of MYND-containing proteins suggests that this family has expanded within the last 30 million years ( Supplementary Information 7). The relatively high zinc concentration of Southern Ocean surface waters 22 may have facilitated the great expansion and functional divergence of zinc-binding MYND domains. The presence of lineage-specific protein families might indicate specific adaptations to the extreme conditions in the Southern Ocean. Some of these protein families appear to have been acquired through horizontal gene transfer from bacteria ( Supplementary Information 8). Those proteins include groups of ice-binding proteins 23 and protonpumping proteorhodopsins 24 (Supplementary Information 9, 10). There is also an unusually large number of genes for chlorophyll a/c light-harvesting complex (LHC) proteins, including 11 members of the Lhcx clade, which is involved in stress response (Extended Data Fig. 2 and Supplementary information 11).
Gene Ontology analysis show significant enrichment of genes in the categories 'catalytic activity' , 'transporter activity' , 'metabolic process' , 'transport' and 'integral to membrane' in the group of diverged alleles compared to the non-diverged alleles (Fisher's exact test, adjusted P < 0.05; Extended Data Table 2 and Supplementary Information 12). We found that similar processes (for example, transport and metabolic process) were enriched in metatranscriptome sequences from Southern Ocean sea ice (Figs 1b, 2a and Supplementary Data 4) with strong homology (BLASTx, E value ≤ 1 × 10 −10 ) to F. cylindrus protein sequences of diverged alleles (Supplementary Information 13). According to these cut-off criteria, 64% of all Bacillariophyta-like metatranscriptome sequences had homology with proteins in F. cylindrus and around 60% of these sequences matched diverged alleles in the genome of F. cylindrus, including sequences from the enriched Gene Ontology categories (Supplementary Information 13).
RNA sequencing (RNA-seq) transcriptome profiling under environmentally relevant growth conditions (darkness, low iron, freezing, elevated temperature and CO 2 ) identified stress-specific responses (Fig. 2b). The broadest transcriptome response (approximately 60% of total genes, including divergent alleles) was observed under prolonged darkness, characteristic of polar winters (Supplementary Information 14 and Supplementary Data 5). Placing F. cylindrus in darkness for seven days downregulated genes involved in photosynthesis, light harvesting and photoprotection relative to their expression under continuous light. By contrast, genes involved in starch, sucrose and lipid metabolism were strongly upregulated in the dark (Extended Data Fig. 3), indicating the utilization of chrysolaminarin and fatty acid storage products. Notably, under prolonged darkness, the percentage of RNA-seq reads that did not map to predicted genes (30%) was higher than under any other tested growth condition.
In allele-specific analyses of transcriptomes, approximately 66% (4,030) of diverged alleles showed greater than fourfold significant differential expression (likelihood ratio test, P < 0.001) relative to optimal nutrient-replete growth (Fig. 2b) and approximately 45% (2,730) of divergent alleles showed greater than fourfold unequal biallelic expression between allele 1 and allele 2 in at least one RNA-seq experiment (likelihood ratio test, P < 0.001; Supplementary Data 6). Additionally, the functional significance of this unequal bi-allelic expression for metabolism was inferred by an individual analysis of both sets of alleles using Gene Ontology. This demonstrated different metabolic signatures between the groups of divergent alleles (Fig. 2c).

Letter reSeArCH
The differential expression of divergent alleles in response to environmental stresses suggests that individual alleles may be under different regulatory controls. Generally, variations in allelic expression have been attributed to differences in non-coding DNA sequences and epigenetic regulation 25 . Notably, nucleotide sequence analysis of gene promoter and coding regions of all diverged allelic pairs revealed a significantly lower (P = 1.0 −23 ) sequence identity in promoter regions (Extended Data Fig. 4 and Supplementary Information 15), which suggests functional diversity in allelic promoter regions 26 .
To test whether the divergent alleles may be the consequence of adaptive evolution to distinct environmental conditions, we divided the allelic pairs into seven subsets according to their ratio of nonsynonymous to synonymous nucleotide substitutions (d N /d S ) (Fig. 3a, b). Alleles with an elevated rate of non-synonymous mutations showed a significantly higher maximum fold change in bi-allelic expression during RNA-seq experiments ( Fig. 3a; median test, adjusted P < 0.05). The highest median log 2 fold-change in bi-allelic expression was 2.73-this was observed for the subset of diverged alleles with d N /d S ≥ 1, which is indicative of positive selection. The lowest median log 2 fold-change was 2.01-this was observed for the subset of alleles with d N /d S 0-0.1, the smallest range (Fig. 3a). This suggests that positive selection has a role in driving the evolution of alleles with strong bi-allelic expression. However, most of the alleles with the highest d N /d S had unknown functions (Extended Data Table 3).
If allelic divergence is important for adaptation to a fluctuating environment, one might predict that recombination would be suppressed. We therefore examined the effect of recombination and genetic drift on the allelic variation, studying a natural population of    (1,000 p.p.m. CO 2 ) and prolonged darkness, relative to optimal growth conditions. Each experimental treatment corresponds to two separate columns for both allelic variants and each single-haplotype gene to a single row. c, Allele set-specific REViGO semantic similarity scatter plots for biological process Gene Ontology terms, showing relative differences between allele 1 and allele 2 in F. cylindrus RNA-seq experiments under low temperatures and low iron relative to optimal growth conditions using the cut-off ≥ 0.3 for relative difference.
F. cylindrus from Southern Ocean sea ice (Fig. 1b). We analysed around 200 high-quality Sanger sequences from alleles of two genes, the ferrichrome ABC transporter (Joint Genome Institute (JGI) protein ID 240308) and large ribosomal protein L10 (JGI protein ID 267462). Recombination analysis identified various intragenic recombinant alleles consistent with reticulate evolution (Fig. 4a, b and Supplementary Information 16). We then analysed the phylo genetic networks of these alleles and compared the branch lengths and the number of splits to networks of simulated populations. In addition, we compared the alleles of 645 genes to homologous alleles from mate-pairs of the temperate diatom Pseudo-nitzschia multistriata, a closely related sexually reproducing species, showing that the alleles in F. cylindrus have an overall higher allelic diversity than those of P. multistriata ( Fig. 4c and Supplementary Information 16). These analyses indicated that the extensive allelic diversity in F. cylindrus is maintained in a vast gene pool with an effective population size N e ≈ 16.5 × 10 7 (assuming a base mutation rate μ = 10 −9 ), and that the recombination rate is about five times the mutation rate (Fig. 4d, e). The observed divergence is thus not the result of genetic introgression after hybridization, but simply the consequence of a high mutation-drift parameter (Θ) in conjunction with positive selection. Furthermore, alleles in the genome of F. cylindrus appear to coalesce shortly after the onset of the last glacial period, (1) (90)   networks (a, b), the recombination rate is between 1 and 10 times the mutation rate (R/μ ≈ 5).

Letter reSeArCH
which began about 110,000 years ago 27 (Extended Data Fig. 5). Thus, our studies suggest that the diversification of alleles took place only recently and is maintained in the vast gene pool of the diatom, which allows it to thrive under the highly variable environmental conditions of the Southern Ocean 2,3,5,7,13 .
Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper.

MEthOdS
No statistical methods were used to predetermine sample size. The experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment. Culture strain and DNA preparation. F. cylindrus (Grunow) Krieger (strain accession CCMP1102) was obtained from the National Centre for Marine Algae and Microbiota. Bacterial contaminants were removed by treatment with ampicillin (50 μ g ml −1 ) and chloramphenicol (1 μ g ml −1 ) and cultures were single cell sorted using flow cytometry. High-molecular-weight DNA for whole-genome sequencing was extracted from an axenic and monoclonal culture as previously described 28 with minor modifications (see Supplementary Information 1); genome size was estimated using quantitative real-time PCR as previously described 29 . Sanger sequencing. All sequencing reads were collected using standard Sanger sequencing protocols on ABI 3730XL capillary sequencing machines at the US Department of Energy Joint Genome Institute. Three different sized libraries were used as templates for the plasmid subclone sequencing process and both ends were sequenced. We obtained 270,371 reads from the 2.5-kb library, 319,392 reads from the 6.3-kb library, and 81,408 reads from a 35.4-kb fosmid library. PacBio sequencing. All sequencing reads were collected using standard PacBio single-molecule real-time (SMRT) sequencing protocols. Two different-sized libraries were created and sequenced on a PacBio RSII instrument using the sixth generation of polymerase and the fourth generation of chemistry (P6-C4 chemistry). A 20-kb fragment length library was sequenced using three SMRT cells with total yield of 1.37 Gb of raw data. Additionally, a 4-kb insert size library was sequenced using four SMRT cells with a yield of 3.85 Gb of raw data. Sanger assembly. A total of 671,171 sequence reads (7.25× final sequence coverage) were assembled using Arachne v.20071016 (ref. 30). The final genome assembly was produced by passing the initial assembly through the Rebuilder module to merge adjacent haplotypes, followed by another complete Arachne assembly process. We obtained 4,622 contigs that were linked into 286 scaffolds, including 105 scaffolds larger than 100 kb. To exclude organelle sequences and contaminating scaffolds from the nuclear genome assembly, each scaffold was screened against bacterial proteins and organelle sequences in the NCBI GenBank database and a set of known microbial proteins using megaBLAST and blastp searches 31 . Additional scaffolds were removed if they consisted of more than 95% 24-mers that occurred four other times in the scaffolds larger than 50 kb or if the scaffold contained only unanchored RNA sequences. We classified additional scaffolds as one chloroplast scaffold, five mitochondrial scaffolds, two scaffolds of < 1 kb length, and seven small repetitive scaffolds. The final nuclear genome assembly contains 4,602 contigs with a contig N50 (that is, the contig size above which 50% of the total length of the sequence is included) of 78.2 kb and 271 scaffolds with a scaffold N50 of 1.3 Mb. The genome completeness was assessed using the Core Eukaryotic Genes Mapping Approach (CEGMA) 15 . The cumulative haploid genome size was estimated at 61.2 Mb, accounting for 46 Mb genomic scaffolds that were collapsed into a single haplotype, 29.8 Mb of genomic scaffolds that could not be collapsed into a single haplotype (that is, 14.9 Mb for collapsed single haplotype; see below), and 0.3 Mb low-coverage scaffolds. This is consistent with independently estimated genome sizes of 57.9 (± 16.9) Mb and 59.7 Mb using qPCR and PacBio sequence data, respectively. Additionally, we experimentally validated the haplotype-resolved genome assembly from whole-genome shotgun Sanger sequences by sequencing a large-insert fosmid library and aligning it to genomic scaffold sequences, using the contiguity information of fosmids to directly phase ascertained collapsed (homologous) and diverged haplotypes. We then assessed nucleotide alignments between annotated protein sequences from the genome assembly scaffolds and the haplotyped Sanger-finished fosmid clones that were not included in the genome assembly. Finally, under the assumption that gene duplicates are more divergent than alleles, we compared the sequence similarity between predicted diverged alleles and duplicates. We independently validated the sequence similarity for the alleles by Sanger sequencing of diverged F. cylindrus alleles from a natural sea ice population (see 'Evolutionary genetic analysis' below). PacBio assembly. PacBio data from seven SMRT cells were used, and after filtering the shortest reads, we obtained 1,971,632 reads and 3.8 Gb of data which gave 63× coverage. We used the diploid aware PacBio assembler FALCON 0.3.0 (ref. 32), which has recently been used to successfully assemble highly heterozygous genomes 33 , to assemble the F. cylindrus genome. The cut off for long reads was 2,000 bp. The FALCON assembly consisted of 745 primary contigs with a total length of 59.7 Mb. The N50 of the primary contigs was 245 kb. The assembler also produced alternate contigs, which represent two diverged haplotypes for those regions. There were 288 alternate contigs, with N50 of 42 kb and total length 9.1 Mb. We used the QUIVER algorithm to polish the PacBio assembly using parts of the smrtanalysis 2.3.0p5 (http://www.pacb.com) pipeline. We assessed the accuracy of this assembly using the Sanger finished haplotyped fosmids, which we aligned with bwa 0.7.12 (ref. 34) using the bwa mem -x pacbio command. The polished assembly was highly accurate, ranging from 99.65% to 100%, with all fosmids aligning to it. One of the fosmids aligned perfectly over 43,010 bp. Out of the remaining 13 fosmids, 8 had an accuracy of > 99.9%. Genome annotation. The F. cylindrus genome assembly was annotated using the Joint Genome Institute (JGI) annotation pipeline 17 . The assembly was masked for repeats using RepeatMasker 35 and the RepBase library 36 , and the most frequent (> 150 times) repeats were recognized by RepeatScout 37 . Protein-coding gene models were predicted using three levels of evidence: ab initio Fgenesh 38 ; homology-based Fgenesh+ 38 and Genewise 39 seeded by BLASTx alignments against the NCBI NR database; and transcriptome-based Fgenesh. For each genomic locus, automated filtering selected the best model based on homology and transcriptome support. tRNAs were predicted using tRNAscan-SE 40 . All predicted proteins were functionally annotated using SignalP 41 60 . The assembly and annotation were released as a public web portal available at http://genome.jgi-psf.org/Fracy1/Fracy1.home.html. Identification and analysis of diverged alleles. To produce a non-redundant single haplotype gene set with only one allele of each gene, we aligned the genomic assembly against itself using BLAT 61 with thresholds of 95% nucleotide identity and ≥ 50% alignment coverage for the smaller scaffold. We obtained alignments of 210 smaller scaffolds against larger scaffolds, with a total length of 15.9 Mb and an alignment coverage over the entire length of the smaller scaffold for 74.3% of the alignments. Syntenic scaffolds that were homologous over the entire length were analysed using Mauve genome alignment software 62 . Gene models of the aligned smaller scaffolds, which formed the best bi-directional blastn hit pairs with corresponding genes on larger scaffolds and > 90% nucleotide identity were removed to produce a final non-redundant protein-coding gene model set. We also defined the set of diverged alleles for allele-specific downstream analyses. Diverged alleles on the larger scaffolds were referred to as allele 1 set and alleles on the smaller scaffolds as allele 2 set.
Additionally, the assembly based on PacBio sequencing was used to validate the allelic divergence observed in the assembly based on Sanger sequencing. For this analysis, the 15 longest scaffolds from the PacBio assembly were used, which accounted for 21 Mb of primary sequence and 2 Mb of alternate sequence. After filtering for genes that were only deviant by ± 1% from the length of the predicted protein-coding gene models, we identified 305 genes that possessed two diverged alleles and 30 genes that were classified as paralogues. Functional Gene Ontology enrichment analysis of diverged alleles. An intragenomic Gene Ontology enrichment analysis was performed to test whether diverged alleles were enriched for functional Gene Ontology classes in comparison to non-diverged alleles. We compared the proportion of diverged allelic pairs associated with a Gene Ontology class in the total set of Gene Ontology annotated diverged allelic pairs against the same proportion calculated for the set of nondiverged alleles using Fisher's exact test and adjusted P values 63 . Promoter analysis of diverged alleles. Putative promoter nucleotide sequences were collected by extracting good quality sequences from different regions relative to the transcription start site (TSS). The collected sequences were divided into promoters (before TSS) and transcripts (after TSS) and ClustalW 64 alignments of both sets were parsed with custom scripts using Bioperl 65 . We calculated the average identity of the alignments and the average percentage identity in 10-bp intervals using a sliding window approach.

Letter reSeArCH
Environmental metatranscriptome signature of F. cylindrus. Sequences from a Southern Ocean, eukaryote-targeted metatranscriptome were quality filtered, clustered and taxonomically classified using PhymmBL 66 with a custom reference database 6 . To identify F. cylindrus-like transcripts a BLASTx search (E value ≤ 1 × 10 −10 ) of sequences classified as Bacillariophyta (PhymmBL confidence score ≥ 0.9) was performed against all predicted gene models in the genome assembly, including gene models for diverged alleles on scaffolds that could not be collapsed into a single haplotype. The total number of hits was then calculated for all genes, including genes present as diverged alleles. We used a functional Gene Ontology analysis to assess the abundance of metatranscriptome reads associated to diverged alleles and investigate the environmental significance of bi-allelic transcript abundance for diverged alleles. Results were visualized with REViGO 67 . Transcriptome sequencing. F. cylindrus batch cultures were grown in three biological replicates under optimal growth conditions (+ 4 °C, nutrient-replete Aquil 68 , 24 h light at 35 μ mol photons m −2 s −1 ), freezing temperatures (− 2 °C), elevated temperatures (+ 11 °C), elevated carbon dioxide (+ 4 °C, 1,000 p.p.m. CO 2 ), low iron (+ 4 °C, iron-depleted Aquil), and prolonged darkness (+ 4 °C, 7 days darkness). Total RNA was extracted using an adaptation of the acid guanidinium thiocyanate-phenol-chloroform method 69 . cDNA libraries for RNA-sequencing were constructed from total RNA using random hexamer primers and sequenced in a single-flow cell lane using multiplex DNA barcodes on an Illumina HiSeq 2000 instrument at Edinburgh Genomics to generate paired-end reads of 101 bases in length. A total of 68,832,506 RNA-seq reads were aligned to the F. cylindrus genome assembly using GSNAP 70 , and HTSeq 71 was used to count unique fragments mapping in each genomic feature. Differential expression analysis. Data were analysed using edgeR 72 for differential expression and goseq 73 for functional Gene Ontology analysis. Results were visualized using REViGO 67 and graphical packages of R statistical software 74 . Statistical testing for genes that were differentially expressed between an experimental treatment and optimal growth reference condition was performed using the generalized linear model (glm) functionality implemented in edgeR. After estimating genewise (tagwise) dispersions and fitting negative binomial models, the glm likelihood ratio test 75 was applied to test for differentially expressed genes and P values adjusted 63 . Testing for differentially expressed diverged alleles was performed analogously. Differential bi-allelic expression was analysed comparing the expression of diverged alleles between experimental treatment and optimal growth reference conditions, and within diverged allelic pairs for each single growth condition. Bi-allelic expression relative to allelic divergence. To test the role of adaptive evolution in allelic divergence we investigated the relationship between the d N /d S and the degree of bi-allelic expression under all experimental conditions (maximum differential expression within diverged allelic pairs). To calculate the d N /d S ratios for diverged allelic pairs their nucleotide transcript sequences were translated into amino acids and aligned with ClustalW2 (ref. 64). Amino acid alignments were mapped back over nucleotide sequences to ensure that nucleotide sequences contained full codons and were in frame. After realignment of adjusted nucleotide sequences the d N /d S was calculated for each allelic pair using codeml (pairwise mode) within PAML 76 . Outlier genes showing abnormally high d N /d S > 10 were discarded for subsequent analysis. The diverged allelic pairs were divided into subsets according to their associated d N /d S ratios and differences in the maximum differential bi-allelic expression between groups were compared using nonparametric statistical testing. Sequencing of diverged F. cylindrus alleles from a natural sea ice population. A total of 196 sequences for alleles of two genes encoding a ferrichrome ABC transporter (JGI protein ID 240308) and the large ribosomal protein L10 (JGI protein ID 267462) were amplified from environmental samples using species-specific PCR primers. Purified PCR fragments were cloned using a TOPO cloning strategy, sequenced using capillary sequencing technology (Sanger method), and manually inspected for sequence quality. The sequence divergence between the two allelic pairs was assessed and compared to the sequence divergence of all alleles and duplicates as predicted for haplotype-resolved genome assemblies based on Sanger and PacBio sequencing. Recombination analysis. We developed a novel approach to identify triplets in the diverged allelic sequences that show evidence of homologous recombination ( Supplementary Information 16). Additionally, we used the R package HybridCheck 77 and pairwise homoplasy index testing 78 to identify sites of recombination. Phylogenetic network analysis. The evolutionary relationships between the diverged alleles encoding for the ferrichrome ABC transporter and the large ribosomal protein L10 was visualized in splits graphs using SplitsTree4 (ref. 79), and branch lengths and number of splits in the observed phylogenetic networks were compared with those from simulated populations using simuPOP 80 . Population genetic simulations assuming a base mutation rate μ = 10 −9 were performed across a range of values for the population mutation parameter theta (θ = 4N e μ) and the recombination rate (R) to estimate effective population size (N e ) and recombination rates relative to the mutation rate (R/μ). To minimize the effects of selection, the substitution patterns at the third codon position were studied only. DNAsp 56 and LAMARC 81 were used to estimate theta from sequences of the diverged alleles. Comparative analysis of allelic nucleotide divergence between a polar and temperate diatom. The nucleotide divergence of alleles for 645 F. cylindrus genes was compared to homologous alleles from RNA-seq transcriptomes of mate pairs of the sexually reproducing temperate diatom P. multistriata. Alleles in P. multistriata were identified using best reciprocal blastn (≥ 90% overall identity, ≥ 75% coverage of both sequences) searches of the two mate pair strains. Homologous alleles between both diatoms were identified using reciprocal tBLASTx (≥ 30% overall identity, ≥ 50% coverage of the query sequence) searches of the theoretical six frame translations of the sequences. The divergence of allelic pairs was calculated as described above using PAML 76 . Coalescence time estimates of diverged F. cylindrus alleles. The coalescence time between the two independently sequenced diverged F. cylindrus alleles from a natural population (ferrichrome ABC transporter and large ribosomal protein L10) were estimated and compared to the coalescence time of diverged allelic pairs of the approximately 6,000 genes in the genome assembly. Coalescence time was estimated using the algorithm implemented in HybridCheck 77 . The coalescence time estimate returned by the algorithm is expressed in terms of generations and was converted into years using an estimated division rate of 12.5 per year for F. cylindrus. Data Availability. The Sanger genome assembly and annotations used in this study are available via the JGI genome portal at http://genome.jgi.doe.gov/ Fragilariopsis_cylindrus and the Whole Genome Shotgun Project has been deposited to DDBJ/EMBL/GenBank under accession number LFJG00000000 (BioProject PRJNA32761). ESTs have been deposited to DDBJ/EMBL/GenBank under accession number GW070125. The PacBio genome assembly has been deposited at DDBJ/EMBL/GenBank under accession number PRJEB15040. F. cylindrus RNA-seq data are available in the ArrayExpress database (http:// www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-5024. Metatranscriptome sequences have been deposited at the Sequence Read Archive under accession number SRR1752079. P. multistriata RNA-sequencing data has been deposited at DDBJ/EMBL/GenBank under accession number PRJNA80045 and is available at http://www.ebi.ac.uk/ena/data/view/SRS190381 (strain Sy373) and http://www.ebi.ac.uk/ena/data/view/SRS190382 (strain Sy379). All other data are available from the author upon reasonable request.