Distribution and conservation of known secondary metabolite biosynthesis gene clusters in the genomes of geographically diverse Microcystis aeruginosa strains

Additional keywords: biogeography, cyanobacteria, cyanobactin, cyanotoxins, microviridin, non-ribosomal peptide, phylogeny, polyketide.

Microcystis aeruginosa is a cosmopolitan bloom-forming species of cyanobacteria notorious for its production of hepatotoxic microcystins. These small cyclic peptides are toxic in very low doses and pose a threat to human health, agriculture and freshwater ecosystems worldwide. Although microcystin is the most toxic metabolite produced by M. aeruginosa, recent genome sequencing efforts have revealed that a variety of other secondary metabolite (SM) pathways are also encoded within the M. aeruginosa pan-genome (Humbert et al. 2013; Fig. 1). Whether the corresponding biologically active compounds pose a threat to water quality is uncertain. However, solely focusing on microcystin when assessing water quality could prove to be an oversight.

Fig. 1.  Schematic of known secondary metabolite biosynthesis gene clusters present in Microcystis aeruginosa and other cyanobacteria. The above clusters (sourced from the Minimum Information about a Biosynthetic Gene cluster (MIBiG) repository) were used as reference clusters in the present study. The asterisk (*) indicates patBC genes (named after their association with the patellamide biosynthesis gene cluster) that are commonly associated with cyanobactin gene clusters but are not essential for cyanobactin biosynthesis (Sivonen et al. 2010). Diagram not to scale. Choi, L-2-carboxy-6-hydroxyoctahydroindole; NRPS, non-ribosomal peptide synthetase; PKS, polyketide synthase.

Interestingly, the ability to produce microcystins appears to be sporadically distributed across the species (Tillett et al. 2001), with toxic and non-toxic M. aeruginosa strains coexisting in many ecosystems. The apparent lack of phylogenetic linkage between microcystin production and traditional taxonomic markers (16S rDNA, phycocyanin, etc.; Tillett et al. 2001) initially hampered molecular differentiation of toxic and non-toxic M. aeruginosa strains. However, the discovery of the microcystin synthetase (mcy) gene cluster in PCC 7806 (Dittmann et al. 1997; Tillett et al. 2000) enabled the development of highly specific and sensitive polymerase chain reaction (PCR)-based methods for assessing the toxigenicity of bloom samples (e.g. Jungblut and Neilan 2006; Al-Tebrineh et al. 2010, 2011; Baker et al. 2013). These methods are now widely used by scientists and water quality authorities around the world.

Although it has long been recognised that M. aeruginosa is capable of producing a wide range of other SMs, including non-ribosomal peptides, hybrid polyketide–non-ribosomal peptides (PK-NRPs), cyanobactins and microviridins (Fig. 2), these compounds have received much less attention than the microcystins in terms of their effects on water quality (Janssen 2019). However, many have demonstrable effects on growth and viability in laboratory test organisms, such as mice and zebrafish, as well as in cell cultures, highlighting their potential to disrupt aquatic ecosystems. Furthermore, the biogeography of SM biosynthesis pathways in M. aeruginosa has not been investigated in detail, and it is unclear whether their distribution is scattered, as for microcystin, or segregated according to the geoclimatic origin of isolates. Characterisation of the distribution, composition and conservation of the various SM biosynthesis gene clusters in M. aeruginosa strains will provide insight into their potential SM profiles and is an important first step towards holistic molecular monitoring of freshwater blooms.

Fig. 2.  Examples of known secondary metabolites produced by Microcystis aeruginosa. Select examples of secondary metabolites previously isolated from M. aeruginosa. Aeruginosin 98, microcystin-LR, microginin FR1, anabaenopeptin A and microviridin B structures were appropriated from Huang and Zimba (2019). Micropeptin K139, aeruginosamide B, microcyclamide 7806B, piricyclamide 7005E1 structures were appropriated from Nishizawa et al. (2011); Leikoski et al. (2013); Portmann et al. (2008) and Leikoski et al. (2012) respectively. For an overview of the different structural variants of these compounds, see Huang and Zimba (2019).

Hybrid PK-NRPs produced by M. aeruginosa include the microcystins, aeruginosins and microginins. The microcystins are a structurally diverse family of non-ribosomal heptapeptides sharing a common monocyclic structure and a unique β-amino acid (Adda) side-chain (Botes et al. 1985). Microcystins inhibit eukaryotic serine/threonine protein phosphatases 1 and 2A (MacKintosh et al. 1990) and are toxic to a variety of animals and plants. Bioaccumulation of the toxin may also occur in resistant organisms (Ferrão-Filho and Kozlowsky-Suzuki 2011). In mammals, the primary target of microcystin is the liver, where disruption of the cytoskeleton causes hepatic haemorrhage and organ failure (Dawson 1998). Chronic subacute exposure to the toxin may additionally initiate and promote tumour development (Fujiki and Suganuma 2011). The gene cluster responsible for microcystin biosynthesis in M. aeruginosa PCC 7806 spans 55 kb and encodes 10 genes (mcyA–J), including three non-ribosomal peptide synthetases (NRPSs), one polyketide synthase (PKS), two hybrid NRPS-PKSs, three auxiliary enzymes and an ABC transporter (Tillett et al. 2000; Fig. 1). The mcy gene cluster has since been identified in several other cyanobacterial genera (Christiansen et al. 2003; Rouhiainen et al. 2004) and is highly conserved in terms of composition and sequence identity compared with other SM biosynthesis clusters.

