In recent decades, the cyanobacterial species Cylindrospermopsis raciborskii has received growing attention because of its increasing geographical distribution and the frequency of toxic blooms (Padisák, 1997; Briand et al., 2002; Kokociński and Soininen, 2012; Sinha et al., 2012; Sukenik et al., 2012). The invasive potential of C. raciborskii is constantly reported and mostly attributed to its ability to tolerate a wide range of climatic conditions and global warming (Neilan et al., 2003; Wiedner et al., 2007; Antunes et al., 2015; Kokociński et al., 2017). Toxic bloom episodes are considered harmful to ecosystems due to their impact on water quality, which decreases the efficiency of energy transfer from primary producers to primary consumers (Lürling and Roessink, 2006) and poses potential health risks to both animals and humans (Humpage et al., 2000; Murphy and Thomas, 2001; Chellappa et al., 2008). In addition to the production of cyanotoxins, C. raciborskii strains may also synthesize other biologically active natural compounds with different activities. An important step in understanding the opportunistic behavior of C. raciborskii is advancing our knowledge of morphological and genetic variation in this cyanobacterial species.

The high phenotypic plasticity of the species C. raciborskii (Hawkins et al., 2001; Komárek and Komárková, 2003; Chonudomkul et al., 2004; Moore et al., 2004; Piccini et al., 2011; Bonilla et al., 2012) affects identification accuracy. The major problems in identifying species arise because the morphology of the trichome varies during population growth, and under certain environmental conditions the differentiating akinete and heterocyte cells are absent. The morphological variability of C. raciborskii strains has also been documented under laboratory culture conditions. This morphological plasticity often leads to misidentification as the morphologically, phylogenetically, and ecologically related genus Raphidiopsis. The genera Cylindrospermopsis and Raphidiopsis belong to the same order Nostocales, and family Aphanizomenonaceae (Komárek et al., 2014), and members of both genera are morphologically very similar and can co-occur in nature. The presence or absence of heterocytes is commonly used as a diacritical morphological feature differentiating these genera. Heterocytes are morphologically distinct cells that develop from vegetative cells and are specialized for nitrogen fixation. In Cylindrospermopsis, heterocytes are always terminal at one or both ends of trichomes, and the absence of this specialized cell is not unusual (Komárek and Anagnostidis, 1989). Several attempts were made to solve the taxonomic relationship between Cylindrospermopsis and Raphidiopsis genera using molecular and phylogenetic approaches (Gugger et al., 2005; Stucken et al., 2009; Alster et al., 2010; Wu et al., 2011; Li et al., 2016). Recently, a proposal to unify Cylindrospermopsis and Raphidiopsis genera under the single genus Rhaphidiopsis was published (Aguillera et al., 2018). The name Rhaphidiopsis was retained as the genus name since it has nomenclatural priority over the name Cylindrospermopsis.

Although the morphological variability in C. raciborskii has been well-documented, only few comparative genomics studies of this species are available (Stucken et al., 2010; Sinha et al., 2014). These few genomic studies are concentrated in Australian C. raciborskii strains, and focused mainly on the toxigenic aspects of the strains. Currently, five draft genomes of C. raciborskii and one of Raphidiopsis brookii are available in public databases. This number of genomes is the minimum recommended to perform a pan-genome analysis (Rouli et al., 2015). A microbial pan-genome is the combination of genes present in all strains (the core genome) and genes that are present in only one or a few strains (known as a dispensable or flexible/accessory genome) (Tettelin et al., 2005). Pan-genome analyses provide a framework for estimating the genomic diversity of a bacterial species but are also helpful to predict, via extrapolation, how many additional genome sequences would be necessary to characterize the entire pan-genome or gene repertoire (Rouli et al., 2015). To shed some light on the knowledge of genetic variation among C. raciborskii strains, we combined all available C. raciborskii genomes with two new genome sequences to perform the first large-scale comparative genomics study for this species, with representatives from different world regions. More specifically, the goals were (1) to define the core genome and accessory genome of C. raciborskii, and (2) to identify toxin, heterocyte differentiation and nitrogen fixation genes, as well as other genomic features associated with the synthesis of natural products.

