Our genomic search for genes encoding H₂ase and the proteins involved in their maturation helped to complete the picture of which strains may possibly synthesize functional H₂ase.

A phylum-wide analysis of the genomic distribution of H₂ase genes among the cyanobacterial genomes in the CyanoGEBA dataset (Shih et al., 2013) showed that only [NiFe] H₂ases are present in these organisms. No obvious homologs of [FeFe] or [Fe] H₂ases were identified. We assumed that only genomes possessing all the hox and hup genes carry a complete set of H₂ase-encoding genes. A complete set of H₂ase-encoding genes was deciphered in 52% of the genomes studied (Figure 1A). Among the 130 genomes analyzed, 49 did not show any H₂ase-encoding genes (Figure 1A), and 13 genomes did not present the complete set of genes required to encode a functional H₂ase (Supplementary Table 1). The lack of H₂ase genes may be attributable to the bacterial habitat, since the proportions of H₂ases-free genomes differ from one ecological niche to another: the highest proportion of H₂ase-free strains was detected in the open ocean (89%), and the remaining 11% carried only hup genes, which suggests that the cyanobacterial contribution to H₂ production in the open ocean is negligible (Figure 1B). The distribution of H₂ase genes and of genes required for their maturation was found to vary in the cyanobacterial phylum, but all the organisms belonging to subsections IV and V have a complete set of genes encoding H₂ases. H₂ oxidation and H⁺ reduction activities seem to be generally conserved in these species (Figure 1C). Since the uptake H₂ase is involved in functional nitrogenase processes, the co-occurrence of H₂ases, and nitrogenase in various environments was investigated by studying the distribution of [FeMo] nitrogenase structural genes (nifH and nifDK), the nifBEN, and the nifUSV genes involved in the synthesis of the [FeMo]-cofactor synthesis. The nifH, nifDK, and nifBEN genes were found in all the cyanobacteria genomes (Supplementary Table 2). The nifBEN genes were found in co-occurrence with nifUSV genes except is six genomes (Supplementary Table 2). Since the nifSU genes have been reported to be dispensable in Anabaena variabilis (Lyons and Thiel, 1995), one might conclude that their absence does not necessarily mean that the strain is not able to fix nitrogen. It is therefore concluded that all the strains listed in Supplementary Table 2, and whose genomes contain nifH, nifDK, and nifBEN genes are potentially nitrogen-fixing.

The data obtained here, indicate that nif genes are present in genomes harboring hup (10%) or hox (8%), or both (15%). Eleven percent of the genomes possess nifH and nifDK without harboring the hox and hup genes (Figure 1A). The co-occurrence of nif and hup genes seems to be significantly more frequent in the genomes of strains belonging to subsection V (Figure 1C).

To further assess the distribution of genes encoding H₂ases among the cyanobacterial phylum, a phylogenetic tree was constructed using 21 concatenated sequences corresponding to the 130 cyanobacterial proteins listed in Supplementary Table 3 (see Methods Section). The hox and hup genes were found to occur more frequently in the late branches of the tree, although the distribution of hox is patchier (Figures 2, 3). The presence of hyp genes was always associated with that of at least one of the Hox or Hup sequences, and these genes therefore occur less frequently in the early branching clades, which is clearly illustrated in the case of clade g (Figures 2, 3). The distribution of the nitrogenase-encoding genes (nifH and nifDK) is in agreement with the phylogenetic tree previuously presented (Bandyopadhyay et al., 2010). These genes are present in four genomes in the early branches of cyanobacterial evolution: clade a [Synechococcus sp. JA-2-3B′a(2–13), Synechococcus sp. JA-3-3Ab], clade b (Pseudanabaena sp. PCC 6802), clade c (Cyanothece sp. PCC 7425), and six genomes of clade d. No hup genes were ever detected in these early clades, which suggests that the nitrogenase may function naturally in the absence of uptake H₂ase. This was previously found to occur in Synechococcus sp. and Cyanothece PCC 7425, which fix N₂ under anaerobic conditions (Bandyopadhyay et al., 2011). In the genomes of the strains Nostocales and Stignematales (subsections V and VI), which belong to clade h, the nif and hup genes were always found to co-occur (Figures 2, 3). It is also worth noting that the co-occurrence of hox, hup and nif genes was observed only in the late branches of the tree.