Aeruginosins are linear non-ribosomal tetrapeptides characterised by a central L-2-carboxy-6-hydroxyoctahydroindole (L-Choi) moiety, a C-terminal arginine derivative (e.g. argininal, argininol, agmatine, 1-amidino-2-ethoxy-3-aminopiperidine) and an N-terminal hydroxyphenyl lactic acid or phenyl lactic acid group (Ersmark et al. 2008; Scherer et al. 2016). Most aeruginosins are potent inhibitors of serine proteases, including the human blood coagulation cascade enzyme thrombin and the digestive enzyme trypsin (Ersmark et al. 2008). Certain variants may also be highly toxic to invertebrates, including crustaceans (Scherer et al. 2016). The gene cluster responsible for aeruginosin biosynthesis in M. aeruginosa PCC 7806 spans 25 kb and encodes 14 genes (aerA–N), including three NRPSs, one hybrid NRPS-PKS, six putative Choi biosynthesis genes, a halogenase, an ABC transporter, an oxidoreductase (OxRed) and several cryptic genes (Ishida et al. 2009; Fig. 1).

The microginins are a family of linear non-ribosomal lipopeptides that contain four to six amino acids, including one or more tyrosines at the C-terminus and the characteristic N-terminal side chain 3-amino-2-hydroxy decanoic acid (Ahda; Strangman and Wright 2016). Microginins inhibit angiotensin-converting enzyme (ACE; Okino et al. 1993), a central component of the renin–angiotensin system that controls blood pressure. The gene cluster responsible for microginin biosynthesis in Planktothrix NIVA-CYA98 spans 21 kb and encodes four genes (micA–E), including a PKS-NRPS, two NRPSs and an ABC transporter (Rounge et al. 2009; Fig. 1).

Non-ribosomal peptides produced by M. aeruginosa include the cyanopeptolins and anabaenopeptins. Cyanopeptolins are a family of cyclic depsipeptides characterised by the presence of 3-amino-6-hydroxy-2-piperidone (Ahp; Martin et al. 1993). Members of the cyanopeptolin family have been shown to inhibit several serine proteases. The gene cluster responsible for cyanopeptolin biosynthesis in Microcystis NIVA-CYA172/5 spans 31 kb and encodes seven genes (mcnA–G), including four NRPSs, a halogenase and an ABC transporter (Tooming-Klunderud et al. 2007). The micropeptin cluster in M. aeruginosa K-139 and PCC 7806 is highly similar in structure and composition to the cyanopeptolin cluster, although it lacks the halogenase (mcnD) and the cryptic enzyme (mcnG; Nishizawa et al. 2011; Fig. 1).

Anabaenopeptins are a family of diverse cyclic hexapeptides characterised by several conserved motifs, including a ureido bond, N-methylation in Position 5, and D-Lys in Position 2 (Christiansen et al. 2011). Anabaenopeptins are potent inhibitors of carboxypeptidase-A (CPA) and thrombin-activatable fibrinolysis inhibitor (TAFIa), enzymes involved in blood coagulation (Schreuder et al. 2016). The anabaenopeptin (apnA–E) gene cluster in Planktothrix agardhii CYA126/8 spans 24 kb and encodes four NRPSs and an ABC transporter (Christiansen et al. 2011; Fig. 1). In some cyanobacterial genomes, the apn gene cluster is colocated with the hph gene cluster, which is believed to provide additional substrates (e.g. L-homotyrosine and L-homophenylalanine) for the biosynthesis of anabaenopeptins (Lima et al. 2017).

M. aeruginosa produces a variety of modified ribosomal peptides known as cyanobactins, which are synthesised by the proteolytic cleavage and cyclisation of precursor peptides. Structural diversity within this family of compounds is achieved by heterocyclisation and oxidation of amino acids, or their post-translational prenylation (Sivonen et al. 2010). Cyanobactins exhibit a diverse range of bioactivities, including cytotoxic, antimalarial and allelopathic activities. Although quite diverse in their size (6–20 kb), composition and arrangement, cyanobactin biosynthesis gene clusters always encode two proteases (homologues of patellamide biosynthesis genes, PatA and PatG) that cleave and cyclise the precursor peptide as well as proteins participating in post-translational modifications (Sivonen et al. 2010). Some cyanobactin biosynthesis clusters (e.g. microcyclamide, mcaA–G; and aeruginosamide, ageA–G) also contain heterocyclase enzymes, whereas others (e.g. piricyclamides, pirA-G) do not (Leikoski et al. 2013).

Another class of post-translationally modified peptides produced by M. aeruginosa are the microviridins, which are tricyclic depsipeptides with an unusual cage-like architecture composed of two lactone rings and one lactam ring (Ishitsuka et al. 1990). Microviridins are potent inhibitors of serine proteases, including elastase, chymotrypsin and trypsin (Murakami et al. 1997). A particular variant, microviridin J, has been shown to be toxic to zooplankton (Rohrlack et al. 2004). The gene cluster responsible for microviridin B biosynthesis in M. aeruginosa NIES-843 contains five genes (mdnA–E), encoding the leader peptide (MdnA), two dedicated ATP-grasp ligases (MdnB and MdnC) a Gcn5-related N-acetyltransferase (GNAT)-type acetyltransferase (MdnD) and an ABC transporter (MdnE; Ziemert et al. 2008a; Philmus et al. 2008; Weiz et al. 2011; Fig. 1). MdnA may be distally encoded in other strains (Gatte-Picchi et al. 2014).