In this study, we used comparative genome analysis to investigate the C. raciborskii pan-genome based on eight genome sequences to comprehend the genetic landscape of this species. The genome dynamics of the C. raciborskii strains showed that the pan-genome size increased steadily without reaching a plateau, while the core genome appeared near saturation, suggesting that further sampling of C. raciborskii genomes is likely to cause a closed state in the core genome. The Cylindrospermopsis genomes from different origins generally showed a significant resemblance, with similar patterns in gene functions. Notably, in all eight genomes, it was only possible to assign functions to approximately 40% of protein sequences, which emphasizes the lack of knowledge of cyanobacterial proteins and/or a more cursory annotation of proteins in this phylum.

Subsequently, increasingly refined examinations showed that specific and recurrent genomic variations occur among the seven strains of C. raciborskii and the single R. brookii D9 strain. The genomic features of these strains appear to be influenced by their geographic origin, since evolutionary analysis differentiated the South American strains (CENA302, CENA303, ITEP-A1, MVCC14, and D9) from the Australian and Singaporean strains (CS-505, CS-508, and CR12). However, the lack of European and African genome sequences of this species did not allow a more precise conclusion about the influence of the environment on their genomic plasticity. Previous studies on phylogeny of several C. raciborskii genes already observed biogeographic variation correlated to genetic variability (Neilan et al., 2003; Gugger et al., 2005; Moreira et al., 2013, 2015; Antunes et al., 2015).

Despite the separation between South American and non-South American genomes, neither the genotypes, subsystem annotation nor enrichment analysis appeared sufficient to explain this separation. By examining the subsystem annotations of the two largest regions of genomic divergence in the BLAST atlas, we once again found that most gene functions could not be assigned (Supplementary Figure S7). While little can be inferred from region (B), region (A) showed that almost half of the genes are involved in N-linked glycosylation, which is a process of protein modification. These processes are not yet fully understood but appear to be frequently present in the three domains of life and control a wide range of cellular processes, such as signal transduction, protein folding, targeting, stability, cell–cell and virus–cell interactions, and host immune responses (Nothaft and Szymanski, 2013). Therefore, although more information is necessary before we can fully determine the differences between these groups (South American and non-South American), the results of our analyses indicate that they have a common primary metabolism and minor differences in very important cell processes, such as cell regulation, biosynthesis, and transport, which allows the exploration of a wide range of environmental conditions.

The identification of significant associations of GO terms between South American and non-South American genomes was not deeply accessed through the GO enrichment analysis. Low resolution may occur in this type of analysis, since functional enrichment analysis is directly associated with sequence annotation level depth (Wang et al., 2011; Glass and Girvan, 2015). The sequence set analyzed in this study gave insufficient functional categories, which influenced the results of this type of functional enrichment analysis. In other words there was an enrichment bias.

Since the group containing CYN-producing strains (Australian CS-505 and Singaporean CR12) was also within the non-South American group (CS-505, CS-508, and CR12), it is possible to evaluate the results of the enrichment analyses of these two groups together. The percentage of genes annotated as involved in the biosynthesis of polyketides might be influenced by the cyr gene cluster, which contains polyketide synthase genes and is present only in Australian CS-505 and CR12 strains. Polyketides boast a wealth of medicinally important activities, including antibiotic, anticancer, antifungal, antiparasitic, and immunosuppressive properties (Staunton and Weissman, 2001). In addition, the Cellular Component ontology showed that almost half of the singletons successfully annotated in the non-South American genomes are involved in the cell membrane. Therefore, these singletons could be directly associated with the adaptation of cellular machinery to particular environmental conditions. The strains not only synthesize secondary metabolites and antibiotics found in the South American genomes, but also show a high percentage of genes involved in cell transport, membrane, and regulation of gene expression.

All genomes of STX-producing strains are from South America, and consequently, these two groups were also analyzed together. Similar to the singletons from the non-South American genomes, the South American singletons were also related to adaptation to a particular environment. However, whereas the first were mainly associated with changes in cellular machinery, the latter were also linked to differences in cellular metabolism and recognition of environmental stimuli and chemicals. This hypothesis is sustained by the fact that the South American genomes presented singletons that participate in water-soluble vitamin biosynthetic processes, similar to the pattern found in cyanobacteria typical of lakes and rivers (Wetzel, 2001), as well as genes related to post-translational/co-translational cofactors, cell movement, and cellular cycle.

