{'fraction_non_alphanumeric': 0.03499140464703599, 'fraction_numerical': 0.018826595685687352, 'mean_word_length': 5.321998247151622, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Diversity and evolution of bacterial bioluminescence genes in the global ocean

Abstract Although bioluminescent bacteria are the most abundant and widely distributed of all light-emitting organisms, the biological role and evolutionary history of bacterial luminescence are still shrouded in mystery. Bioluminescence has so far been observed in the genomes of three families of Gammaproteobacteria in the form of canonical lux operons that adopt the CDAB(F)E(G) gene order. LuxA and luxB encode the two subunits of bacterial luciferase responsible for light-emission. Our deep exploration of public marine environmental databases considerably expands this view by providing a catalog of new lux homolog sequences, including 401 previously unknown luciferase-related genes. It also reveals a broader diversity of the lux operon organization, which we observed in previously undescribed configurations such as CEDA, CAED and AxxCE. This expanded operon diversity provides clues for deciphering lux operon evolution and propagation within the bacterial domain. Leveraging quantitative tracking of marine bacterial genes afforded by planetary scale metagenomic sampling, our study also reveals that the novel lux genes and operons described herein are more abundant in the global ocean than the canonical CDAB(F)E(G) operon.


INTRODUCTION
Marine biodiversity and evolution are intimately related with biogeography and ecology (1)(2)(3). The Tara Oceans expedition recently provided a global picture of the complex interactions between marine micro-organisms and their environment (4)(5)(6). Bioluminescence, the chemical emission of visible light, is produced by a remarkable diversity of or-ganisms and is particularly widespread in marine species (7)(8)(9). The luciferase enzymes that catalyze the emission of photons have evolved independently over 30 times, by convergence from non-luminescent enzymes (10,11). Although bioluminescent bacteria are the most abundant and widely distributed of all light-emitting organisms (7,12), certain functional and evolutionary aspects of bacterial luminescence still remain enigmatic, such as its biological role which remains a matter of debate (13). Early on bioluminescence was proposed to have evolved from ancient oxygendetoxifying mechanisms (14)(15)(16). It has also been argued that stimulation of DNA repair through the activation of DNA photolyase may confer an advantage to luminous bacteria (17), although this hypothesis is still controversial (18). Yet another hypothesis is that bioluminescence is a visual attractant for zooplankton and fish that both provide ingested bacteria with growth medium and means for dispersal (19). Symbiosis with squid or fish is also an intriguing feature of specific bioluminescent bacteria (20,21).
To date, most of the few culturable light-emitting bacterial species that have been characterized fall within the Gammaproteobacteria class. These bacteria cluster phylogenetically in three families (Vibrionaceae, Shewanellaceae and Enterobacteriaceae) which all carry a highly conserved lux operon (12). Since its first identification 40 years ago, the canonical luxCDAB(F)E(G) organization has been systematically observed in all the bacterial bioluminescent genomes (22)(23)(24). The luxA and luxB genes encode the alpha and beta subunits of the luciferase heterodimer that emits light by the oxidation of FMNH 2 and a long chain aldehyde. Whereas they both adopt a TIM-barrel fold (25), LuxA specifically displays a disordered loop playing a critical role in light emission (26). LuxC, D and E together form a fatty acid reductase complex responsible for the synthesis of the long chain aldehyde substrate (24,27).
Despite a highly conserved core, some variations have been observed in the lux operon organization. Small differences in gene content, for instance the presence of an op-tional riboflavin genes or luxF in Photobacterium species, have been observed (28,29). LuxG, which reduces FMN into FMNH2 is absent in Photorhabdus spp. whose operon also contains multiple insertions of ERIC sequences (30). Natural merodiploidy of the lux-rib operon has been also noticed in some strains in Photobacterium leiognathi (31). While the phylogeny of lux genes generally supports a vertical inheritance, multiple examples of instability of the lux locus and horizontal gene transfers (HGT) have been reported in different clades, at various taxonomic levels (32)(33)(34). Also, mutations or loss of the lux operon are frequently observed in non-luminous strains and appear to correlate with some environmental parameters (35)(36)(37)(38).
Although bioluminescent bacteria are cosmopolite in the oceans and occupy a great diversity of ecological niches, including surface and deep waters, (39)(40)(41)(42)(43), many studies have revealed an intricate relationship between bacterial bioluminescent phenotype, lux operon diversity, environmental parameters and life style (29,(44)(45)(46)(47)(48)(49)(50). Given the apparent ubiquity of bioluminescence in the ocean and the ease with which light emitting bacteria can be isolated from seawater, it has come as a surprise that bacterial bioluminescence has so far escaped detection by previous marine metagenomic studies (51). According to the pioneering authors, the unexpected absence of lux genes might have been explained by sampling protocols which filtered out size classes of potential interest, and to sequencing depth which might have been insufficient to catch bioluminescent bacteria if these were of low abundance (51).
In the present report, we surveyed the distribution of bacterial lux-related genes in a compilation of publicly available large-scale metagenomic environmental databases (Tara Oceans 2009-2013, Malaspina 2010, GOS and OSD2014) giving special care to screen the largest possible variety of organismal size sampling fractions. Spanning a wide spectrum of marine bacteria diversity, including a majority of unculturable species, our study reveals new insights about distribution, diversity and evolution of marine luxrelated genes and their operon organization at a planetary scale.  Table S1). A similar procedure was used for building the Lux B, C, D, E, G and F reference datasets.

Diversity of lux operon
Genbank bacterial genomes containing the reference and marine lux-like sequences are download from the NCBI web site (https://www.ncbi.nlm.nih.gov/genbank/). A syntheny graph representing the operon structural organizations was done using Easyfig (60).

OM-RGC Lux homologs search
The Tara Oceans OM-RGC dataset (5,61) was screened with each of the Lux HMM profiles obtained above, using hmmsearch with an E-value threshold of 1.0E-10 (Supplementary Table S1). Further filtering based on alignment lengths eliminated incomplete Lux sequences. The length thresholds were set to 340, 300, 430, 275, 300 and 200 aa for the LuxA, B, C, D, E and F homologs, respectively.

OM-RGC LuxA homolog structural filtering
Bacterial luciferases and monooxygenases share a highly conserved TIM-barrel fold (25) that renders discrimination from each other difficult from primary sequence alignments alone. We therefore developed a specific procedure to help luciferase/monooxygenase discrimination based on threedimensional (3D) structure modeling and comparison. The 3D coordinates of close structural homologs of bacterial luciferases were retrieved from PDB (52) and structurally superimposed with PyMOL (The PyMOL Molecular Graphics System, Version 1.2r3pre, Schrödinger, LLC.) (Supplementary Table S2 and Figure S2a-i). This analysis selected a set of amino acids previously described to be specific of the luciferase active site (62). In a first filtration step, truncated sequences and sequences lacking the catalytic His44 were excluded from our selection. Moreover, a disordered loop from 262-290 (V. harveyi LuxA coordinates) that covers the active site is a specific signature of LuxA proteins (63). Thus, we used Disopred v3.1 (64) to systematically calculate the disorder content in all the aligned LuxA homologs. We then filtered out sequences displaying a disorder score below 0.2 in this region (defined by the multiple sequence alignment, see below) from our final LuxA homolog dataset (Supplementary Table S3).

Phylogenetic analyses of Lux-related sequences
Multiple sequence alignments (MSA) for phylogenetic analyses were obtained with Clustal Omega (65) using default settings. In order to generate a phylogenetic tree that integrates LuxA-related, LuxB and selected monooxygenases (Supplementary Table S2), we first aligned each of these three sequence groups independently, and then processed each of the three MSAs with MaxAlign (66) to remove sequences with excess INDELS that reduce alignment area. TrimAl v1.2 with default setting (67) was then used for automated alignment trimming (removal of poorly aligned positions in the MSAs). The phylogenetic analysis was performed with RAxML (68) (parameters: -f a -# 10000 -m PROTGAMMAAUTO -p 12345 -x 12345). Trees were midpoint rooted and visualized with FigTree (https://github. com/rambaut/figtree/) and Evolview v2 (69). Alignments were visualized and inspected with Jalview (70). A similar procedure was used for the phylogenetic analyses of the other Lux proteins.

Abundance of lux genes in the global ocean
We computed the abundance of the canonical and noncanonical lux operon genes (Supplementary Table S4 and File S2) by metagenomic read recruitment from large scale marine metagenomic data: Tara Table S5). This procedure, which searches individual reads that align to a query sequence in each of the available raw unassembled metagenomes, is expected to be more sensitive than searching for homologs in the assembled metagenomes (e.g. OM-RGC). Indeed, the raw read mining procedure has the potential to detect even singleton homologous reads, whereas such singletons often don't make it past the assembly stage of reduced datasets such as the OM-RGC. Metagenomic read counts per gene (RPKM MG ) correspond to the number of mapped reads per query gene divided by the total number of reads sequenced for each sample (in millions) divided by gene length (in kilobases). The reported relative abundance values correspond to: log(1+(RPK M * 10 9 )) log (2) . Metagenomic short reads were recruited with Bowtie2 2.1.0 aligner with default parameters (78). We discarded alignments below 80% identity and below 50% of query gene horizontal coverage.

Geographic distribution of specific bacterial genomes
We analyzed the global ocean distribution of Vibrio Campbellii ATCC BAA-1116 (NCBI:txid338187) and Alcanivorax bacterium GenoA1 TS13 700 (NCBI:txid2072747, Pacific ocean: South China 700m, 19.909 N, 115.243 E) genomes by metagenomic read recruitment from the same marine metagenomes as described above (Tara Oceans, Malaspina, GOS and OSD). For each genome, we computed relative genomic abundances as the number of reads mapped onto the genome normalized by the total number of reads sequenced for each sample (Supplementary Table  S6) (79). We identified and removed from the read recruitment counts 267 and 528 outlier genes on V. Campbellii and Alcanivorax spp. respectively, to reduce apparent biases due to atypical genomic regions with low sequence specificity. We generated the world maps with R version 3.

Biogeography and taxonomic distribution of the OM-RGC luxA homologs
We used the Ocean Gene Atlas web service (80) (http://taraoceans.mio.osupytheas.fr/) to describe the geographic and taxonomic distributions of the 3 groups of OM-RGC LuxA homologs (OLAHs). Abundance and location were visualized using the Tara Oceans metagenomic samples corresponding to size fractions 0.22-1.6 and 0.22-3 m sampled either at surface or mesopelagic depths.

Diversity of lux gene homologs
In order to survey luxA, B, C, D and E homologs in the global ocean, we first inventoried all Lux sequences deposited in public reference databases (NR, RefSeq, UniProtKB/Swiss-Prot) to produce a set of reference lux gene sequences (Supplementary Table S1). As expected for genes forming an operon, a similar number of occurrences of each lux genes was found in public reference databases. Surprisingly, this search also revealed a set of lux homologs in the reference genomes of six bacteria that had not previously been described as bioluminescent: Alcanivorax spp., Rhizobacter spp. (Gammaproteobacteria); Enhygromyxa salina, Plesiocystis pacifica (Deltaproteobacteria); Leptospira santarosai (Spirochaetes); Actinomyces spp. (Actinobacteria) (Supplementary Table S1). These six bacterial reference genomes possess luxA, C, D and E of the lux operon, with the notable exception of luxB. This thus revealed the presence of lux-like operons in three bacterial classes never described as capable of producing light: Deltaproteobacteria, Spirochaetes and Actinobacteria. Genes in these luxB-less operons will henceforth be referred to as lux-like genes.
We then used this set of reference sequences to search for lux homologs in the most recent microbial gene catalog of Tara Oceans (OM-RGC), which represents assembled genetic material recovered from the free-living mostly bacterial size fractions (0.2-3 m) collected in the global ocean (Supplementary File S1). Consistent with the current knowledge of bacterial bioluminescence distribution, we detected the presence in the OM-RGC of lux genes from several well-known bioluminescent bacteria such as Vibrio campbeli, Vibrio cholerae, Aliivibrio fischeri, Shewanella woodyi and P. leiognathi (Figure 1 and Supplementary Table  S1).
An initial search for LuxA proteins provided a set of 800 marine candidates having both high similarity score and size close to our reference LuxA sequences (Supplementary Figure S1). Further filtering according to the presence of a correctly positioned structurally disordered region as well as the presence of specific active site residues generated a set of 401 OM-RGC LuxA homologs (OLAHs) (Supplementary Table S3).
By combining LuxA, LuxA-like and LuxB reference sequences together with the 401 marine OLAHs and some selected monooxygenases of known structure, we inferred a phylogenetic tree of LuxA homologs (Figure 1). Phylogeny and disorder content of LuxA sequences are highly congruent (Supplementary Figure S3) as indicated by the clustering of bona fide reference LuxA sequences that display a similar disordered loop. The LuxA reference sequences are grouped together in a monophyletic branch well separated from the the LuxB monophyletic branch, both of which are then connected to a group of homologs we annotated as LuxA-like sequences. The bulk of the marine homologs then cluster into three branches, forming OLAH groups 1-3 ( Figure 1). However, OLAH groups 1 and 2, as well as LuxA-like sequences all share a similar consensus in the catalytic site with LuxA reference sequences (Supplementary Figure S4). The monooxygenases and OLAH group 3 appear to form a separate branch, coherent with a distinct active site logo signature (Supplementary Figure  S4). Interestingly, the few LuxA-like sequences belonging to bacteria that have not been previously described as luminescent organisms form a separate but close branch to the LuxA reference group (bootstrap value of 100). We should note that the LuxA sequences of Alcanivorax spp. and P. pacifica found in the OM-RGC dataset are grouped with the LuxA-like branch ( Figure 1). Two monooxygenase sequences (pdbid: 3i7g and 4uwm) are found within OLAH groups 2 and 3, respectively (Figure 1 and Supplementary  Table S2). Phylogenetic analysis of marine LuxB, C, D, E homologs similarly show distinct groups consisting of Lux references, Lux-like and OM-RGC distant relatives (Supplementary Figure S5).

Diversity of lux operons
Reference and marine lux-like sequences differed not only in their phylogeny, but also in their operon organization. Indeed, we observed operon structural organizations very different from the canonical CDAB(F)E(G) order ( Figure  2 and Supplementary Figure S6). In particular, our study revealed a novel CEDA order in which luxB is missing, in both Gammaproteobacteria and Deltaproteobacteria. Moreover, the Leptospira spp. that belongs to the more distant Spirochaete class encodes a CAED operon, while the Actinomyces AxxCE lux-like operon structure is even more atypical. Importantly, we observed a congruence between the operon organization and the phylogenetic distance in the LuxA sequence tree (Figure 2). Thus, our findings extend the known diversity of the lux operon organizations and reveal the existence of non-canonical lux-like operons in new taxons.
The similarities between homologous genes in the different operons have been also analyzed (Figure 2 and Supplementary Figure S7). As expected, each lux gene is highly similar to its relative in close taxons. Interestingly, the genomes of the two distinct myxobacteria genera Enhygromyxa and Pleisiocistis (81,82), exhibit full lux operons more than 90% identical. High similarities are also found between the Alcanivorax and Rhizobacter lux like-operons (61.1%). These findings suggest a lateral transfer of the CEDA lux-like operon between these species.

Abundance and distribution of lux-related genes in the global ocean
In order to investigate the distribution of the lux-related genes (in canonical or non-canonical operons), all reference and OLAH sequences were used as queries to recruit raw metagenomic reads from public large scale marine metagenomes (Tara Oceans, Malaspina, GOS and OSD) (Supplementary Table S5; see materials and methods). In addition, to differentiate lux genes likely harboured in free living bacteria from those likely associated with particles or larger organisms (e.g. through symbiosis), we analyzed their relative abundances in samples obtained by filtering through filters of pore sizes of either 0.22-3 m (free-living) or 0.8-2000 m (associated). The abundance of lux genes are summarized in Figure 3 and detailed in Supplementary Table S7. We have also estimated the abundance of the canonical lux operons of all the well-known bioluminescent bacteria (Figure 2). Rather surprisingly, the luxoperon of V. campbelii was the only detectable known bioluminescent bacterium, which furthermore was detected in Tara Oceans samples only, thus suggesting its higher abundance relative to other bioluminescent bacteria. Reference lux genes were not observed in either free-living or associated size fractions collected during the Malaspina expedition (which systematically sampled bathypelagic layers), nor the GOS expedition, nor on OSD2014 solstice day (mostly coastal waters). An important finding was that the reference V. campbelii luxA gene was preferentially detected in associated size fractions (i.e. non free-living) collected in just five Tara Oceans coastal and surface samples (TARA005, TARA007, TARA008 in the Mediterranean sea, and TARA123, TARA125 close to the Marquise islands) (Supplementary Table S7). This shows that the V. Campbellii canonical lux operons are mainly associated with particles or larger planktonic species. In contrast, Figure 3 shows that the genes of the Alcanivorax non-canonical CEDA operon are abundantly distributed in both Tara Oceans and Malaspina samples, in all size fractions (from 0.2 to 2000 m), both in surface and deep ocean depths.
In order to further investigate these intriguing findings, we extended this lux only metagenomic read recruitment analysis to the distribution of Vibrio spp. and Alcanivorax spp. full length genomes in the global ocean. Figure  4 summarizes their geographic, water column and freeliving/associated distributions (quantitative data are reported in Supplementary Table S6). Vibrio genomes are widely distributed in temperate regions, in surface, meso- LuxA-like in red; OM-RGC LuxA Homologs (OLAH) group 1 (128 sequences), 2 (209 sequences) and 3 (63 sequences) from light to dark blue. Reference monooxygenases clustered together in a separate branch are colored purple; the two purple circles correspond to PDB monooxygenase sequences, 4uwm (13) and 2i7g (14) that branch in OLAH groups 3 and 2, respectively. The pink circles correspond to the OM-RGC assigned as LuxA reference sequences of Vibrio campbellii (1), Vibrio cholerae (2), Aliivibrio fischeri (3), Shewanella woodyi (4) and Photobacterium leiognathi (5).The red circles correspond to the sequences we annotate as LuxA-like of Plesiocystis pacifica(6) and Alcanivoracaceae bacterium (7).The yellow circles correspond to the OM-RGC assigned as LuxB reference sequences of S. woodyi (8), A. fischeri (9), P. leiognathi (10), V. campbellii (11), V. cholerae (12). The origin of the sequences used for tree inference are presented in Supplementary Table S3, and OLAH amino acid sequences are provided in fasta format (Supplementary File S1). and bathypelagic depths and are not observed in polarregions ( Figure 4A). In contrast to the preference of individual lux genes for associated fractions (Figure 3), Vibrio full genomes were detected in both free-living and associated fractions, albeit significantly more abundant in the free-living bacterial fractions than in the associated fractions (Mann-Whitney U-test P < 0.0002) ( Figure 4B). Consequently, our study suggests that free-living Vibrios tend to be devoid of the lux operon.
Alcanivorax genomes showed a very different distribution ( Figure 4C): they were highly abundant (more than 10 times the abundance of vibrio genomes), and except for one station (West of the North Atlantic gyre), they were exclusively localized in tropical and sub-tropical zones of the Pacific Ocean. They were present in a wide range of depths from surface to bathypelagic layers. Similarly to Vibrio, Alcanivorax genomes were present in both the free-living and the associated-size fractions ( Figure 4D), which, contrary to observations with Vibrios, was consistent with the distribution of the Alcanivorax lux-like operon (Figure 3).
To gain further insight in the distribution of sequences corresponding to OLAH group 1 (Figures 5A and B) and OLAH groups 2-3 ( Supplementary Figures S6 and 7) of the LuxA phylogenetic tree (Figure 1 and Supplementary Figure S3), we also analyzed their ocean distributions in the Tara Oceans samples. Group 1 sequences showed a particularly interesting distribution along the water column. Al-though they were observed in all oceans and latitudes, they were systematically much more abundant in mesopelagic (about 200-1000 m) than in surface depths ( Figure 5C). These sequences clustered into three known groups of bacteria: Entotheonella spp., Actinobacteria and Proteobacteria ( Figure 5D). Interestingly, the differences in abundances between surface and mesopelagic depths gradually decreased from group 1 to group 3 ( Figure 5C; Supplementary Figures  S8c and 9c) and a more diversified taxonomic distribution was also observed in group 3 ( Figure 5D; Supplementary  Figures S7d and 8d).

DISCUSSION
Many studies have shown that bacterial luminescence is related to environmental parameters and/or life-style, such as symbiotic associations (35)(36)(37)(38). Bacterial bioluminescence has so far been observed as a canonical lux operon that follows the CDAB(F)E(G) gene order, in genomes of gramnegative bacteria that group phylogenetically in three families of Gammaproteobacteria: the Vibrionaceae, Enterobacteriaceae and Shewanellaceae (12). However, minor variations around this common theme have been observed in various species (12,(28)(29)(30)(31). Our deep exploration of public genomic and metagenomic marine databases expands the lux operon catalog (Figure 2). It also reveals a larger distribution of lux operons among bacterial taxons than previously  Supplementary Table S5.  Table S6). Stations with at least two distinct lux genes are circled in black. Boxplots showing the Log10 of the relative abundance of (B) Vibrio spp. and (D) Alcanivoracacae at different depths. Blue boxes represents free living size fractions and red boxes the associated size fractions. Genomes and samples used for metagenomic read recruitment are listed in Supplementary Tables S4 and 5. known, as well as a highly variable lux operon organization repertoire (Figure 2). We describe the existence of a new lux-like CEDA operon, lacking the luxB gene, in both Gammaproteobacteria and Deltaproteobacteria (83). In Gammaproteobacteria, the CEDA lux-like operon configuration is observed in Pseudomoniales (Rhizobacter) and Oceanospirales (Alcanivorax) (84). Since the canonical CDAB(F)E(G) form of the lux operon is highly conserved in three Gammaproteobacteria families (7,12), it is somewhat surprising to find a distinct lux gene order among the Gammaproteobacteria. Remarkably, the finding of two different lux operon organizations within the Gammaproteobacteria appears to be congruent with the phylogeny of this clade. Indeed, according to Gammaproteobacteria species phylogeny (83), the CDAB(F)E(G) and CEDA lux-like operon configuration correspond to the divergence of two distinct groups, the PO (Pseudomonadales and Oceanospirillales) and VAAP/Entero (Vibrionales, Alteromonadales, Aeromonadales and Pasteurellales/Enterobacterials), respectively ( Figure 2). Despite an absence of homology around the operon in the genomes of the two Gammaproteobacteria families harbouring the CEDA lux-like operon, the lux-like genes display a high similarity, suggesting a lateral transfer of the lux operon between the two remote Rhizobacter and Alcanivorax genera (Supplementary Figure S7). Such   Supplementary Table S3. lateral transfer might have occured in rhizosphere microbiomes shared by these two species (85). In addition to marine environments (84) alcanivorax spp. has been indeed observed in coastal areas where it thrives in oil-polluted soils and rhizospheres (85).
In Deltaproteobacteria, we have also found the noncanonical CEDA lux operon form in two myxoccocale marine bacteria (86), Enhygromyxa salina (81) and P. pacifica (82). The close homology between the CEDA operon of these two distant species that share the same ecological niche strongly supports a lateral transfer origin (81,82) (Supplementary Figure S7). Remarkably, two distant bacterial classes share the CEDA configuration: the Gammaproteobacteria and the Deltaproteobacteria. Our findings cannot yet distinguish between an HGT or/and a vertical inheritance of this operon configuration from a common ancestor of these two classes. In the latter case, the CEDA configuration would correspond to an ancestral lux operon that predated the luxA/luxB duplication, congruent with the LuxA inferred phylogenetic tree. On the another hand, an HGT scenario is not uncommon between distant taxons when two species share a common niche (87).
We also show that Leptospira spp. belonging to a bacterial taxon phylogenetically distant from proteobacteria, the Spirochaetes, also carries an even more distinctive lux operon configuration: the CAED form. Leptospira genera occupy a great variety of habitats (88) and display highly plastic genomes consistent with frequent HGT events (89). A yet more distant form AxxCE is observed in Actinomyces, a gram-positive Actinobacteria (90). Interestingly, it has been reported that some Mycobacteria spp. contain a locus, mel2, very similar to the lux operon (91)(92)(93). Strikingly, the mel2 locus organization follows a similar AxxCE order as found here in Actinomyces. Consistent with the absence of the luxB gene from the Actinobacter lux operon, the analysis of the relatedness of MelF to luxA and luxB placed MelF on an independent branch equally related to both of them (91). Our results suggest that the non-canonical lux-like operon of gram-positive bacteria may correspond to an ancestral form and support the hypothesis of Kasai et al. suggesting that the lux operon originated in gram-positive bacteria (32).
In order to better understand the link between the genetic diversity, evolution and the ecology of the lux operon, we have searched independently each lux genes in the global ocean and collected information about their geographic distribution, their localization in the water column and their presence in either free-living or associated size fractions. Two distinct scanning procedures of the metagenomic data provided complementary results.
First, we investigated the presence of lux genes in the assembled metagenomic data (OM-RGC) by protein homology search (Figure 1), thus providing a qualitative view of the presence of lux genes in the global ocean. Previous analyses did not recover lux sequences from assembled metagenomic DNA of planktonic marine microorganisms (51), thus suggesting that luminous Photobacterium, Aliivibrio and Vibrio were present at abundances so low they went undetected with afforded sequencing efforts (29). In the present study however, the higher sequencing deepness detected several well-known bioluminescent bacteria such as V. campbellii, V. cholera, A. fischeri, S. woodyi and P. leiognathi. We infer that the other well-known marine bioluminescent bacteria that are not detected here are either not represented in the Tara Oceans sampling locations, or in abundances too low to allow their assembly. Interestingly, the OM-RGC assembly does contain lux-like genes belonging to the non-canonical CEDA lux-like operon of Alcanivorax spp. (84) and P. pacifica (82), two species that have been recently characterized as marine bacteria.
Second, we quantitatively explored the biogeographic abundances of known reference lux operons by metagenomic read recruitment, a method that provides a measure of geographic abundance of lux operons in the global ocean. Only V. campbellii and Alcanivorax spp. lux genes were detected, indicating either that they have a higher abundance than other known bacteria containing lux or the lux-like genes, or that the niche habitats of other bioluminescent bacteria were not sampled by Tara Oceans, Malaspina, GOS and OSD. The comparison of their full genome distributions in the global ocean shows contrasting stability of the lux operon in V. campbellii compared to Alcanivorax spp.. The Vibrios appear to contain an accessory lux operon preferentially found in symbiotic or particle associated strains localized in coastal environments (Figures 3 and 4). This is consistent with previous studies showing that non-luminescent Vibrios lacking the lux operon are free-living (35,36,38). In contrast, most of the detected Alcanivorax strains stably possess the CEDA-lux operon in their genomes and display a homogenous distribution restricted to the Pacific Ocean. Alcanivorax spp. are ∼10× time more abundant than V. campbellli, and are slightly more frequently observed in free living fractions (Figures 3 and 4). The two distinct patterns of bacterial/lux operon relationships suggest that the different lux operon configuration may have different functions and genomic stability. Although several Alcanivorax spp. have been characterized as hydrocarbon degrading bacteria, to our knowledge, their bioluminescent phenotypes are not documented.
To obtain finer grain resolution of bacterial bioluminescence evolution and biogeography, we carried out a comprehensible analysis of LuxA homologs, carefully selected on the basis of both sequence homology, structural and disorder properties. Interestingly, disorder content is congruent with phylogeny and shows distinct patterns in the three well separated OLAH groups (Figure 1 and Supplementary Figure S3). We believe that this set of sequences may provide new insights to decipher luciferase evolution pathways. LuxA belongs to the group C family of two component monooxygenases, which adopt a very similar TIM-barrel fold and probably share a common ancestor (94)(95)(96). However, a comparison of LuxA and the closest non-luminous monooxygenase structures (pdb id: 4us5 and 4uwm) (97,98) revealed that their active sites have both different sequences and structures (Supplementary Table S2 and Figure S2ai). A first interesting result is the discovery of two OLAH groups (1 and 2) that share many common features with reference LuxA that seem to distinguish them from monooxygenases ( Supplementary Figures S3 and 4). Remarkably, the phylogeny and the ocean distribution of OLAH sequences also show a congruent pattern. Indeed, sequences of the cosmopolite OLAH group 1 are systematically much more abundant in mesopelagic layers ( Figure 5A-C). In contrast, the two OLAH groups 2 and 3 that include known reference monooxygenases are uniformly distributed along the water column (Supplementary Figures S8a-d and S9a-d). The fact that the group 1 is preferentially found in dark environments is coherent with the hypothesis that they could represent a new class of luciferases. An alternative hypothesis explaining higher abundance deep in the ocean could be a specific adaptation to high-pressure of luminescent or nonluminescent proteins. Interestingly, more than one third of OLAH group 1 sequences are assigned to Entotheonella. These bacteria have been recently described to form symbiotic associations with sponges and are characterized by an over-representation of FMN-monooxygenase genes (99), that are closely related to luciferases. Such reservoirs of FMN-monooxygenase genes may have contributed to the evolution of luciferases.
How and why bacterial luciferases have evolved to catalyze light-emission is still an open question. To date, the diversity of bioluminescent bacteria is predominantly based on the study of culturable Gammaproteobacteria harbouring a highly conserved lux operon organization. Our metagenomic study provides new insights about lux operon diversity and distribution in the global ocean, as well as clues for future exploration on luciferase evolution and functions in marine bacteria including unknown and uncultivable species.

SUPPLEMENTARY DATA
Supplementary Data are available at NARGAB Online.
The bulk of genomic computations were performed using the OSU Pythéas cluster with the help of the SIP members and the Airain HPC machine provided through GENCI-[TGCC/CINES/IDRIS] (grants t2011076389, t2012076389, t2013036389, t2014036389, t2015036389 and t2016036389). Tara Oceans (which includes both the Tara Oceans and Tara Oceans Polar Circle expeditions) would not exist without the leadership of the Tara Expeditions Foundation and the continuous support of 23 institutes (http: //oceans.taraexpeditions.org). We further thank the commitment of the following sponsors: CNRS (in particular Groupement de Recherche GDR3280 and the Research Federation for the study of Global Ocean Systems Ecology and Evolution, FR2022/Tara Oceans-GOSEE), European Molecular Biology Laboratory (EMBL), Genoscope/CEA, The French Ministry of Research and the French Government 'Investissements d'Avenir' programmes, FRANCE GENOMIQUE, MEMO LIFE and PSL* Research University. We also thank the support and commitment of agnès b. and Etienne Bourgois, the Prince Albert II de Monaco Foundation, the Veolia Foundation, Region Bretagne, Lorient Agglomeration, Serge Ferrari, Worldcourier and KAUST. The global sampling effort was enabled by countless scientists and crew who sampled aboard the Tara from 2009-2013, and we thank MERCATOR-CORIOLIS and ACRI-ST for providing daily satellite data during the expeditions. We are also grateful to the countries who graciously granted sampling permissions. The authors declare that all data reported herein are fully and freely available from the date of publication, with no restrictions, and that all of the analyses, publications and ownership of data are free from legal entanglement or restriction by the various nations whose waters the Tara Oceans expeditions sampled in. This article is contribution number 101 of Tara Oceans. Conflict of interest statement. We declare no conflict of interest and confirm the material is original research and has not been previously published or submitted for publication elsewhere.