Building a locally diploid genome and transcriptome of the diatom Fragilariopsis cylindrus

The genome of the cold-adapted diatom Fragilariopsis cylindrus is characterized by highly diverged haplotypes that intersperse its homozygous genome. Here, we describe how a combination of PacBio DNA and Illumina RNA sequencing can be used to resolve this complex genomic landscape locally into the highly diverged haplotypes, and how to map various environmentally controlled transcripts onto individual haplotypes. We assembled PacBio sequence data with the FALCON assembler and created a haplotype resolved annotation of the assembly using annotations of a Sanger sequenced F. cylindrus genome. RNA-seq datasets from six different growth conditions were used to resolve allele-specifc gene expression in F. cylindrus. This approach enables to study differential expression of alleles in a complex genomic landscape and provides a useful tool to study how diverged haplotypes in diploid organisms are used for adaptation and evolution to highly variable environments.


Background & Summary
How extreme and highly variable environmental conditions impact the evolution and adaptation of eukaryotic organisms has not been extensively studied, yet. Our recent work on the diploid polar diatom Fragilariopsis cylindrus 1 , which has successfully adapted to the sea-ice environment of the Southern Ocean, showed that approximately 25% of its diploid genome is characterized by highly diverged haplotypes that are interspersed throughout its otherwise homozygous genome. Alleles from these diverged haplotypes were differentially expressed across variable environmental conditions. Rather than facilitating an adaptive response through the differential expression of sub-or neo-functionalised genes 2 , the allelic variation within genes of this polar diatom is differentially expressed 1 , which we refer to as environment-dependent Differential Allelic Expression (DAE). Genes with the most pronounced DAE had the largest ratio of non-synonomous to synonomous nucleotide substitutions (d N /d S ), suggesting a correlation between diversifying selection and allelic differentiation, hence a potential role of the divergent alleles for adaptation to environmental fluctuations (e.g., sea-ice formation and melting) in the Southern Ocean.
The most challenging part of this work was to build the locally diploid genome and transcriptome as the genome was randomly interspersed with diverged haplotype contigs assembled with ARACHNE 3 based on Sanger sequencing. To resolve the exact location of these diverged haplotype contigs, we used the diploid aware FALCON assembler 4 for long PacBio sequence reads (Data Citation 1). To resolve genes that have different haplotypes, we lifted the annotation off the ARACHNE assembly, and mapped it to the FALCON assembly, which enabled us to differentiate between alleles and paralogs.
These datasets were mapped onto the FALCON assembly to resolve haplotype-specifc gene expression in the complex landscape of the F. cylindrus genome (Fig. 1). single cell sorted into 96 well plates (InFlux flow cytometer, BD Biosciences, San Jose, CA, USA) and one of these uniclonal and axenic cultures was used for whole-genome sequencing. Genomic DNA for genome sequencing was extracted using a cetyltrimethylammonium bromide (CTAB) method modified from Doyle and Doyle 6 .

DNA library preparation and PacBio sequencing
For genome assembly, we created a 20 kb fragment length library and sequenced 3 SMRT cells with the P6C4 chemistry on a PacBio RSII instrument (Earlham Institute, Norwich, UK). All sequencing reads were collected using standard PacBio single-molecule real-time (SMRT) sequencing protocols (Data Citation 1). The total yield was 1.37 Gb of data. The final N50 of read length varied between 8,215 to 8,898 bp. We also created a 4 kb insert size library and sequenced 4 SMRT cells with the P6C4 chemistry on a PacBio RSII instrument, to generate more accurate Read of Insert (RoI) consensus sequences from multi-pass reads of the same library molecules (Data Citation 1). The total amount of data was 3.85 Gb and the N50 ranged from 2,558 to 2,680 bp between the SMRTcells. We combined the data from these 7 SMRT cells and the total coverage of the raw data was 87x. After filtering the shortest reads, we had 3.8 Gb of data which gave 63x coverage.