The strains CENA303 and CS-508 did not show either the cyr or sxt gene clusters, and therefore they formed their own group (non-toxic strains). The enrichment analyzes of this group resulted in only four annotated genes, and thus not much information can be inferred. This result illustrates the current lack of information on cyanobacteria and the difficulty of functional annotation in this phylum.

This comparative genomic analysis showed that the five South American (Brazil and Uruguay) genomes are slightly smaller and more conserved than the non-South American (Australia and Asia) genomes, suggesting that genomes from South America underwent gene loss events. For example, the Brazilian strains CENA303 and D9 have smaller numbers of coding genes than the other strains and have substantial differences in their genomes compared with the Australian CS-505 and Singaporean CR12. Signs of genome reduction have already been reported for the R. brookii D9 strain (Stucken et al., 2010). Genome reduction is an evolutionary strategy strongly related to genome efficiency and selective pressure, and can distinctly affect various cellular processes, directing selection for the maintenance of specific genes. Usually, events of gene loss and genome shrinkage are related to the accessory genome, i.e., genes that are not conserved among genomes, and do not involve crucial housekeeping functions (Bolotin and Hershberg, 2016; Paul et al., 2016). Housekeeping genes are conserved, fitted in the core genome, and they are considered good molecular markers for phylogeny (Wu and Eisen, 2008; Calteau et al., 2014; Uyeda et al., 2016). In this study, the similarity matrix indicated that the presence of additional genotype traits results in more distantly related genomes. Accessory genes contribute to increase the distinction among genomes, however, they are not conserved and they are not good molecular markers for phylogeny (Sinha et al., 2014).

Evolutionary relationships among the eight strains elucidated by this comparative genomic study demonstrated that the strain deposited in the GenBank as R. brookii D9 is highly related to C. raciborskii strains, as already evidenced in previous reports (Stucken et al., 2009, 2010; Hoff-Risseti et al., 2013; Aguillera et al., 2018). Previous studies identified a different morphotype within the culture of C. raciborskii T2 (also known as SPC338) as Raphidiopsis likely due to the absence of heterocytes in their filaments (Castro et al., 2004; Stucken et al., 2009). However, the lack of this specialized cell is not unusual in Cylindrospermopsis strains (Komárek and Anagnostidis, 1989). Misidentification between Cylindrospermopsis and Raphidiopsis strains is common, since they are morphologically and genetically close; specifically, both genera are included in the order Nostocales and the family Aphanizomenonaceae (Komárek et al., 2014). Indeed, this study also showed that C. raciborskii CENA303 filaments do not contain heterocytes. Even after complete depletion of inorganic nitrogen sources in culture media, vegetative cells of CENA303 were incapable of initiating heterocyte differentiation. The CENA303 and D9 strains are phylogenetically close and do not possess nif genes, although they contain regulatory genes (hetR and hetN) involved in heterocyte differentiation (Herrero et al., 2013). The presence of these genes suggests that selective losses of the nif and hgl clusters have occurred in the CENA303 and D9 strains, which lost the capacity to differentiate heterocytes.

Given the metabolic potential for the production of natural products, this comparative genome study emphasized gene clusters involved in the synthesis of the antifungal hassallidin and the toxins cylindrospermopsin and saxitoxin. Much of the natural products synthesized by cyanobacteria are attributed to the versatility of NRPS and PKS enzymes that catalyze chain elongation from simple building blocks to create a diverse array of molecules (Welker and von Döhren, 2006). The extent and distribution of secondary metabolite synthesis in cyanobacteria appears to be correlated to genome size (Shih et al., 2013). The C. raciborskii genome is the smallest known genome in multicellular cyanobacteria (Stucken et al., 2010) and therefore a low frequency of NRPS and PKS enzyme systems for the production of secondary metabolites is expected, which was observed in this study.

