The signs of adaptive mutations identified in the chloroplast genome of the algae endosymbiont of Baikal sponge.

Sampling and sequencing

The technology for collection of sponges by scuba divers, with the use of labelled transects, has developed since the detection of sponge disease and is described by Khanaev et al. [2018]. The samples collected by divers were immediately placed in containers with Baikal water over ice and transported to the lab, maintaining a constant water temperature.

The three samples of freshwater sponge L. baicalensis were collected in June 2016 in the Bol'shiye Koty area (51° 90´ 69 N´´, 105° 07´ 05 E´´) at a depth of 10 m within the same transect. One sample was obtained from a sponge that was healthy in appearance (exhibiting a green colour), one sample was taken from a diseased sponge and one from dead rotten sponge tissues.

Illumina pair-end reads were obtained by DNA metagenome sequencing at Novogene Inc. (Illumina PE 150), for all three samples; in addition, they were processed by a conventional bioinformatics pipeline, which included the filtering of sequencing errors and the assembly of contigs. The extraction and sequencing of RNA samples was also performed at Novogene Inc. This technique was possible only for healthy and diseased sponge tissues; not enough RNA was extracted from the rotten tissues.

Assembly and annotation of the chloroplast genome fragment

Chloroplast species are abundant enough in the genetic material available from metagenomic sequencing of sponge tissues. A template-based assembly is suitable for obtaining the sequence of a chloroplast genome or at least a substantial fragment of the genome. The sequence of the chloroplast genome of Choricystis parasitica (NC_025539) was used as a template for assembly and comparative analysis. This genome is a circular DNA with a length of 94206 base pairs.

To get a refined template from the scaffolds, obtained by de novo assembly of each separate metagenomic sample, scaffolds were selected which had a similarity to the template genome. Then, from cleaned reads, from both DNA and RNA sequencing, reads were selected which had a similarity to any of the selected scaffolds. The reads selected in all samples were merged to a single volume, and both reads in a pair were treated as unpaired. The ‘lightweight’ assembler Inchworm from the TrinityRnaSeq package was then applied to the volume of filtered reads, adjusted to the maximal size of k-mer (K = 31) and minimal sensitivity to sequencing errors. The contigs obtained after the second assembly were compared with the reference genome. This made it possible to select a single contig of 55638 bp in length, which was with confidence identified as a fragment of the chloroplast genome. The stages of the pipeline are listed in Table 2, to explain the data flow at each stage, and to provide the software specifications.

To check the obtained result, steps 3 and 5 were repeated with other settings: refined templates were obtained by running both blastn with an evalue of 1e-8 instead of 1e-40, and tblastx; the Inchworm assembly was run with default settings instead of the setup with the lowest sensitivity. The assembled contigs were not completely identical to the contig selected as a result, but the difference in any case was fewer than 10 substitutions.

Open reading frames were identified in the obtained contig and most of the identified proteins were annotated following the annotations of the reference genome. TrnaSCAN software was used to identify 13 transport RNA genes in the putative chloroplast sequence, and the locations of 18S rRNA and 23S rRNA were identified by direct alignment with reference rRNAs.

Identification of polymorphisms

In order to distinguish the polymorphic sites from traces of sequencing errors, their selection in the genome was implemented following a previously published approach [Truong et al., 2017]. Each RNA and DNA sample was represented as pair-end reads and separately aligned to the assembled fragment of the chloroplast genome using Bowtie2. The alignments were then processed using the Samtools 1.7 software pipeline with conventional settings, and the indexed archives of alignments were used as the input of the algorithm for the identification of polymorphic sites.

When describing the algorithm, s represents each position on the alignment of the reads, N_s is defined as the total number of reads and T_s is defined as the number of reads supporting the most abundant allele. Given the sequencing error rate E, the non-polymorphic null hypothesis was rejected if the probability that the number of reads equivalent to N_s − T_s coming from the non-dominant allele was < α = 0.05. This value was estimated using the probability mass function of a binomial distribution with N_s trials and a successful rate of 1 − E. The error rate was set to 0.01 for Illumina sequencing. Bases with a quality below 30 were removed, and reads with an average identity to the reference below 99% were ignored before applying the statistical test. Failing to reject the null hypothesis reflected the absence of alternative alleles or the inability to distinguish between low-coverage potential alternative alleles and sequencing noise.