Genome assembly with FALCON
We used the diploid aware PacBio assembler, FALCON 0.3.0 (ref. 4) to assemble the F. cylindrus genome. Raw subreads fasta files the 7 SMRTcells were used as input data for FALCON which was run with default parameters, except the minimum cut off for long reads we set to 2,000 bp. See, fc_di_UV.cfg, (github. com/paajanen/FC_scientific_data/) specifying all the parameters used for the Falcon assembly. We note that the length of the genome was not a parameter that the assembler required. The assembly was submitted in an environment containing PYTHON 2.7.10 using the command fc_run.py fc_di_UV.cfg.
The output of the FALCON assembler was divided into two parts. The haploid assembly resulted in primary contigs from which we deduced a genome size of 59.7 Mb. However, the assembler also produced alternate contigs, which represent two diverged haplotypes for those regions (Fig. 2). The assembly metrics of the two sets of contigs are shown in Table 2. We used QUIVER to polish the PacBio assembly using part of the SMRTANALYSIS 2.3.0p5 (pacb. com) pipeline using default settings and the combined PacBio data from 7 SMRT cells in native h5 format of PacBio. In preparation for polishing, we had to rename the scaffolds from the FALCON assembly, by concatenating the fasta names, using a custom script, called rename_falcon_github.py (github. com/paajanen/FC_scientific_data/). We uploaded both the renamed primary contigs and the renamed alternate contigs separately as a reference the polishing using referenceUploader --skipIndexUpdate -c -n "reference" -p /path/to/ working_directory -f /path/to/reference.fasta --saw = "sawriter -blt 8 -welter" --samIdx = "samtools faidx" --jobId = "Anonymous" --verbose. The polishing step was submitted with command smrtpipe.py --params = polishing_params.xml and the parameter file polishing_params.xml is available at github.com/paajanen/FC_scientific_data/.