Although the aforementioned compounds have been isolated from several M. aeruginosa cultures and bloom samples, the corresponding biochemical pathways and biosynthesis gene clusters have only been examined in a few cyanobacterial strains, in some cases from distantly related genera. In order to better understand the corresponding genetic diversity and evolutionary history of these pathways in M. aeruginosa, we examined the distribution, composition and conservation of known SM biosynthesis gene clusters in 27 M. aeruginosa strains isolated from different geoclimatic regions. A comprehensive multilocus phylogenetic analysis was also undertaken to understand the biogeography of SM production in this taxon. The results are discussed with regard to predicting SM profiles in this species and designing future molecular diagnostic tests for potentially harmful strains.

Climate zones were allocated to M. aeruginosa strains based on an updated version of the Köppen–Geiger climate classification system (Kottek et al. 2006) accessed through the kgc package (C. Bryant, N. R. Wheeler, F. Rubel, and R. H. French, see https://CRAN.R-project.org/package=kgc, accessed September 2018), which classifies regions of the Earth by their relative heat and humidity through the year. Climate zones are divided into five main groups: A (tropical), B (dry), C (temperate), D (continental) and E (polar). The groups are then subdivided based on seasonal precipitation and temperature. Geographical coordinates used to determine Köppen–Geiger climate zones were based on reports in the literature for each isolate.

M. aeuriginosa genome assemblies with RefSeq (Haft et al. 2018) membership (n = 27) were downloaded from RefSeq and GenBank (Benson et al. 2018) between 2 February 2017 and 23 June 2018 in fna (used in average nucleotide identity, ANI, and phylogenetic analyses) and gbff (used in the prediction of secondary metabolite biosynthesis gene clusters) format respectively. General genome features were also obtained from GenBank. Of the two available NIES-298 whole-genome sequences, the assembly with the higher genome size and gene count (BEYQ01) was used.

ANI values were calculated according to the orthologous average nucleotide identity (OrthoANI) algorithm (Lee et al. 2016) implemented with USEARCH (OrthoANIu; Yoon et al. 2017).

Putative SM biosynthesis gene clusters encoded within the 27 M. aeruginosa genomes were identified and analysed using Antibiotics and Secondary Metabolite Analysis Shell (antiSMASH, ver. 4.10, see https://antismash.secondarymetabolites.org, accessed Sept 2018; Blin et al. 2017), with ClusterFinder enabled (minimum cluster size 4; minimum probability 60%). Results were verified by comparing the inferred peptide sequence identity (ID; minimum 50%), coverage (minimum 50%) and organisation of genes within putative clusters to 1481 previously characterised SM biosynthesis clusters in the Minimum Information about a Biosynthetic Gene cluster (MIBiG) dataset (Medema et al. 2015; Fig. 1), GenBank (by BLASTP and TBLASTN; Altschul et al. 1990) and the literature.

Phylogenetic analysis of the bacterial core-gene set was undertaken using the up-to-date bacterial core gene (UBCG) set and pipeline for phylogenomic tree reconstruction programme, UBCG (ver. 3.0, see https://www.ezbiocloud.net/tools/ubcg, accessed September 2018; Na et al. 2018), configured with default settings (Jones, Taylor and Thorton model; Jones et al. 1992): where varying rates of evolution across sites were modelled using the ‘fixed number of rate categories’ (‘CAT’) model (Stamatakis 2006); codon-based alignment (i.e. where the output is nucleotide sequences, but the alignment is carried out using amino acid sequences); removal of alignment positions composed of more than 50% gap characters; 1000 resamples during tree construction phase. The UBCG pipeline was run with Prodigal for gene finding (ver. 2.6.3, see https://github.com/hyattpd/prodigal/releases/, accessed September 2018; Hyatt et al. 2010), hmmsearch for identification of UBCGs using HMMER (ver. 3.1b2, see http://eddylab.org/software/hmmer/hmmer-3.1b2.tar.gz, accessed September 2018; Eddy 2011), Multiple Alignment using Fast Fourier Transform (MAFFT) for alignments (ver. 7.402, see https://mafft.cbrc.jp/alignment/software/linux.html, accessed September 2018; Katoh and Standley 2013) and FastTree 2 (ver. 2.1.10, Double precision, number SSE3, see http://www.microbesonline.org/fasttree/, accessed September 2018; Price et al. 2010) to construct single gene trees and an unrooted maximum-likelihood (ML) tree from the concatenated alignment of the nucleotide sequences of 91 genes representing the bacterial core-gene set (length = 84 294 nucleotide positions; see Fig. S1, available as Supplementary material to this paper). The ML tree was subsequently rooted by applying the Minimum Ancestor Deviation function (ver. 1.1, https://www.mikrobio.uni-kiel.de/de/ag-dagan/ressourcen, accessed September 2018; Tria et al. 2017) with default settings, then visualised with the APE package (ver. 5.0, https://cran.r-project.org/web/packages/ape/, accessed September 2018; Paradis et al. 2004).

To characterise the phylogenetic distribution of SM gene clusters, Phylo.D (caper ver. 1.0.1, D. Orme, R. Freckleton, G. Thomas, T. Petzoldt, S. Fritz, N. Isaac, and W. Pearse, see https://CRAN.R-project.org/package=caper, accessed September 2018) was used to estimate the extent of phylogenetic clumping. Phylo.D (run with 10 000 permutations per trait model) calculates the D value (Fritz and Purvis 2010), a measure of phylogenetic signal in a binary trait, and tests the estimated D value for significant departure from both random association and the clumping expected under a Brownian evolution threshold model.

Putative genomic islands within the 27 M. aeruginosa genomes were identified using AlienHunter (ver. 1.7, see https://www.sanger.ac.uk/science/tools/alien-hunter, accessed September 2018) with default settings, which predicts putative horizontal gene transfer (HGT) events with the implementation of interpolated variable order motifs (IVOMs; Vernikos and Parkhill 2006). The IVOM approach exploits compositional biases at various levels (e.g. codon, dinucleotide and amino acid bias, structural constraints) and more reliably captures the local composition of a sequence than fixed-order methods. The relative position of SM biosynthesis clusters compared to IVOMs within the four closed M. aeruginosa genomes (NIES-2481, NIES-2549, NIES-843 and PCC 7806SL) was also examined.

The genomes ranged in size from 4.3 to 5.9 Mb (mean = ~5 Mb) with GC contents ranging from 42.1 to 43.2% (mean = ~43%). They were assembled across 1 to 1363 scaffolds (mean = ~235) and contained 3839–5190 protein-coding regions (mean = ~4417; Table 1). In pairwise genome comparisons (Fig. S1), 85.2% of ANI values (which ranged from 94.3 to 99.9%; median 95.3%) were above the nominal prokaryote species threshold of 95% (Goris et al. 2007; Fig. 3).

Table 1.  Microcystis aeruginosa strains used in this study and their corresponding genome features
Aw, tropical wet; Cfa, humid subtropical; Cfb, temperate oceanic; Cwa, monsoon-influenced humid subtropical; Cwb, subtropical highland; Dfb, warm-summer humid continental; GC%, guanine–cytosine content; INSDC, International Nucleotide Sequence Database Collaboration accession number; IVOM, interpolated variable order motif; Level, genome completion level; K-G classification, Köppen–Geiger climate zone; Location, isolation site

Fig. 3.  Orthologous average nucleotide identity (OrthoANI) boxplot. The dashed horizontal line marks the nominal prokaryote species threshold of 95% (Goris et al. 2007). The box shows the interquartile range, with the median value indicated by the horizontal line; whiskers show the range. Individual symbols are outliers.

Homologues of nine known SM biosynthesis gene clusters were distributed among the 27 M. aeruginosa genomes, including three PK-NRP biosynthesis clusters (aeruginosin, aer; microcystin, mcy; and microginin, mic), two NRP biosynthesis clusters (anabaenopeptin, apn; and micropeptin, mcn), three cyanobactin biosynthesis clusters (aeruginosamide, age; microcyclamide, mca; and piricyclamide, pir) and a microviridin (mdn) biosynthesis cluster. The individual M. aeruginosa strains encoded 0–7 known SM biosynthesis clusters (mean = 4), which occupied ~0–3% of their total genomes (mean ≥2%). PCC 9432 encoded the greatest number of known clusters, whereas NIES-44 encoded none, although fragments of the pir cluster were found in the NIES-44 genome. The most commonly identified SM biosynthesis clusters were the aer and mdn clusters, which were identified in 21 different strains. In contrast, the age cluster was only found in a single strain (PCC 9432). The most common class of SM biosynthesis cluster identified among the 27 strains was the mixed PK-NRP class, which was identified at approximately twice the frequency of the NRP, cyanobactin and microviridin classes (Tables 2, S1).

Table 2.  Distribution of known secondary metabolite biosynthesis gene clusters in 27 Microcystis aeruginosa strains
+, gene cluster identified with high homology (>80% ID) to Minimum Information about a Biosynthetic Gene cluster (MIBiG) reference gene clusters: aer, aeruginosin; mcy, microcystin; mic, microginin; apn, anabaenopeptin; mcn, micropeptin or cyanopeptolin; age, aeruginosamide; mca, microcyclamide; pir, piricyclamide; and mdn, microviridin B. P, partial microcystin cluster (mcyA–C); NRPS, non-ribosomal peptide synthetase; PKS, polyketide synthase

The complete microcystin (mcyA–J) biosynthesis gene cluster was identified in 15 M. aeruginosa strains, including the reference strain PCC 7806SL (MIBiG accession BGC000107). A partial microcystin synthetase cluster (mcyA–C) was identified in PCC 9717. Genes within the mcy cluster were highly conserved with inferred peptide sequence IDs ≥91% (mean = ~98%) compared with their homologous counterparts in strain PCC 7806SL.

The aeruginosin (aer) biosynthesis gene cluster was identified in 21 strains, including the reference strain NIES-98 (MIBiG accession BGC0000298). The consensus cluster within this group comprised aerA, aerB, aerD, aerE, aerF, aerG1, aerK, aerL, aerN and OxRed genes. A putative transposase gene was also identified just upstream of aerA. The full complement of aer genes (aerA, aerB, aerD, aerE, aerF, aerG1, aerG2, aerK, aerL, aerM, aerN) was only detected in three strains (NIES-98, PCC 7005 and TAIHU98). An additional open reading frame (ORF), namely aerJ, putatively encoding a halogenase, was detected in seven strains. The inferred peptide sequences of aer genes were mostly well conserved (mean = ~93% ID) compared with the NIES-98 sequences. However, the sequences of aerB and aerL were poorly conserved in some strains (≥58%; mean = ~76% ID).

The closest match to the microginin (mic) gene cluster in the MIBiG repository was the puwainaphycin cluster from Cylindrospermum alatosporum CCALA 988 (accession BGC0001125). Therefore, 27 M. aeruginosa genomes were screened for the mic cluster by BLASTP (Altschul et al. 1990) analysis with MicA, MicC, MicD and MicE from PCC 7941 used as query sequences (GenBank accession CCI09456–9). Homologues of all mic genes were identified in seven additional strains with an average inferred peptide sequence identity of ~90% compared with the reference sequences. Among these sequences, the PKS-NRPS micA and ABC transporter micE were the most conserved (≥90%; mean = ~95% ID), whereas the NRPSs micC and micD were less conserved (≥64%; mean = ~84% ID).

The anabaenopeptin (apn) biosynthesis gene cluster was identified in three strains (PCC 7005, PCC 9432 and PCC 9701). All identified clusters encoded the full complement of apn genes (apnA, apnB, apnC, apnD, apnE). Sequence conservation was generally high compared with the Planktothrix agardhii NIVA-CYA 126/8 reference sequences (MIBiG accession BGC0000301; mean = ~88% ID). The NRPS genes apnA and apnB and the ABC-transporter gene apnE were particularly well conserved in terms of inferred peptide sequence (≥92%; mean = ~93% ID), whereas the NRPS genes apnC and apnD were less conserved (≥76%; mean = ~80% ID). As seen in P. agardhii, but in contrast with Anabaena sp. 90 (Rouhiainen et al. 2010), the M. aeruginosa apn gene clusters only had one copy of apnA.

The cyanopeptolin (mcn) biosynthesis gene cluster was identified in 19 strains. Only three of these clusters (NaRes975, PCC 7941 and PCC 9808) encoded the full complement of mcn genes (mcnA–G). However, most of the mcn clusters identified (17/19) encoded homologues of mcnB, mcnC, mcnE, mcnF and mcnG. Approximately half the mcn clusters identified (10/19) encoded the NRPS mcnA, whereas only four clusters encoded the halogenase mcnD. Interestingly, four of the strains that lacked mcnA (CHAOHU 1326, NIES-2481, NIES-2549 and NIES-298) had elongated (~2× compared with the reference) mcnB genes encoding an additional peptidyl carrier protein domain and an unusual FkbH domain. The inferred peptide sequences of these elongated mcnB homologues were poorly conserved (mean = ~65% ID) compared with the Microcystis sp. NIVA-CYA 172/5 reference sequence (MIBiG accession BGC0000332), but highly similar to each other. Sequence conservation among the remaining mcn genes was generally high (mean = ~92% ID).

The microcyclamide (mca) biosynthesis gene cluster was identified in eight strains, including the reference strain PCC 7806SL (MIBiG accession BGC0000474). All clusters encoded the full complement of mca genes (mcaA–G), except TAIHU98, which lacked the subtilisin-like protein mcaA. Overall inferred peptide sequence conservation was high compared with the reference cluster (mean = ~95% ID). However, the sequence of the microcyclamide precursor gene mcaE was only moderately conserved (mean = ~82% ID).

An aeruginosamide (ageA, ageB, ageD, ageF1, ageF2, ageG) biosynthesis gene cluster was identified in only one strain, the reference strain PCC 9432 (MIBiG accession BGC0000483). The aer cluster has similar organisation and identity to the microcyclamide (mca) cluster, except it also encodes a methyltransferase gene ageF1.

Piricyclamide-like biosynthesis gene clusters were identified in 12 strains, including the reference strain PCC 7005 (MIBiG accession BGC0001167). Only one other strain (PCC 9807) encoded the full complement of pir genes (pirA, pirB, pirC, pirE2, pirE3, pirE4, pirF, pirG). All piricyclamide-like biosynthesis clusters contained homologues of the N- and C-terminal cyanobactin protease genes pirA and pirG and the hypothetical protein pirB. Most contained the hypothetical protein pirC, the prenyl transferase gene pirF and at least one cyanobactin precursor (pirE2, pirE3, pirE4). Overall inferred peptide sequence identity averaged ~86% compared with the reference cluster from PCC 7005. Sequence conservation was high among pirA, pirB, pirG and pirE4 genes (mean = ~94% ID), but lower among pirC, pirF, pirE2 and pirE3 genes (mean = ~75% ID).

The microviridin (mdn) biosynthesis gene cluster was identified in 21 strains, including the reference strain NIES-298 (MIBiG accession BGC0000592). All strains encoded mdnB, mdnC, mdnD and mdnE, except PCC 7806SL, which lacked the acetyltransferase gene mdnD. Seven strains encoded the microviridin precursor gene mdnA. The remaining 14 strains encoded a gene with low inferred peptide sequence homology (mean = ~51%) to mdnA. This gene was located outside the mdn cluster in nine strains. The sequences of mdnB–E genes were highly conserved in all clusters (mean = ~97% ID). Fragments of the mdn cluster were also identified in PCC 9806 and TAIHU98.

Maximum-likelihood phylogenetic analysis of the concatenated alignment of 91 bacterial core genes grouped the 27 M. aeruginosa genomes into 25 clades exhibiting a branching pattern broadly concordant with that found in the phylogenetic analysis of 10 strains by (Humbert et al. 2013), with generally high local support values (Fig. 4). The number of single gene trees supporting a given branch in the UBCG tree (designated the gene support index, GSI) ranged from 11 to 91, with lower GSI values observed in early branch splits. The 16S rRNA locus was also examined as a candidate single-locus phylogenetic marker, but because there was only a 0.5% sequence variation across aligned positions, this locus was not able to resolve closely related M. aeruginosa strains (data not shown).

Fig. 4.  Rooted maximum likelihood up-to-date bacterial core gene (UBCG) tree inferred using the concatenated alignment of 91 bacterial core genes. A total of 84 294 nucleotide positions was used. The number of single gene trees supporting a branch was calculated and designated the gene support index (GSI; Na et al. 2018). GSI and percentage local support values are given at branching points as GSI/percentage local support value. Köppen–Geiger climate zones and continents of origin are given in bold: Cfa, humid subtropical; Dfb, warm-summer humid continental; Cfb, temperate oceanic; Cwa, monsoon-influenced humid subtropical; Aw, tropical wet; Cwb, subtropical highland. AS, Asia; EU, Europe; NA, North America; SA, South and Central America; AF, Africa; OC, Oceania. SC1–SC4, marked with vertical bars, delineate the major phylogenetic subclades.

Representative biosynthesis gene clusters from the four SM classes were detected in each of the major phylogenetic subclades (SC1–4). However, the distribution of individual SM biosynthesis clusters varied across the subclades. For example, all strains from SC1 lacked mcy, age and mca clusters; all strains from SC2 possessed the aer, mcy and mcn clusters, but lacked the mic, apn and age clusters; all strains from SC3 possessed mcy and mdn clusters, but lacked apn, age and mca clusters; and all strains from SC4 possessed the aer cluster (Fig. 4; Table 2).

Phylo.D analysis revealed that the phylogenetic distribution of all PK-NRP biosynthesis clusters, particularly mcy, was clumped, the distribution of the apn cluster was random and the distribution of the age cluster (n = 1) was overdispersed. The phylogenetic patterning of the other gene clusters was less marked (Tables 3, 4).

Table 3.  Phylogenetic distribution of secondary metabolite (SM) biosynthesis gene clusters in Microcystis aeruginosa strains (with DIANCHI905 excluded)
Data show D values (ordered from lowest to highest) and associated probabilities for random and Brownian phylogenetic patterns (Fritz and Purvis 2010), constrained to avoid tips with zero branch lengths by excluding DIANCHI905 from the analysis. The dummy variables ‘clumped’ and ‘dispersed’ were added to the trait matrix as interpretative aids and 10 000 permutations were run per trait model. aer, aeruginosin; age, aeruginosamide; apn, anabaenopeptin; mca, microcyclamide; mcn, micropeptin or cyanopeptolin; mcy, microcystin; mdn, microviridin B; mic, microginin; NA, not applicable; NRPS, non-ribosomal peptide synthetase; pir, piricyclamide; PKS, polyketide synthase

Table 4.  Phylogenetic distribution of secondary metabolite (SM) biosynthesis gene clusters in Microcystis aeruginosa strains (with PCC 7806SL excluded)
Data show D values (ordered from lowest to highest) and associated probabilities for random and Brownian phylogenetic patterns (Fritz and Purvis 2010), constrained to avoid tips with zero branch lengths by excluding PCC 7806SL from the analysis. The dummy variables ‘clumped’ and ‘dispersed’ were added to the trait matrix as interpretative aids and 10 000 permutations were run per trait model. aer, aeruginosin; age, aeruginosamide; apn, anabaenopeptin; mca, microcyclamide; mcn, micropeptin or cyanopeptolin; mcy, microcystin; mdn, microviridin B; mic, microginin; NA, not applicable; NRPS, non-ribosomal peptide synthetase; pir, piricyclamide; PKS, polyketide synthase

The 27 M. aeruginosa strains, from six different continents, were allocated to six different Köppen–Geiger climates: 12 to Cfa (humid subtropical), 7 to Dfb (warm-summer humid continental), 4 to Cfb (temperate oceanic), 2 to Cwa (monsoon-influenced humid subtropical), 1 to Aw (tropical wet) and 1 to Cwb (subtropical highland).

Representative gene clusters from the four SM classes were detected in each of the six continental groups, except the South American group (n = 1), which lacked the NRP class. However, the distribution of individual SM biosynthesis clusters varied across the geographical groups. For example, the African strains (n = 2) both possessed aer, mcy, mcn, pir and mdn gene clusters, but lacked mic, apn, age and mca clusters. The Asian strains (n = 15) all lacked apn and age clusters. The Australian strain (n = 1) possessed aer, mcy, mcn, pir and mdn gene clusters, but lacked mic, apn and age clusters. The European strains (n = 3) all possessed mcn and mdn clusters, but lacked the age cluster. When Köppen–Geiger climate zones (with n ≥ 2) were considered, all strains within the Cfa group (n = 12) lacked apn and age clusters; all strains within the Cfb group (n = 5) lacked the age cluster, but possessed the mcn and mdn clusters; and all strains within the Cwa group (n = 2) possessed the aer, mcy, mcn and mdn clusters, but lacked the mic, apn, age and mca clusters (Fig. 4; Table 2).

The number of IVOMs (Vernikos and Parkhill 2006) in the 27 M. aeruginosa genomes ranged from 127 to 246 (mean = 190; Table 1). A positive correlation between the number of IVOMs and genome size was observed (Fig. 5), but there was no correlation between the number of IVOMs and the number of SM biosynthesis clusters identified (data not shown). Examination of the relative position of SM biosynthesis gene clusters in the four closed M. aeruginosa genomes (NIES-2481, NIES-2549, NIES-843 and PCC 7806SL) revealed that the mcy, aer, mcn, mca and pir gene clusters completely or partially overlap with IVOMs (37–100%; mean = 84.2%). In contrast, the mdn cluster did not overlap with IVOMs (Table 5).

Fig. 5.  Correlation between interpolated variable order motifs and genome size. IVOM, interpolated variable order motif; Pr(>|t|), P-value for the t-test.

Table 5.  Relative position of secondary metabolite (SM) biosynthesis gene clusters and interpolated variable order motifs (IVOMs) in closed Microcystis aeruginosa genomes
Coordinates for SM cluster and IVOM start and stop positions are given. Overlap, percentage of cluster that overlaps with IVOM; aer, aeruginosin; mca, microcyclamide; mcn, micropeptin or cyanopeptolin; mcy, microcystin; mdn, microviridin B; pir, piricyclamide

Automated genome mining pipelines are powerful tools for identifying SM pathways in microorganisms, particularly when biochemical and genetic data are available for similar pathways in related organisms (Blin et al. 2017). In the present study, antiSMASH and BLAST analyses were used to screen 27 M. aeruginosa genomes for 1481 known SM pathways. Nine known SM biosynthesis gene clusters from four different metabolite classes were identified. However, the only cyanotoxin gene cluster identified within the 27 M. aeruginosa genomes was the mcy cluster, suggesting that the production of non-microcystin cyanotoxins by this taxon is absent or rare. This finding is in contrast with previous reports that the M. aeruginosa is also capable of producing anatoxin-a (strains TAC80, TAC87, TAC117 and TAC121; Park et al. 1993) and paralytic shellfish poison analogues (Sant’Anna et al. 2011).

Although the number of known SM pathways detected in each genome varied, the fact that SM biosynthesis gene clusters comprised on average 2% of the total genome size suggests that the corresponding compounds play an important ecological role in M. aeruginosa. The high frequency of aeruginosin, microviridin, cyanopeptolin and microcystin biosynthesis gene clusters further suggests that these compounds are particularly significant to the success of the species. However, the lack or complete absence of SM biosynthesis clusters in some strains (e.g. NIES-44) seems to contradict this. It is likely that the different SM cluster profiles observed in this study correspond to different M. aeruginosa ecotypes that have evolved to suit specific niche environments. Previous growth competition studies on toxic v. non-toxic M. aeruginosa strains support this theory, demonstrating that microcystin production, although not essential for survival, may be advantageous under certain growth conditions, such as nutrient limitation or high light (Kaebernick et al. 2000; Kardinaal et al. 2007; Zilliges et al. 2011).

Although gene sequence identity was generally high between homologous counterparts in each SM gene cluster, the composition of genes in each cluster sometimes varied. Gene identity and composition were particularly well conserved among microcystin, microginin, anabaenopeptin and microcyclamide pathways, suggesting that all mcy, mic, apn and mca genes play key roles in the biosynthesis or activity of their corresponding compounds, and major mutations are not tolerated. Previous studies have shown this to be the case. For example, mcyB- and mcyH-knockout mutants were unable to produce microcystin (Dittmann et al. 1997; Pearson et al. 2004), whereas mcyJ-knockout mutants produced toxin lacking the O-methylation on the C₉ hydroxyl unit and had reduced inhibitory activity against protein phosphatases (Christiansen et al. 2003). Based on these examples, it is anticipated that strain PCC 9717, lacking mcyD–H, is unable to synthesise microcystin.

Interestingly, although microviridin biosynthesis genes were well conserved, homologues of the precursor peptide gene mdnA were poorly conserved in terms of sequence conservation and located outside the mdn gene cluster in some strains. The distal location of mdnA, although recognised previously (Ziemert et al. 2010; Humbert et al. 2013), is unusual given that most of SM biosynthesis genes are arranged in operons. Sequence variations in mdnA have recently been shown to underpin chemical diversity in the microviridin family (Gatte-Picchi et al. 2014). Furthermore, the two grasp ligases encoded by mdnB and mdnC could potentially modify a variety of other leader peptides scattered throughout the genome, as has been observed for the prochlorosin synthetase of Prochlorococcus MIT9313, which processes up to 29 different precursor peptides (Zhang et al. 2014). The absence of the acetyltransferase gene mdnD in PCC 7806SL suggests that novel microviridin variants lacking the N-acetylation could also be produced by this strain. Together, these results suggest that several different microviridins are produced by the M. aeruginosa group analysed in this study.

The identity and composition of core biosynthesis genes in the aeruginosin biosynthesis cluster was well conserved, but many strains lacked the halogenase gene aerJ and the NRPS genes aerG2 and aerM. The loss of aerJ suggests that these strains may produce non-chlorinated variants of aeruginosin. A previous study by Cadel-Six et al. (2008) has shown this to be the case for numerous M. aeruginosa strains, including several analysed in this study. The absence of the additional NRPS genes further suggests that some strains may produce aeruginosin variants with different amino acid backbones.

The composition of genes in the cyanopeptolin and piricyclamide biosynthesis clusters was highly variable, and no consensus gene sets could be established for these pathways. The loss of different genes in different strains suggests that a wide variety of cyanopeptolin- and piricyclamide-like compounds could be produced by M. aeruginosa. For example, strains lacking mcnD are likely to produce non-chlorinated cyanopeptolins, whereas strains lacking mcnD and mcnG are likely to produce micropeptins. The fusion of an FkbH domain to mcnB was observed in several cyanopeptolin biosynthesis gene clusters. This is interesting because most of the reported FkbH-like proteins are involved in polyketide biosynthesis, with cyanopetolin-1138 biosynthetic machinery as the only example of an FkbH-like protein found in a pure NRPS setting (Auerbach et al. 2018). The mosaic structure of these mcnB genes suggests that they may be products of horizontal gene transfer (HGT; Rounge et al. 2007).

Compositional and sequence diversity among the piricyclamide biosynthesis clusters was not surprising considering the size and complex evolutionary history of the cyanobactin family (Leikoski et al. 2013). However, the loss of multiple pir genes in some M. aeruginosa strains could be a sign of inactive pathways. For example, the absence of all but two piricyclamide genes, namely pirE3 and pirG, in the NIES-44 genome suggests that this strain is unable to produce cyanobactins.

Highly conserved genes within compositionally conserved SM biosynthesis clusters are good candidate targets for molecular diagnostic tests. PCR-based tests targeting core genes within the microcystin, nodularin, cylindrospermopsin and saxitoxin pathways have already proven to be reliable methods for detecting toxic cyanobacteria (Al-Tebrineh et al. 2010, 2011; Baker et al. 2013). Although mass spectrometry methods are still considered the ‘gold standard’ for detecting cyanotoxins, molecular methods are becoming more widely accepted and are often preferred by water quality managers and researchers because they are quicker, more economical and user friendly and provide evidence for toxigenic potential. Similar tests could be designed for the diagnosis of the other potentially harmful SM biosynthesis gene clusters identified in this study, but the choice of target genes may be limited in the case of the less conserved clusters, such as mcn and pir.

In addition to investigating the composition and conservation of SM biosynthesis gene clusters in M. aeruginosa, this study sought to determine whether phylogeny and geography or climate zone were good indicators of SM profiles. Although phylogeny appeared to be correlated with the distribution of certain SM biosynthesis clusters within certain subclades (particularly PK-NRP biosynthesis clusters), these relationships are not strong enough to allow reliable prediction of SM profiles based on phylogenetic markers alone. Similarly, although geoclimatic forces seem to affect SM gene cluster profiles in this group, further studies with larger sample sizes are required to establish the significance of these results.

The complex distribution of SM biosynthesis clusters in cyanobacteria has been mostly attributed to gene loss events, rather than HGT. For example, phylogenetic studies on mcy, aer and mcn gene clusters suggest that these pathways existed in an ancient common ancestor and were lost in subsequent lineages (Rantala et al. 2004; Rounge et al. 2007). Although the lateral transfer of whole clusters seems unlikely, HGT could explain some of the sequence variations observed in earlier studies and herein.

The 27 M. aeruginosa genomes analysed in this study contained a large number of IVOMs, suggesting that frequent HGT events occur within this species. Examination of the relative position of SM biosynthesis clusters in the four closed genomes revealed that the aer, mcy, mcn, mca and pir clusters overlap partially or completely with IVOMs. This suggests that HGT events may have led to the acquisition or recombination of genes within these clusters. These results are in agreement with previous phylogenetic studies that suggest HGT events led to DNA polymorphisms in the mcy and aer and cyanobactin biosynthesis clusters (Mikalsen et al. 2003; Tanabe et al. 2004; Ishida et al. 2009; Leikoski et al. 2009). The occurrence of transposase genes within, and proximal to, SM biosynthesis clusters, including aer (present study), mcy (Tillett et al. 2000), mcn (Nishizawa et al. 2011), mca (Ziemert et al. 2008b) and pir (Leikoski et al. 2013), lends further support to this theory.

Recombination events in cyanobacterial genomes can give rise to diverse SM variants. Knowing which variants are produced by a given strain can be important for both drug discovery and water management purposes because different variants often have different activities and toxicities. For example, non-sulfated variants of saxitoxin are ~10-fold more toxic than disulfated variants (Wiese et al. 2010). Similarly, the toxicity of microcystin-RR is approximately 10-fold higher than that of microcystin-LR (Sivonen and Jones 1999). Whole-genome approaches, capable of predicting SM profiles by the analysis of conserved domains within biosynthesis and tailoring enzymes, could foreseeably become routine features of water quality monitoring as rapid sequencing technologies become more accessible.

Cyanobacterial SMs exhibit a wide range of bioactivities, including the inhibition of eukaryotic protein phosphatases. These activities have the potential to disrupt aquatic ecosystems and affect water quality. This study has revealed the distribution, composition and conservation of 9 different SM biosynthesis gene clusters in 27 M. aeruginosa genomes, highlighting the potential chemical diversity inherent within this species. Furthermore, the combined results suggest that multiple factors, including geography and climate, gene loss and HGT, have played a role in shaping the evolution of these pathways. This research has also laid the foundations for future molecular screening tests for predicting SM profiles in uncharacterised, but potentially harmful, Microcystis cultures and bloom samples.

The authors declare that they have no conflicts of interest.

Funding was received from Melbourne Water Corporation as part of the Nuisance and Harmful Algae Science–Practice Partnership and from the Australian Research Council (Grant LP130100311).