Therefore, the number of polymorphic sites can be counted for each gene. Another property of each gene is the number of polymorphic sites where the count of alternative alleles is higher than the count of the dominant allele (T_s < N_s − T_s). This property could be used, in addition to the fraction of mutations in each gene, to detect the intensity of mutations in the DNA and RNA obtained in the samples.

Phylogenetic analysis

The chloroplast genome sequences of Picocystis salinarum (NC_024828), Myrmecia israelensis (KM462861), Botryococcus braunii (KM462884), Coccomyxa subellipsoidea (NC_015084), Hydrodictyon reticulatum (NC_034655), Mychonastes jurisii (NC_028579) and Chlorella vulgaris (NC_001865) were used to reconstruct the phylogenetic tree for the 16S ribosomal RNA (rrs gene) and ATP synthase subunit beta (atpB gene).

The 16S rRNA reference sequences from the Greengenes gg_13_7 database were used to clarify the relations between subspecies of the algae in Lake Baikal. These sequences were identified in the studies of Belikov et al. [2019] and Feranchuk et al. [2018] as the 16S rRNA templates which are extensively represented in the sponge microbiomes collected in different locations of the lake. Both templates, 4365343 (gb|JF495289.1) and 1118847 (gb|GU936925.1), were directed to the C. parasitica rRNA (SAG:17.98, gb|KM462878) as the closest annotated species, by a blast search at NCBI site (98% and 96%).

The nucleotide sequences of the selected genes were aligned using Mafft. A straightforward distance-based approach was used for phylogenetic analysis of the 16S rRNA and atpB genes of the chloroplasts of alga species from the Chlorophyta division. These trees were constructed using Fastme, with a neighbour-joining method to select tree topology and the Jukes–Cantor measure to calculate the distances between genes. A maximum-likelihood algorithm was selected to be used for a tree which was based on 16S rRNA fragments of the Choricystis clade, to consider the contribution of the nucleotide substitution model to branch lengths in the constructed trees. A comparison of 16S rRNA genes for the selected strains was performed using PhyML with the default parameters (--model HKY85 -d nt).

Remarks about the confidence of the assumptions

Indirect signs of adaptive mutations were detected in the integral distribution of nucleotides in polymorphic sites of the chloroplast genome. The polymorphisms detected in the samples could occur due to sequencing errors or to custom variations in the populations of alga substrains in closely located sponges, not only due to an adaptation to stress.

The samples were collected at the same place, and the alga strain was expected to be the same in all three samples. But micro-populations of algae can anyway be separated in any sponge, and their development can anyway change in response to stress. So, signs of an acceleration of mutations can be demonstrated only at a qualitative level. That is, assuming that there are several substrains of algae in the three samples, which are considered as a single environment, from the distribution of nucleotides at polymorphic sites one can estimate the relative abundance of dominant and alternative substrains in each sample. The observed difference in these abundances can be naturally explained, if the hypothesis about acceleration of mutations is assumed. It can also be explained by another reason, or just accepted as a random case, but to exclude the assumption about adaptive mutations would be a kind of ‘reduction’, an unnecessary simplification of the reality.

This approach should be treated not as a confident explanation of the observed event, but rather as a proof of the possibility. But even an approved possibility of extremely high mutation rates is of value, to clarify the ‘enigmas’ which are noticed in the evolution of algae, and to suggest the ways in which the community of sponges could save itself in a time of crisis.

The questions which should be answered to specify the degree of confidence in the presented results are: how adequate is the assembly of the chloroplast genome? To what extent is the estimated proportion of alternative alleles in polymorphic positions caused by an adaptation, compared with reasons like the natural diversity of genotypes, and with the biases introduced at the experimental and calculation stages?

