Dataset: The Long Term Phytoplankton Evolution Experiment: Genomic Analysis

FastQ files (NCBI BioSamples SAMN34542194 through SAMN34542251) were analyzed using breseq with default settings. Synechocystis PCC6803 cultures were assembled against the chromosomal reference genome (NCBI accession number NC_000911) as well as its four plasmids (NC_005229, NC_005230, NC_005231, and NC_005232). Prochlorococcus MIT9312 and Synechococcus CC9311 were assembled against their RefSeq genomes (respectively, CP000111 and NC_008319). The T. oceanica CCMP1005 and E. huxleyi CCMP371 reference genomes were still in draft form and were assembled in breseq with the contig flag activated; all contigs were downloaded as a single GenBank format file from NCBI BioProject PRJNA36595 (CCMP1005) or BioSample SRX112492 (CCMP371, also known as strain 12-1). The T. oceanica CCMP1005 sequence file contained the T. oceanica chloroplast sequence but E. huxleyi CCMP371 did not; we therefore used the chloroplast sequence of another E. huxleyi strain, CCMP373, also with downstream curation, to analyze mutational changes to those sequences. For both T. oceanica and E. huxleyi, breseq was set to not attempt to predict structural variants due to the much greater computational demand required to assemble these large genomes relative to the bacteria. Alteromonas genomes were assembled against the most recent EZ55 genome (40), including both the chromosomal (CABDXN010000001) and plasmid (CABDXN010000002) sequences in the same GenBank format file. For Prochlorococcus and Synechococcus lineages, EZ55 genomes were assembled in the same breseq run as the phytoplankton genome. However, because of the computational demand of assembling T. oceanica and E. huxleyi genomes, the Alteromonas portion of these cultures was assembled separately to allow for structural variant prediction. All reference genomes are included with this archive in the format used for the analysis.

Any specific mutational call present in 100% of replicate lineages for a given organism was considered to have occurred prior to the initiation of the experiment (i.e., it was ancestral) and was removed from further consideration. Because of additional genomic complexity in the eukaryotic genomes of T. oceanica and E. huxleyi, we undertook additional efforts to curate these datasets. First, coverage for the eukaryote genomes was lower than for the bacterial genomes, and some loci had insufficient coverage to make predictions about mutations; all these loci were removed from breseq output files prior to further analysis. Also, because the GenBank files we used as reference sequences represented haploid sequences, we sought to discover potentially heterozygous loci in the ancestral genomes. To do this, we re-assembled the raw read files from the original genomic sequencing runs for both organisms against the GenBank reference sequences using breseq (SRX112492 for E. huxleyi and all sequencing runs within BioProject PRJNA36595 for T. oceanica). We identified several variants present in 100% of reads compared to the published GenBank sequences, perhaps suggesting differences between breseq’s handling of data and the programs used to produce the published sequences; all mutations assigned by breseq to these loci were therefore removed from consideration as untrustworthy. All variants present in 0% < n < 100% of sequences in the reference genome reassemblies were assumed to represent loci that were heterozygous in the ancestral population. Mutations mapped to these loci were only considered further if they became fixed in evolved lineages; if they remained at intermediate frequencies they were removed from subsequent analyses.

Mutational analysis: We used the application gdtools from the breseq package to convert breseq output genome difference files into long-format data files for subsequent analysis; these are provided in .tsv format with this archive. Custom python scripts were used to bin all mutations within a given gene in each lineage, producing mutational count tables of genes vs. lineages, also provided with this archive. For bacterial genomes, we produced count tables either with or without considering insertions, deletions, and intergenic mutations within 50 bp upstream of a gene’s start codon (i.e., putative cis-acting promoter or other regulatory elements) in addition to mutations within a gene’s coding region. Eukaryotic count tables only considered coding sequences; all intergenic and intron regions were removed from analysis. We used python scripts to analyze three metrics of directional selection: dN/dS (56), transition:transversion ratio (57), and nonsense mutation ratio (58). All three metrics exclusively used SNP data from coding regions of annotated genes. Effects of treatment groups on these metrics were analyzed using linear models in R.

We sought to test for convergent evolution of gene targets by determining which, if any, genes received more mutations than would be expected by chance under a completely random mutational model. This analysis was complicated by the fact that larger genes represented larger mutational targets than smaller ones, precluding the use of a simple Poisson model. Instead, we used a bootstrapped Monte Carlo procedure to produce randomized genomes with the same number of coding sequence mutations observed in our real dataset, only distributed randomly across the genome. Given the total coding genome size g as the sum of the lengths of all coding sequences in the genome, and the total number of observed mutations n, the probability p of any given base pair receiving a mutation is p = n/g, and the probability l of a gene of length l receiving a mutation is thus l = p x l. For each lineage, we produced 100 randomized matrices where each gene in the genome was assigned several mutations drawn from a Poisson distribution with average rate of occurrence l. Each random matrix was compared to the real matrix of mutations by first converting each into an empirical cumulative distribution function (ECDF) and then applying either a Kolmogorov-Smirnov or dts test. The Kolmogorov-Smirnov (KS) test only considers the single value in the ECDF where the gap between the two samples is greatest, whereas the dts test considers the entire distribution. Because the greatest difference between our observed and Monte Carlo matrices was the fact that the real dataset had many more genes with large numbers of mutations than were ever observed in any simulation, the dts test was generally more sensitive in detecting significant differences. These procedures were conducted separately for nonsynonymous and synonymous mutations and revealed that the distribution of mutations was significantly different from random expectations.