The genes encoding the bidirectional H₂ases (hoxY, hoxH) and the hoxU, hoxE, and hoxF genes encoding the diaphorase part are widely distributed among the cyanobacterial phylum and are particularly abundant in the genomes of organisms belonging to subsections II, III, and IV (Figures 1B,C and Supplementary Table 4). All the hox genes listed in Supplementary Table 4 potentially encode soluble H₂ases belonging to subgroup 3d (Vignais et al., 2001). In the large subunit (HoxH), the sequences of the L1 and L2 motifs typical of each [NiFe] group show a high level of conservation. Only a few amino-acid substitutions were observed in the L1 motif in seven genomes of strains from various habitats (terrestrial, coastal, and freshwater strains) (Supplementary Table 5). The Cysteine residues involved in the coordination of metal ions are strictly conserved in all the HoxH and HoxY sequences. The three subunits in the diaphorase of the bidirectional H₂ase part (HoxE, HoxF, and HoxU) also contain the conserved cysteine residues potentially required for the coordination of [2Fe2S] and [4Fe4S] clusters. These cysteine residues are largely conserved, since the only few exceptions observed were HoxF and HoxU proteins in Synechococcus sp. CB0205, P. hollandica PCC 9006, and Cyanobium sp. PCC 7001 (Supplementary Table 1). These genomes also lack some of the genes involved in the maturation process (Supplementary Table 1). The bidirectional H₂ase in these strains may therefore not be active. The last step in the maturation of the bidirectional H₂ases involves the HoxW endopeptidase. The co-occurrence of the hoxW gene and the hox structural genes (HYUEF) was observed in all the genomes analyzed (Supplementary Table 4). Based on the difference between the patterns of expression of the structural hox genes and hoxW, it has been suggested that the endopeptidase HoxW might have multiple functions in cyanobacteria (Wünschiers et al., 2003). The results of the present study confirm this assumption, since hoxW homologs were found to exist in four genomes containing no hoxYHUEF genes (Supplementary Table 4). In addition, the presence of multicopies of the hoxW gene observed in three genomes provides a further argument supporting this hypothesis (Supplementary Table 4).

Seven different patterns of organization were observed among the structural hox genes (Figure 4A, Supplementary Figure 1). In Group 1, the hoxE, hoxF, hoxU, hoxY, and hoxH genes are clustered together and show the same orientation, whereas the hoxW gene occupies another position in the genome. Group 1 includes 26 genomes belonging to all the subsections except subsection V. Group 2 includes 20 genomes belonging to subsections I, II and IV, and all the structural hox genes (EFUYHW) are clustered together in the same orientation (Figure 4A). Group 3 contains two genomes belonging to subsections I and V: the hoxE, hoxF, hoxU, hoxY are clustered together and in the same orientation, whereas the hoxW and hoxH are located in another part of the genomes. The hox genes are more widely scattered in Groups 4-6: hoxE and hoxF are clustered together and the other hox genes are either clustered or scattered in various combinations. All the hox genes hoxF, hoxU, hoxY, hoxH, and hox W are clustered together in Fiscarella sp. PCC 9605 (Group 7), whereas hoxE is located in another part of the genome. The organization of the hox genes is generally not conserved throughout the tree of species, where the seven groups are randomly distributed among the eight clades (Figures 2, 3).

HupL and HupS homologs encoding the large and small subunits of the uptake H₂ase, respectively, were identified only in genomes of diazotrophic strains belonging to subsections I, III, IV, and V (Figure 1C and Supplementary Table 6). The strains carrying uptake H₂ase genes are widely distributed in various habitats. They are absent only in the genomes of strains collected from salt lakes (Figure 1B). The amino acid sequences of HupS and HupL show a high degree of conservation: the L1 and L2 motifs typical of H₂ases belonging to group 2a (Vignais and Billoud, 2007) were found to be conserved in all the Hup sequences analyzed. These motifs include the cysteine residues involved in the coordination of the [NiFe] in the case of HupL and [FeS] in that of HupS (Vignais and Billoud, 2007). In the genomes of Calothrix sp. PCC 7103, Tolypothrix sp. PCC 9009, Rivularia sp. PCC 7116, Cyanobacterium sp. UCYN A2, and Calothrix sp. PCC 6303, since the HupS sequence shows deletions and substitutions of the residues involved in the binding of [FeS] cluster, these enzymes may be inactive or might have different enzymatic characteristics (Supplementary Figure 2). The specific peptidase HupW was identified in all the genomes carrying hupSL genes (Supplementary Table 6). The HupW sequences consistently showed well-conserved residues thought to contribute importantly to the specific interactions between the peptidase and its cognate H₂ase subunit (Devine et al., 2009).

In all the genomes analyzed, the hupS and hupL genes form clusters. The organization of the five groups of hup genes depends on the location of the hupW gene and the disruption (or otherwise) of hupS or hupL genes by the xisC gene (Figure 4B). The distribution of these clustering groups varies in the tree of species (Figures 2, 3). Groups 1 or 2 are mostly present throughout the late branches of the tree (clades f, e, g, and h), whereas groups 3, 4, and 5 occur only in clade h (Figure 4B).