The fragment of the genome obtained in the assembly is 55683 nucleotides in length, compared to 94206 nucleotides in the circular DNA of the C. parasitica chloroplast. An auxiliary run of the assembly pipeline with different settings showed that the fragment considered as the final result does not include at least 5000 nucleotides. But in all cases, the differences were fewer than 10 nucleotides in the whole fragment.

The order of annotated genes in the reference genome corresponded to the order of annotated open reading frames in the assembled fragment, with a difference in a single event of reordering, such that the orientation of the segment from rpoB to tufA genes was reversed. The segment in the reference genome between the rrs and tufA genes was missing in the assembled fragment; instead, the psaB gene from the middle of that segment followed just after the rrs gene. The homology between nucleotide sequences of genes was from 80% up to 98%.

The average proportion of polymorphic positions in all genes was 98.4%, which is comparable with the previously reported 97.8% for bacterial genomes [Truong et al., 2017]. These values are also consistent with the results reported by Feranchuk et al. [2018]. It was shown that on limiting the selection of 16S rRNA fragments of the chloroplast, sequenced by 454 technology, to 95% identity, variation was at most 2% between any of the fragments.

Coverage values for some coding frames were not sufficient to evaluate the distribution properties of alternative alleles. For all the annotated tRNA genes, coverage was at an adequate level, but no polymorphic positions were detected in any of these genes.

However, for 42 of the 56 annotated genes, the coverage was sufficient to detect polymorphisms with adequate significance. Within the purposes of the research, acceleration of mutations may be detected in the distribution of allele frequencies in polymorphisms only at a qualitative level. For qualitative analysis, for each gene and each sample, a coverage of reads was obtained where any of the alternative nucleotides was substituted in a polymorphic position, and those positions were counted for each gene where the coverage of reads with a dominant allele was less than half of the total coverage. These counts were then used to evaluate the abundance of a dominant strain, relative to minor substrains. The fractions of these numbers for each gene are comparable to the fractions for the whole genome, with acceptable dispersion.

The annotation of the alga species under study, with respect to close species of algae, including the species of sponge symbiotic algae in other parts of the lake, was also refined. This made it possible to compare the variations observed between strains in the three samples collected at the same location with those in strains of algae from sponge samples collected in other locations of the lake. The phylogenetic analysis explained above is reported in the last subsection of the results.

Distribution of polymorphic sites and its interpretation

The fraction of alternative alleles in polymorphic positions is shown in Figure 1, for a whole set of genes and for each separate gene. Within this fraction, the part is separated where a crucially low (< 50%) abundance of baseline nucleotides was detected in the polymorphic positions.

Figure 1. Relative proportion of polymorphic sites in the chloroplast genome of sponge symbionts, for the five samples studied.

The results for each of the annotated genes are shown on the right. The left column presents the integrated results for each sample. The proportions of polymorphic sites and sites with high levels of alternative alleles (‘mutations’) are shown as pie charts, relative to the total number of polymorphic sites in all samples. The proportion of sites which are polymorphic in some other samples, but not in the given sample, are shown in light blue. The legend on the left shows the colour scheme used to represent three types of site. The circle radius represents the total number of sequencing reads aligned to the gene segment and used to identify polymorphic sites. The scale of the circle radii is transformed for better appearance, to compensate for the high variation in the number of aligned reads.

The separation into three fractions, which is shown in Figure 1, allows one to guess, at a qualitative level, the abundance of the dominant strain, relative to the second most dominant strain and minor strains. Also, comparison of distributions for DNA and RNA allows one to guess how intensive is the development in each of the three separated groups, assuming that RNA synthesis is a sign of any development. This qualitative interpretation is shown in Figure 2.

Figure 2. Qualitatively assessed abundance of chloroplast substrains in three samples of sponge.