Hassallidins are antifungal cyclic glycosylated lipopeptides that were first identified as being produced by cyanobacteria of the genus Hassalia (Neuhof et al., 2005, 2006). Other cyanobacteria were also shown to be producers of this class of peptides and the gene cluster in the Anabaena sp. SYKE748A strain was characterized as spanning 59 kbp and consisting of 26 ORFS (hasA-hasZ) (Vestola et al., 2014). The biosynthesis of hassallidins involves many genes, and the role of some of these in the biosynthetic pathway is not yet clear. While NRPS genes (hasO, hasN, hasV, and hasY) encode the nine modules responsible for the incorporation of amino acids into the hassallidin peptide backbone, the glycosylation genes (hasD, hasQ, hasT, and hasX) are then responsible for the incorporation of sugars, and correspondingly, the lipidation genes (hasG, hasH, and hasL) for the addition of lipids (Vestola et al., 2014). The Australian CS-505 strain is a known producer of hassallidins (Vestola et al., 2014) and the has gene cluster was recently described (Pancrace et al., 2017). Clusters showing similar genetic organization to that found in the CS-505 strain were identified in the genomes of the CENA303, CS-508, and CR12 strains. Predictions of the A-domains of NRPS genes indicate that compared with the Anabaena sp. SYKE748A strain, which incorporates a methylthreonine at position nine in the molecule, the Cylindrospermopsis strains incorporate a threonine. In addition, the number of glycosylation and lipidation genes varies among strains and may act in more than one step in hassallidin assembly. Therefore, in silico analyses indicated that the CENA303, CS-508, and CR12 strains are potential producers of new hassallidin variants. The IS insertion found inside the hasR gene in the CENA303 strain may irreversibly inactivate this gene. Nevertheless, HasR has been described as an acyltransferase that acts on the acetylation of hassallidin sugars (Vestola et al., 2014). Thus, its inactivation cannot affect the production of the molecule but merely the variants produced. As half of the strains included in the present study are potentially hassallidin producers, and the two others contain fragments of the gene cluster, this molecule appears to be widespread among C. raciborskii.

The shift in toxin production according to geographic origins had already been reported in C. raciborskii strains (Lagos et al., 1999; Piccini et al., 2011; Hoff-Risseti et al., 2013). In South America, C. raciborskii strains produce the neurotoxin saxitoxin and some of its derivatives, while in Australia and Asia, the isolates synthesize the alkaloid cylindrospermopsin. Both toxins are synthesized by the NRPS and PKS enzyme complexes (Kellmann et al., 2008; Mihali et al., 2008). The results of this study showed that the genomes of South American strains contain the sxt gene cluster, while Australian and Singaporean strains possessed the cyr gene cluster, thus corroborating information available in the literature. Furthermore, among the 26 genes proposed to be involved in STX production (Kellmann et al., 2008), 21 were found in the genomes of the STX-producing strains CENA302, D9, and MVCC14. Since the Brazilian strains CENA302 and D9 are known as STX producers (Stucken et al., 2010; Hoff-Risseti et al., 2013), the genes missing (sxtJ or sxtK, sxtV, sxtZ, sxtY, and sxtW) should be not essential for STX production or may be involved in generating toxin variants.

Vestiges of cyr genes were found in the genomes of South American STX-producing C. raciborski strains. However, the annotated gene fragments were different from the gene fragments amplified by PCR using specific primers in previous studies (Hoff-Risseti et al., 2013; Piccini et al., 2013). No vestiges of sxt genes were found in the genomes of the CYN-producing strains. The finding that STX-producing C. raciborskii strains also carry fragments of cyr genes suggests either a remnant or an otherwise ancestral intermediate of a functional CYN gene cluster (Hoff-Risseti et al., 2013).

This study provides new insights into the genotypic and phenotypic plasticity of the species C. raciborskii, including the presence of nitrogen-fixing and non-nitrogen-fixing strains. Both strains may co-occur in nature in large populations, such as those formed during bloom events. Environmental conditions usually found during blooms likely favor gene loss events found in the smaller genomes of non-nitrogen-fixing C. raciborskii strains, without compromising their survival. This comparative genomic analysis also successfully estimated the genetic diversity among C. raciborskii strains and the R. brookii D9 strain revealing a conserved core genome with major differences at the accessory genome levels. The biosynthetic capacity for producing the secondary metabolites of the investigated genomes was low, which is expected due to its small genome size.

MF conceived the study. CH-R, PS, and DA performed the experiments. VA, RP, DA, CH-R, PS, AV, and MF analyzed the data. All authors were involved in writing the paper and had final approval of the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.