Almost all the cyanobacterial genomes harboring structural H₂ase genes (hox, hup, or both) also harbor the hypABCDEF genes known to encode proteins involved in the maturation of the H₂ase (Supplementary Table 7), apart from the genomes of Chroococcidiopsis thermalis PCC 7203, Synechococcus elongatus PCC 7942, Synechococcus sp. CB0101, and Synechococcus sp. PCC 7336, from which some hyp genes are missing. (Supplementary Table 1). Whether the maturation of the H₂ase in these strains involves different mechanisms, or whether the maturation process is not efficient in these case is still an open question.

Since little is known about the process of H₂ase maturation in cyanobacteria, we analyzed the amino acid composition of the Hyp proteins in the light of the data available in the literature on other organisms. All the information based on the resolution of the crystal structure of the HypF protein of Caldanaerobacter subterraneus (Shomura and Higuchi, 2012), that of the HypECDA of Thermococcus kodakarensis (Watanabe et al., 2009, 2012; Tominaga et al., 2013) and that of the HypB of Archaeoglobus fulgidus, Bradyrhizobium japonicum, and Escherichia coli (Olson and Maier, 2000; Chan et al., 2012; Douglas et al., 2013) are summarized in Supplementary Table 7 and Supplementary Figures 3–8. The fact that the cyanobacterial Hyp sequences showed highly conserved residues reported to contribute importantly to the Hyp features (Supplementary Table 7 and Supplementary Figures 3–8) suggests that the process of maturation of the H₂ase enzymes in cyanobacteria might be similar to that described in other organisms (Hansel et al., 2001; Shomura and Higuchi, 2012; Watanabe et al., 2012; Douglas et al., 2013; Tominaga et al., 2013). The hyp genes are either clustered together in various combinations or scattered throughout the genome without any correlations being detected with the diazotrophic ability of the strains or their habitat or their classification (Figure 4C, Supplementary Figure 1). The hyp genes can be classified into 11 main classes depending on their patterns of organization. The genomes in class 1 carry all the hyp genes in a single cluster, while those in class 2 carry five clustered hyp genes and one gene located in another part of the genome, for example. Many rearrangements of the hyp clusters have occurred during the evolution of cyanobacteria, and the number of clusters increases in the late branches of the tree. In clade h, the genes are all clustered together and show a similar pattern of organization (Figures 2, 3, Supplementary Figure 1).

A search for homologs of O₂-tolerant H₂ases encoding genes in all the cyanobacterial genomes available in the NCBI database yielded positive findings in five genomes (Table 2, Supplementary Table 8). A blast analysis using the MBH H₂ase Hyd1 from E.coli (Group 1, accession number: 3UQY PDB) showed a match with a protein from Lyngbya confervoides BDU141951 (Chandrababunaidu et al., 2015). Multiple sequence alignments indicated that the six cysteine residues (C17, C19, C20, C115, C120, and C149 in E. coli HydI) involved in the coordination of the proximal [4Fe3S] as well as the proline residue (residue 242 in HydI), both of which are typical of this class of O₂-tolerant enzymes, are conserved in the protein of L. confervoides BDU141951 (Figure 5).

The Hox enzyme of Aphanocapsa montana BDHKU210001 (Bhattacharyya et al., 2015) showed similarities with the SH H₂ase of R. eutropha (Group 3d, accession numer: AAP85843.1). The HoxH, HoxY, and HoxU proteins showed 51, 50, and 45% identity, respectively, with their respective homologs in the R. eutropha enzyme. Homologs of the Pyrococcus furiosus H₂ase SH (Group 3b) were identified in Leptolyngbya boryana PCC 6306, Cyanothece sp. PCC 7425 and Mastigocoleus testarum BC008. The sequences encoding the four subunits α (pf0894), β (pf0894), γ (pf0892), and δ (pf0893) showed an average rate of identity of 33% with those of P. furiosus. In the small subunit, the four cysteine residues serving as ligands in the coordination of the [4Fe-4S] cluster in the small subunit are conserved. In conclusion, three of the four O₂-tolerant enzymes described so far are present in the cyanobacterial phylum. Three of the strains potentially producing these enzymes are marine (Aphanocapsa montana BDHKU210001, Lyngbya confervoides BDU141951 and Mastigocoleus testarum BC008), and the other two originate from freshwater environments (Cyanothece sp. PCC 7425, and Leptolyngbya boryana PCC 6306).