When a sponge is healthy, the major strain dominates, but mutations do not affect the rate of development. This corresponds to a minimal fraction of polymorphisms, and the dominant allele is anyway most frequent in any of the polymorphisms. When a sponge suffers from the disease, it causes, by assumption, an increase in mutation rates and, in turn, an increased abundance and re-distribution of minor strains. This corresponds to an increased fraction of alternative alleles in the genome and, to a lesser extent, in the transcriptome. When a sponge is completely destroyed, just the ‘white noise’ of mutations is found in the genetic material which is left from its chloroplasts.

The precise quantitative estimates of these abundances, to our best knowledge, by no means can be obtained from the available volume of data. The abundance of substrains in the three sponges before the development of disease also cannot be reconstructed.

However, for the sample from the diseased sponge, the fraction of alternative alleles was much higher than in the other two samples. So, as a hypothesis which expands the neutral model, it can be suggested that in the alga cells of a diseased sponge, mutations become accelerated to a rate much higher than any rate compatible with the consistence of metabolic relations. In other words, the living cells in the diseased but alive tissue are desperately trying to survive. Subpopulations of mutated strains develop slowly and die earlier than a dominant lineage – a natural consequence of so many fast mutations in the genome.

But the proposed hypothesis means that the observed acceleration of mutations was a determined event, and some mechanism encoded in the genome triggers this event. But for what reason did this mechanism arise and was it kept conserved in evolution? It is possible just to guess that reason, and several suggestions for that guess are proposed the Discussion section.

Comparative analysis and taxonomic annotation of symbiotic algae

The phylogenetic trees for the two selected genes, 16S rRNA and ATP synthase beta (Figure 3 A, B), in general confirm the conventional relationship between Chlorophyta algae [Lemieux et al., 2014; Lemieux et al., 2015]. The trees in Figure 3 support the assumption that the symbiotic algae of the sponge L. baicalensis are close in taxonomy to the genus Choricystis.

Figure 3.

Phylogenetic trees for two chloroplast genes, ATP synthase subunit beta (A) and 16S rRNA (B, C). The trees in A and B are constructed for species within the Chlorophyta division; tree C is constructed for species within the genus Choricystis. Bars located at the node for the studied chloroplast genome in trees A and B represent the relative number of polymorphic positions in all five samples, at a 1 : 1 ratio.

The bars in Figure 3 A and B, which show the proportion of polymorphic positions in the genes of symbiotic algae in metagenomic samples, are comparable in size with the scale of distances between genera. The timescale which separates the origins of the close genera in Figure 3 is much larger than the timescale which could adequately interpret the separation of the chloroplast strains detected in the metagenome. So, the observed relation between timescales may be interpreted as being caused by adaptation events, and in that case it can be used to estimate the order of magnitude of the acceleration rate.

Tree C in Figure 3 demonstrates in a more precise way the relation between the species and subspecies of symbiotic algae. The tree is composed an rrs gene for the reference species of C. parasitica, the two abundant chloroplast rRNA templates for the sponges from other parts of the lake, and three consensus sequences for the three sponges under study. From these three sequences, one represents the 16S rRNA gene from the reference assembly, and two sequences represent specifically the sample of diseased sponge, as a DNA and as an RNA. The consensus sequences for the latter two datasets were created from the distribution of polymorphic positions.

The chart in Figure 3C illustrates in more detail the estimates of acceleration value provided below. The following assumptions can be introduced: the fraction of polymorphic positions in the whole chloroplast genome is 1.6%; the age of Baikal is about 20 million years, but let the separation of the alga species in Baikal from the rest of the Choricystis lineage be about 1 million years ago, and the genetic distance from the alga species under study to C. parasitica be about 20%; the fraction of mutations in the diseased sample is assumed to be 12.5% in polymorphic positions and 0.2% in the whole genome; the disease was developed in one season.

Under these assumptions, the relative increase of mutation rate in diseased sponge is estimated to be 10000. Such high acceleration is obviously incompatible with life but, as a comparable precedent, the mutations in cells of some tumours are accelerated 200 times [Bielas et al., 2006] and even up to 10000 times [Berger et al., 2012].