Annotation of the PacBio genome
We extracted the gene models from the annotation provided by the Joint Genome Institute (JGI) (Fracy1_GeneModels_FilteredModels1_nt.fasta from the JGI website at http://genome.jgi.doe.gov/pages/ dynamicOrganismDownload.jsf?organism = Fracy1) of F. cylindrus and mapped them to the FALCON assembly using GMAP v 20160816 (ref. 8) with the commands gmap_build -D./ -d database assembly.fasta gmap -D /annotation/ -d database -f samse -n 0 -t 16 Fracy1_GeneModels_FilteredModels1_nt.fasta > mapped.sam 2 > mapped. log We used samtools 1.3 (ref. 9) to create a bam alignment file, and used the commands samtools view -bS, samtools sort and samtools index to sort and index it. From the bam file we created a gtf annotation file by using the tool bam2gtf.py from Mikado version 1.0.0 (https://github.com/lucventurini/mikado).

Fosmid sequencing and assembly
The fosmids were Sanger sequenced from paired end sub clones and then computational and experimentally finished using Consed 10 with direct primers walking on sub clones (Data Citation 3).
Fluorescence microscopy combined with 4',6-diamidino-2-phenylindole (DAPI) fluorescent nucleic acid staining was used to confirm axenic cultures before the beginning of culture experiments.
Experimental batch cultures were grown in three biological replicates in chemically defined Aquil artificial seawater media. Fragilariopsis cylindrus cultures were subjected to six different treatments including (1)   F. cylindrus stock cultures from exponential growth phase were used to inoculate 2 l experimental batch cultures with an initial cell count of 50,000 cells ml − 1 . During experimental treatments (except elevated CO 2 treatment), cultures were bubbled with filtered ambient air (Swinnex unit equipped with 25 mm Whatman GF/F filter) passed through milliQ-H2O and manually shaken before subsampling to ensure sufficient CO 2 supply and mixing. Subsamples were taken on a daily basis throughout the experiments to determine physiological parameters including specific growth rate and maximum quantum yield of photosystem II (F v /F m ) as a proxy for cell fitness 14 . Cell counts were determined using automated cell counting with a Multisizer 3 particle counter (Beckman Coulter, Brea, CA, USA) equipped with a 100 μm aperture capillary.
Whilst the experimental treatment of F. cylindrus with elevated carbon dioxide was instantly applied to cell cultures, experimental cultures grown under prolonged darkness, freezing temperatures and elevated temperatures were first grown to early exponential phase at optimal growth conditions before shifting to the final experimental conditions. These experimental treatments were initiated during early exponential phase when cultures had a cell density of approximately 300,000 cells per ml. For low iron treatments, F. cylindrus was grown in iron-free Aquil media that had been passed through a Chelex cation exchange column (Chelex 100 Resin, biotechnology grade sodium form, 100-200 dry mesh size, 150-300 μm wet bead size, Bio-Rad Laboratories, Hercules, CA, USA). Cells from iron-replete stock cultures were transferred into iron-free Aquil media and allowed to grow for several days prior to experimentation to ensure iron limitation as performed previously 15 . Preparation of iron-free Aquil media and handling of low iron cultures were carried out using standard trace metal clean techniques as described for trace metal studies 11,16,17 . Accordingly, 2 l aliquots of Aquil seawater were supplemented with macronutrients (NO 3 , PO 4 and Si(OH) 4 in accordance with Aquil medium concentrations), passed through a Chelex cation exchange column, filter-sterilised (nitrocellulose membrane filter, 47 mm 0.22 μm GSWP, Millipore, MA, USA) and placed into 10% hydrochloric acid-cleaned, milli-Q H 2 O-rinsed 2.5 l polycarbonate bottles. Trace metal concentrations were buffered using 100 μmol l − 1 of ethylenediaminetetraacetic acid (EDTA), which reacts with metal ions (including Fe 3+ ) to metal chelates that are not directly available to phytoplankton, rendering potential iron contaminations insignificant. Dispensed chelexed and filter-sterilised Aquil seawater was supplemented with filter-sterilised (25 mm 0.2 μm syringe filter) EDTA-trace metals (minus iron) and vitamins (B 12 , thiamine and biotin), and allowed to equilibrate chemically overnight at final growth conditions before inoculation of cells. Experimental cultures were sampled for RNA preparations when they reached mid-exponential phase (~500,000 cells ml − 1 ) after 5-10 days of acclimation to the experimental treatment by gentle filtration of cultures (~300 psi vacuum pressure) onto 1.2 μm membrane filters (Isopore membrane, Millipore, MA, USA), placement in 2 ml cryogenic centrifuge tubes and flash-freezing in liquid nitrogen. Finally, the limiting effect of experimental treatments on F. cylindrus was confirmed according to La Roche et al. 18 , which is based on addition of the limiting nutrient to reconstitute optimal growth conditions leading to the recovery of physiological parameters that are depressed by the experimental treatment.

RNA preparation
Total RNA was extracted using a guanidinium thiocyanate-phenol chloroform extraction according to Chomczynski   experimental conditions were run in a single lane of a flowcell using multiplex DNA barcodes, generating paired-end reads of 101 bases length (Data Citation 2).

Code availability
Code is available at github.com/Paajanen/FC_scientific_data including the parameter files for generating the assembly, and the annotation.gtf file.

Data Records
The PacBio genome assembly and sequence data (Data Citation 1) has been deposited at DDBJ/EMBL/ GenBank under accession number PRJEB15040.

RNA quality
Purity of RNA was checked using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and RNA integrity was assessed using 2% denaturating formaldehyde gels and an Agilent 2100 Bioanalyzer electrophoresis system (Agilent, Santa Clara, CA, USA).

Read data quality control
We assessed that the read quality was long enough for assembly by observing the read length distribution of the 20 kb PacBio sequencing library (Fig. 3).

Assembly validation
Assembly size. The haploid assembly resulted in primary contigs from which we deduced a genome size of 59.7 Mb. The genome size was independently estimated as 57.9 (±16.9) Mb using RT-qPCR 1 .
Local accuracy. We assessed the accuracy of this assembly using Sanger finished haplotyped fosmids (Data Citation 3), which we aligned with bwa 0.7.12 (ref. 23) using the command bwa mem -x pacbio.
All fosmids aligned to the assembly. One of the fosmids aligned perfectly over 43,010 bp with no SNPs or  indels, while the worst fosmid aligned for 14 kb region. The lengths of the alignments were calculated using samtools 1.3 (ref. 9) and awk 24 to count positions that were mapped. The counts of indels and SNPs in each of the fosmids were calculated manually and are summarized in Table 3. We calculated the percentage accuracy by adding the number of SNPs and indels and divided it by the length of alignment.
Mapping the reads back to the assembly. We created circular consensus (ccs) reads from the 4 kb insert size library using the SMRTANALYSIS 2.3.0p5 pipeline and a cut-off of 500 bp as minimum read length with the command smrtpipe.py --params = ReadsOfInsert.params.xml xml:input.xml > smrtpipe.log where input.xml was created from a file input.fofn containing paths to the h5-files from the SMRTcells from the 4 kb libraries fofnToSmrtpipeInput.py input.fofn > input.xml and ReadsOfInsert.params.xml (github.com/paajanen/FC_scientific_data). The circular consensus reads had an average number of 10 passes, and the mean read quality of the inserts was >0.9862 for all four SMRT cells. The total number of reads was 119,511 with a mean of 2,429 bp. The total number of bases of ccs reads was 290,305,468 bp, resulting in a coverage of 4.8x over a genome of size ca. 60 Mb. We aligned these reads to the primary and alternate contigs of the assembly using bwa 0.7.12 (ref. 23

Quality of annotation
We searched the gtf annotation file with the custom python script coverage_identity.py (github. com/paajanen/FC_scientific_data) and verified that 24,635 out of the 26,418 transcripts (93.3%) from the JGI annotation map to the FALCON assembly with at least 95% coverage and 95% identity.

Allelic pairs in the FALCON assembly
Using the 24,535 transcripts from the JGI annotation that map to the FALCON assembly, we updated the list of allelic pairs identified in the ARACHNE assembly 1 , using the custom script SD. find_allelic_pairs.py (github.com/paajanen/FC_scientific_data). Out of the 6,071 allelic pairs identified in the ARACHNE assembly 5,400 (89%) are found in the FALCON assembly. A total of 1,894 of the 5,400 allelic pairs in the FALCON assembly map in to alternate haplotypes, while the rest are ambiguous, and are either paralogs or located on smaller scaffolds that are not easily distinguished as alternate haplotypes. Additionally, we identified 1,615 transcripts on alternate haplotype contigs that are not assigned as allelic pairs in the ARACHNE assembly. Some of these have identical copies at the corresponding location of the primary contig, possibly because polishing the assembly has polished out the differences. Lists of JGI protein identifiers for (1) confirmed alternate allelic pairs in the ARACHNE assembly, (2) transcripts that have alternate alleles in both, ARACHNE and FALCON assembly, and (3) new possible transcripts with alternate allele are in Supplementary Table 1.

RNA-seq data validation
Metrics of the RNA-seq dataset. A total of 421 million reads were generated, with an average of 3.8 million read pairs per sample (elevated temperature, elevated carbon dioxide, darkness, optimal growth, low iron, low tempertature). The reads were analysed with FasQC software, which did not flag any sample as bad quality. The average read quality was 36 (mean) and 37 (median) for each sample. The samples had negligible amounts of adapter residues. The read length was 101 bp.
Alignment validation. The percentage of reads aligning uniquely varies from 68 to 74%, depending on the experiment (Table 1). We also counted the number of mapped reads for each experiment with htseq version 0.6.1 (ref. 25) using the command htseq-count --mode = intersection-nonempty --stranded = no --type = exon Aligned.sam annotation.gtf.

Usage Notes
The genome assembly (ref. 1) sequencing reads described in this work were intended to confirm the published assembly 1 and to clarify the diverged haplotypes. We have not annotated the PacBio FALCON assembly from scratch, but rather used the annotation of the ARACHNE assembly produced by JGI to map to the FALCON assembly. Therefore, the ARACHNE assembly based on Sanger sequencing might be more suitable for detailed genetic studies, while the FALCON assembly should be used for longer range information and precise location of allelic diversity. RNA-seq data can be used to map to either assembly based on individual needs.