The maturation process of the MBH-O₂ tolerant H₂ase of Ralstonia eutropha has been shown to involve some hox specific genes in addition to the hyp genes (Bernhard et al., 1996; Schubert et al., 2007; Ludwig et al., 2009; Fritsch et al., 2011a). The peptidase specific for this enzyme is encoded by the hoxM gene. hoxO and hoxQ genes encode for specific chaperones and hoxZ for a b-type cytochrome (Bernhard et al., 1996; Schubert et al., 2007). The maturation process of the MBH-O₂ tolerant H₂ase of R. eutropha has been shown to also involve the Hox LRTV proteins (Fritsch et al., 2011b). Homologs of the hoxZMLOQRTV genes were searched in the genome of the cyanobacterium Lyngbya confervoides BDU14195, and as a control in genomes of other organisms known to harbor the MBH-O₂ tolerant enzyme (E. coli (Evans et al., 2013), Alteromonas macleodii (Vargas et al., 2011), Hydrogenovibrio marinus DSM 11271 (Shomura et al., 2011), Rubrivivax gelatinosus (Maness et al., 2002), and Salmonella enterica (Bowman et al., 2014). The result of this analysis showed that while the hoxZ, hoxM, hoxL, hoxO, and hoxQ were conserved in all non-cyanobacterial genomes analyzed, only the hoxZ, and hoxM genes were identified in Lyngbya confervoides BDU14195 (Supplementary Table 9). The ability of this cyanobacterium to produce an active MBH-O₂ tolerant enzyme is therefore questionable. Since the maturation process of the other O₂-tolerant H₂ases found in cyanobacteria has not been reported to require any specific proteins other than the Hyp, it is possible that Aphanocapsa montana BDHKU210001, Cyanothece sp. PCC 7425 and Mastigocoleus testarum BC008 might produce active O₂-tolerant H₂ases. The genome of Leptolyngbya boryana PCC 6306 was found to contain only the hypAB genes, this strains can therefore regarded as inable to build an active O₂-tolerant H₂ase (Supplementary Table 1).

The genome set analyzed in this study includes 126 cyanobacterial genomes of the CyanoGeba dataset (Shih et al., 2013; Calteau et al., 2014), and genomes of Aphanocapsa montana BDHKU210001, Cyanobacterium sp. UCYN-A2, Lyngbya confervoides BDU141951, Mastigocoleus testarum BC008 which are present in the JGI database (https://img.jgi.doe.gov/cgi-bin/mer/main.cgi). In the case of H₂ases not generally found to occur in cyanobacteria (the [FeFe] H₂ases, and [NiFe] H₂ases other than Hox and Hup), the analysis also included cyanobacterial genomes present in the NCBI database (https://blast.ncbi.nlm.nih.gov/Blast.cgi). The complete list of the genomes analyzed and their accession numbers is given in Supplementary Table 3.

The cyanobacterial genomes present in the databases cited above were analyzed using the sequences listed in Supplementary Table 10 as queries. The e-values were adapted to the legnt of the sequences analyzed. A BLASTp (Altschul et al., 1990) analysis was conducted with a specific threshold e-value for each protein, in order to limit the number of paralogs found and therefore to avoid false positives (Supplementary Table 8). Best Reciprocal Blast Hits method and context genomic analysis were used to discriminate false positive and to choose the best the e-value threshold. Sequence alignments were carried out with Clustal-W and displayed with GeneDoc (Thompson et al., 1994; Nicholas et al., 1997). Phylogenetic analysis was performed using the Neighbor-Joining (NJ) method (Saitou and Nei, 1987) implemented in Clustalw to identify eventual false positive.

The species tree was generated by concatenating 21 conserved proteins selected from the phylogenetic markers proposed for use with bacterial genome trees (Wu and Eisen, 2008). The 21 selected proteins are: DnaG, Pgk, PyrG, RplB, RplC, RplD, RplE, RplF, RplL, RplM, RplN, RplP, RplT, RpoB, RpsC, RpsE, RpsI, RpsK, RpsS, SmpB, and Tsf. The sequences of these proteins from Anabaena variabilis ATCC 29413 were used as queries in BlastP analyses. The genomes of Chloroflexus auranticus J-10, Rhodobacter sphaeroides 2.4.1, Heliobacterium modesticaldum Ice1, and Chlorobium tepidum TLS were used as outgroups to root the tree as previously used (Calteau et al., 2014). Multiple sequence alignments of the proteins were performed using MUSCLE 3.8.31 (Edgar, 2004). The alignments were concatenated and the phylogenetic tree was generated with PhyML 3.3.2 (BioNJalgorithm/default parameters) (Guindon et al., 2009).

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.