Because genomes are known to include mutational hot-spots, we performed additional curation of multiply-mutated genes to find genes that were most likely targets of natural selection instead of accelerated mutational rates. We only considered genes as convergently evolved if they i) accumulated more mutations across replicate evolved lineages than any gene received in any randomized bootstrap trial or ii) were observed in at least half of all replicate evolved lineages. Additionally, reasoning that mutational hot spots should not be biased for or against silent mutations, the gene had to meet these criteria for the nonsynonymous mutation dataset but not also the synonymous mutation dataset. We applied these tests separately for each treatment group (e.g., pCO2 treatments for phytoplankton genomes, or pCO2 × partner for Alteromonas genomes). We further supported this curation step by performing a gene-by-gene linear model test in R to discover genes significantly more mutated under one treatment than another. Genes that met the curation criteria in one treatment group only and that also demonstrated statistical enrichment in that group were considered optimal candidates for genes under natural selection for a given treatment.

Genes remaining after this curation process were then analyzed for function. First, we binned genes into KEGG pathway groups using Over Representation Analysis (ORA) with the function enrichKEGG in clusterProfiler in R. The argument 'organism' was set to the KEGG organism codes 'pmi' and 'syn' for Prochlorococcus and Synechocystis, respectively, whereas it was to set to 'ko' (KEGG Orthology) for Synechococcus, T. oceanica, E. huxleyi and Alteromonas. KO identifiers to individual genes were assigned for each organism using the KEGG automatic annotation server website by the bi-directional best hit (BBH) method. P-values correspond to the comparison between the gene ratio (i.e., the number of genes that match that gene set divided by the number of genes in the ‘hit’ database) and the background ratio (i.e., the number of genes in the gene set divided by the total number of unique genes in the genome database). Under specific contrast comparisons for Alteromonas considering the mixed effect of co-cultures, KO identifiers were examined directly using the KEGG database website.

Second, the degree of association between mutating genes in EZ55 was analyzed using the Species Pairwise Association Analysis (SPAA) algorithm in R, reasoning that presence vs. absence of mutations in a given gene in a culture was mathematically comparable to presence vs. absence of species in an ecosystem. The mutation data was transformed into a binary presence vs. absence matrix, with genes having at least 1 mutation in each lineage coded as 1 and those without mutations coded as 0. The logistic regression function was used to calculate the odds ratio for mutations in each gene as described in, using cutoffs of 0.5 and 0.9 for the minimum and maximum mutation frequencies, respectively. SPAA estimated Spearman's rank correlation coefficient for each gene pair, measuring the monotonic predictive relationship, i.e., the likelihood of a mutation in one gene predicting a mutation in the other. After manually examining the data, we further filtered the pairwise predictions to consider only the most positively and negatively correlated pairs (coefficients between 0.5 to 0.9 for positive relationships and -0.3 to -0.2 for negative relationships). The curated matrix was imported in Cytoscape 3.8.2 for network visualization. The network was analyzed as a directed graph to obtain the indegree and outdegree of each node, i.e., the number of other genes to which each gene is connected as a predictive target or source, respectively. The target gene of each set of genes was illustrated as an arrow pointed towards the target gene, and the gradient of the edge connecting genes was set to reflect either a positive or negative correlation.

Finally, we manually examined several genes of interest identified from these analyses using BLAST against the NCBI nr database to attempt to provide superior annotations for hypothetical proteins or other vaguely described gene products.

Copy number analysis: We used R scripts to extract coverage data from all breseq-assembled genomes. Estimated copy numbers for plasmids were initially obtained simply by dividing average plasmid coverage by average chromosome coverage. However, inspection of coverage maps for the Alteromonas plasmid revealed highly uneven coverage. Closer analysis revealed three plasmid regions with different patterns of coverage. One region was generally present at the same or greater coverage as the chromosome, one was often completely absent, and a third often existed at an intermediate coverage level. We reasoned that these differences may reflect homologous recombination events leading to the excision of parts of the plasmid and/or insertion of plasmid regions into the Alteromonas chromosome and therefore looked for homologous regions between the plasmid and chromosome sequences using a dot plot analysis via the D-GENIES web interface. This analysis revealed several points of homology, including one approximately 9kb sequence, thus suggesting two things. First, a reduced size plasmid likely evolved as a subpopulation in several lineages, with approximately 50kb deleted following an unknown but reproducible event. Second, an approximately 70kb section likely integrated into the chromosome of most lineages, and in several lineages became the only surviving portion of the plasmid, with no evidence remaining of subpopulations carrying the remaining 150kb of the plasmid. Because of these trends, we calculated average coverage of the plasmid at three different regions: 170kb to 180kb to estimate the "insertable" portion of the plasmid, 120kb to 130kb to estimate the full, free plasmids, and 40kb to 50kb to represent the possible reduced-size plasmid. We tested the effects of phytoplankton partner and pCO2 treatment on Alteromonas plasmid copy number using pairwise Mann-Whitney tests in R with manually calculated Holm-Bonferroni corrections for multiple tests. We also compared the likelihood of plasmid loss (defined as less than 1% coverage of plasmid region 120kb to 130kb) in Alteromonas paired with cyanobacteria vs. eukaryotic phytoplankton using Fisher's exact test in R.

Examination of coverage maps further revealed the presence of elevated copy numbers for some chromosomal regions, suggesting the presence of duplications. We therefore investigated these more closely using R scripts, retrieving any regions where the average coverage was 5× the coverage (or 5× the standard deviation of coverage for low-coverage eukaryote genomes) across the chromosome. We excluded from further analysis all duplications in eukaryotic genomes where repetitive elements led to tens of thousands of qualifying duplications, as well as duplications falling outside of coding sequences or cis-acting promoter regions. After this process, only one duplication remained of interest: a promoter region mutation in the apolipoprotein N-acyltransferase gene in Prochlorococcus MIT9312, clearly visible in coverage maps and present in most lineages with up to 120× duplication.