Sophie Mathieu

Mélanie Touvrey

Laurent Poulet

Sophie Drouillard

Selvi Thirumalai

Laure Ulaganathan

Miroslav Ségurel

William Cygler

Helbert

Thirumalai Selvi Ulaganathan

Laure Ségurel

Miroslav Cygler

William Helbert

The porphyran degradation system of the human gut microbiota is complete, phylogenetically diverse and geographically structured across Asian populations

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The impact of diet on the diversity and composition of the bacterial communities inhabiting the human gut is now well documented 1,2,3 . For example, several studies have highlighted the lower bacterial diversity of populations having westernized diet, as compared to hunters-gatherers or populations with more diversified food 4,5 . Conversely, it is also fascinating to observe that the gut microbiota can adapt and is able to acquire new functional capacities, notably by lateral transfer of genetic materials between strains of the microbiota but also from strains associated with the food ingested by the host 6,7 . One of the most illustrative example of how the food participates in shaping the human gut microbiota, at the molecular level, is the discovery of the lateral transfer of genes encoding enzymes involved in the degradation of porphyran, the cell wall polysaccharide of the red algae Porphyra sp. used to prepare the maki-sushi, from a marine strain living at the surface of the algae to the gut bacteria Bacteroidetes plebeius found in the Japanese gut microbiota 8 . Since this discovery, other algal polysaccharides degradation systems (i.e. agarose, alginate, carrageenan), likely acquired by lateral transfer from marine organisms, have also been studied 9,10,11 .

Porphyran is a agar type polysaccharide made of two disaccharides repetition: agarobiose and porphyranobiose decorated up to 64% by methyl groups 12,13 (Figure 1A). In B. plebieus, the genes encoding the enzymes for the degradation of porphyran are co-localized and co-regulated in a so-called "polysaccharide utilizing loci" (PUL). Many enzymes of the PUL (i.e. glycoside hydrolases and sulfatases) have been characterized, explaining the degradation process of the agarose components and the non-methylated fractions of porphyran by B. plebieus 8,14,15 (Figure 1B), the enzymes involved in the degradation of methylated porphyran segments remaining to be discovered.

Even though the fronds of Porphyra sp. are nowadays mostly consumed by Asian populations for its taste and nutritional value 16 , the most ancient proof that humans harvested and stored Porphyra sp. was found in the archaeological site of Monte Verde (Chile) where population occupied 12,300 years ago 17 . The traditional use of Porphyra spp. by the first nations living along the Pacific coast of North America was also documented by ethnobotanists [START_REF] Goddard | Life and culture of the Hupa[END_REF][START_REF] Coulter | Plants used by the Indians of the Mendocino county, California[END_REF]20,[START_REF] Turner | The ethnobotany of edible seaweed (and related species; Rhodophyta: Bangiales) and use by First Nations on the Pacific coast of Canada[END_REF] . In South America, pre-Inca and Inca populations also incorporated red algae in their diet [START_REF] Michanek | Seaweed resources for pharmaceutical uses[END_REF] . The very ancient know-how and use of Porphyra sp., preserved for millennia along the Pacific Ocean's coast, suggests that several lateral transfer events of genetic material dedicated to porphyran degradation from marine bacteria to gut bacteria might have occurred in ancient human populations. The origin of the algal degradation system in the human gut and the timing of these transfers have yet to be estimated. Are there traces of ancestral algal catabolism acquisitions still conserved in the contemporary human gut? How recent are the lateral transfers currently observed in Asian populations?

We observed that the B. plebeius PUL dedicated to the degradation of porphyran encoded enzymes have a larger spectrum of specificities than previously reported 8,14,15 and can notably process the methylated fraction of porphyran. The discovery of novel enzymes that can accommodate the methyl group of porphyran substrate in their active sites demonstrates that the complete degradation of all the components of porphyran (e.g. agarose, methylated and nonmethylated porphyran fractions) can be achieved autonomously by the enzymes encoded by the PUL. Moreover, we have found that this self-sufficient PUL is present not in one but in many bacterial strains of the gut suggesting that several lateral transfers have occurred between bacteria of the gut. We noticed several mutations and recombination events between individuals suggesting the acquisition of the porphyran PUL was not recent. To further characterize the genetic differences between populations living in different geographical areas, we examined genomic and metagenomic data available to identify and locate individuals having the porphyran degradation machinery.

Results

Porphyran PUL, the PUL dedicated to the porphyran degradation, identified in B. plebieus contained genes encoding two -porphyranases (BpGH16B, Bacple_01689; BpGH86A, Bacple_01693) and two -agarase (BpGH16A; Bacple_01670; BpGH86B, Bacple_01694) 8,14 (Figure 1B). The oligosaccharides resulting from degradation of porphyran by these enzymes are then degraded further by exo-acting enzymes including sulfatase (BpS1_11, Bacple_01701), -L-galactosidase (BpGH29, Bacple_01702) and -Dgalactosidase (BpGH2C, Bacple_01706) produced by B. plebeius or another agarolytic strain such as Bacteroidetes unifromis 15 (Figure 1B). The presence of a gene (BpGH50, Bacple_01683) encoding an enzyme grouped in the GH50 family known to contain agarases, was also studied but its substrate specificity was not demonstrated 11 . Altogether, our knowledge on the catabolism of porphyran is limited to its non-methylated fraction. We continued the detailed analysis of the porphyran PUL of B. plebeius and studied the uncharacterized glycoside hydrolases BpGH16C (Bacple_01703) and the BpGH2C, completed the structural data of the sulfatase BpS1_11, and demonstrated the porphyran PUL encodes a portfolio of enzymes catalyzing to the complete degradation of the methylated fraction of the polysaccharide.

Biochemical and molecular characterization of the endo-acting 6O-methyl-porphyranase

BpGH16C (Bacple_01703)

The gene Bacple_01703 encodes a glycoside hydrolase of the GH16 family (BpGH16C), which includes characterized agarases in sub-families GH16_15-16, and porphyranases in sub-families GH16_11-12 [START_REF] Viborg | A subfamily roadmap of the evolutionarily diverse glycoside hydrolase family 16 (gh16)[END_REF] . BpGH16C was classified in the GH16_14 subfamily that contains one characterized -agarase (Vibrio sp. strain PO-303, BAF62129.1) [START_REF] Dong | Molecular cloning, expression, and characterization of a beta-agarase gene, agaD, from a marine bacterium, Vibrio sp. strain PO-303[END_REF] and one endo-6O-methyl--porphyranase (Wenyingzhuangia fucanilytica, WP_068825734.1) [START_REF] Zhang | Characterization of a novel porphyranase pccommodating methyl-galactoses at its iubsites[END_REF] .

We assayed BpGH16C on various marine polysaccharides including carrageenan, agarose and porphyran. Only porphyran was degraded, leading to the production of a series of oligosaccharides characteristic of endo-acting glycoside hydrolase (Figure 2A). The structural characterization by NMR of the end-products purified by chromatography [START_REF] Correc | Structural analysis of the degradation products of porphyran digested by Zobellia galactanivorans βporphyranase A[END_REF] revealed the occurrence of a methyl group on the -linked D-galactose residue (Figure 2B), demonstrating that BpGH16C is an endo-6O-methyl--porphyranase able to accommodate the methyl group of porphyran in its active site. For comparison, we have examined the substrate specificities of one another predicted GH16_14 endo-6O-methyl--porphyranase from Paraglaciecola atlantica T6c (Patl_0824, ABG39352.1) presenting 51.1% identity with the BpGH16C and one GH16_12 -porphyranase (Patl_0805, ABG39333.1) homologous to BpGH16B. Analyses of the degradation products confirmed the different porphyranase specificities (Figure 2A). We determined the crystal structure of BpGH16C at 1.9 Å resolution (PDB ID 8EP4) (Figure 3A). The enzyme has a β-sandwich jelly-roll fold with two stacking antiparallel β-sheets typical of proteins from GH16 family [START_REF] Viborg | A subfamily roadmap of the evolutionarily diverse glycoside hydrolase family 16 (gh16)[END_REF] and like in most GH16, a Ca 2+ is bound on the convex side of the β-sandwich [START_REF] Naretto | The agar-specific hydrolase ZgAgaC from the marine bacterium Zobellia galactanivorans defines a new GH16 protein subfamily[END_REF] . One molecule of Hepes from the crystallization buffer bound in the center of the cleft, with the 2-hydroxyethyl end buried into the molecule interior and the sulfate group exposed on protein surface (Figure 3A). To get a better understanding of substrate binding and catalysis, we crystallized the E145L inactive mutant with the tetrasaccharide L-α-6O-sulfate-Gal-(13)-D-β-Gal-(14)-L-α-6O-sulfate-Gal-(13)-D-a/β-Gal and determined its structure at 1.8 Å resolution (PDB ID 8EW1). An electron density was present within the groove in all the three independent molecules in the asymmetric unit fitted a D-β-Gal-(14)-Lα-6O-sulfate-Gal-(13)-D-β-Gal trisaccharide (Figure 3B). The electron density for the D-Gal on the nonreducing end is weaker than for the first two residues and no density is observed for the fourth residue indicating its mobility in the crystal. The residues correspond to the -1, -2 and -3 positions according to the standard nomenclature [START_REF] Davies | Nomenclature for sugar-binding subsites in glycosyl hydrolases[END_REF] .

The trisaccharide sits edge-on in the groove with C5 substituent of the ring in the -1 site directed toward the bottom of the cleft while O1 and O2 hydroxyls point toward the solvent.

Trp134 stacks against hydrophobic side of the -1 residue while Arg67 and Glu238 hydrogen bond to the O4 hydroxyl. Residue in the -2 subsite is stacked between Trp64 and the edges of Trp143 and Phe169. Finally, the sulfate group of the -2 residue is hydrogen bonded to Asn61, Ser136 and Asn234 through a bridging water molecule. The D-Gal on the non-reducing end (-3 site) has somewhat different orientation in molecules A/B and C and makes only one bridging hydrogen bond through a bridging water (Figure 3C). The structure of the complex with a trisaccharide adds also to our understanding of the catalytic mechanism. Glu150 is indeed the closest to the O1 hydroxyl that formed the bond to the next sugar eliminated by hydrolysis, while Glu145 approaches from the opposite side of the ring relative the O1 hydroxyl (Figure 3C) and is likely helping in proper orientation of the substrate. Asp147 is directed to the bridging O1 and might be helping in catalysis.

Structural alignment of BpGH16C with three other GH16 porphyranases: porphyranase A (PDB ID 3ILF) and porphyranase B from Zobellia galactanivorans (PDB ID 3JUU) 8 and another porphyranase from B. plebeius (PDB ID 4AWD) 14 shows that the Gly132 is replaced by either a Ser or a Thr, making this pocket too small to accommodate an additional methyl group on the C6-O (Figure 3D). Sequence alignment of -porphyranases with BpGH16C and the two other characterized 6O-methyl--porphyranases: Patl_0824 of (P. atlantica T6c, this study) and AXE80_06940 (W. fucanilytica) [START_REF] Zhang | Characterization of a novel porphyranase pccommodating methyl-galactoses at its iubsites[END_REF] confirmed that the replacement of the Ser/Thr by Gly was correlated to the accommodation of the methyl group in the -1 sub-site of the catalytic site (Figure S1). The BpGH16C is the first crystallized 6-OMe-porphyranase that tolerates the 6-methyl group on the galactose ring. Indeed, there is an extra space in the active site at the end of the pocket that accommodates C6-OH, with its bottom formed by Gly132 (Figure 3E).

Complete degradation of oligo-methyl-porphyrans

Following the same strategy reported by Robb and co-workers 15 , we have incubated the methylated oligosaccharides with the BpS1_11 sulfatase, the BpGH29 -L-galactosidasee and the predicted BpGH2C-D-galactosidase. Purified methylated and non-methylated oligoporphyrans were incubated with the BpS1_11 sulfatase and analyses of the products by 1 H NMR showed a strong chemical shift (about 0.1 ppm) of the signal attributed to the proton at position 6 of the galactose located at the non-reducing end, demonstrating the removal of the sulfate ester group independently of the degree of methylation of the oligosaccharides (Figure S2).

Crystal structure of the BuS1_11 sulfatase from B. uniformis complexed with a disaccharides 15 and the B. plebeius sulfatase BpS1_11 on its own (this study, PDB ID 7SNJ, 1.64 Å resolution) and complexed with a tetrasaccharide (this study, PDB ID 7SNO, 2.1 Å resolution) revealed that there is room in the active site to accommodate oligo-methylporphyrans (Figure S3). The porphyran sulfatase BpS1_11 belongs to the formyl-glycine dependent sulfatases, requiring maturation of a cysteine or a serine into the formyl-glycine residue as an active catalytic residue. A large N-terminal domain of BS1_11 contains active site residues and forms main part of the substrate binding site. A smaller C-terminal domain completes the substrate binding site.

To map the details of substrate binding, we determined the structure of H214N inactive

mutant with L-α-6O-sulfate-Gal-(13)-D-β-Gal-(14)-L-α-6O-sulfate-Gal-(13)-D-a/β-Gal
tetrasaccharide substrate (PDB ID 7SNO). The electron density map clearly showed the entire tetrasaccharide bound to BpS1_11 (H214N). Following the recently proposed nomenclature [START_REF] Hettle | The molecular basis of polysaccharide sulfatase activity and a nomenclature for catalytic subsites in this class of enzyme[END_REF] , the tetarasaccharide binds to the 0, +1, +2 and +3 positions (Figure S3). The 6S-L-Gal in position 0, which the sulfate is removed, has all three hydroxyl groups hydrogen bonded to the protein side chains: O2 to Arg271, O3 to Asp215 and Arg271, O3 to Asp345 and Arg347. The 6-sulfate group is facing Ser83 (formyl-glycine in mature enzyme) and one of its oxygens is a ligand to the Ca 2+ ion. D-Gal in +1 position hydrogen bonds C6OH to NE2 of His133 and its O2 is bonded to the O2 of L-Gal (position 0) through a bridging water molecule. The +2 6S-L-Gal makes only one contact, between the 6-sulfate oxygen and His428. The +3 D-Gal extends far from the binding site and makes no contacts with the protein. However, in the crystal, this sugar hydrogen bonds from O1 and O4 as well as O5 within the ring to the guanidino group of Arg68 from a symmetry-related molecule. These interactions stabilize the +2 and +3 sugars in the crystal. Importantly, the 6-O group of +1 D-Gal is not tightly constraint by the enzyme and there is sufficient space for the additional methyl group present in the 6-O-Me-D-Gal (Figure S3).

We confirmed the exo-based cycle of porphyran depolymerization 15 using tetrasaccharides as starting substrate (Figure S4). After desulfation of the tetra-saccharides, BpGH29 exo--L-galactosidase was active on the methylated and non-methylated substrate demonstrating that the methyl group did not hinder the activity of the enzyme. Independently of Robb and co-workers, 15 we have determined the structure of BpGH29 (PDB ID 7SNK). The superposition with their structure (PDB ID 7LJJ) shows root-mean-squares deviation of 0.6 Å, indicating that the structures that were crystallized under different conditions are virtually identical, and confirms that ordering of the loop 408-426 and small rearrangement of other surrounding loops result from the substrate binding to the active site. Finally, we observed that the methylated or non-methylated D-galactose residue located at the non-reducing end of the trisaccharide was cleaved by the BpGH2C -D-galactosidase (Figure S4). Altogether, the enzymology and crystallography experiments conducted in this work revealed the pathway leading to the degradation of the methylated fraction of porphyran.

Our observations combined with previous investigations demonstrate that all the methylated and non-methylated fractions of porphyran are digested by the enzymes encoded by porphyran PUL of B. plebeius. Therefore, this PUL can catabolize autonomously the chemically complex porphyran, without the help of enzyme(s) external to this PUL.

Genetic diversity and bacterial occurrence of the porphyran utilizing loci

Previous transcriptomic experiments conducted on B. plebeius grown in the presence of porphyran in the culture medium interrogated the genes Bacple_01668 to Bacple_01706 14 .

Interestingly, the transcriptomic results suggested that at least two groups of genes were differently regulated. The cluster Bacple_01692 to Bacple_01699 genesnamed PUL-PorA in Figure 1B was moderately up-regulated in the presence of porphyran. This contrasts with the two neighboring clusters of genes: Bacple_01668 to Bacple_01689 (PUL-PorC) and Bacple_01700 to Bacple_01706 (PUL-PorB, Figure 1B), which were highly up-regulated (10fold more than PUL-PorA). The PUL-PorB genes encoding the enzymes studied herein (Figure 1B) were strongly up-regulated, in agreement with their implication in porphyran metabolism 14 .

Until now, B. plebeius DSM 17135 is the only gut bacteria shown to carry the porphyran PUL. In order to investigate whether other human gut bacteria can also digest porphyran, we used BLASTn to identify homologs of the entire coding sequences of the B. plebeius DSM 17135 PUL-PorB among the genomes of isolated human gut bacteria (more than 10,000 non redundant genomes). We found that the PULs were present in the genomes of 26 Phocaeicola/Bacteroides strains (corresponding to at least 8 different identified species), all coming from strains isolated from the gut microbiome of Asian individuals (Table S1). The same species/strains isolated from other human populations worldwide had no homologs of porphyran PULs. The identified PULs were present in different species, including B. dorei, B.

eggerthii, B. ovatus, B. plebeius, B. stercoris, B. uniformis, B. vulgatus and B. xylanisolvens,

but were absent in other common Bacteroides sister species from the human gut, such as B. intestinalis or B. fragilis. This phylogenetic diversity probably reflects numerous lateral transfer events between bacteria of the gut microbiome. An important consequence of this observation is that the taxonomic composition of an individual human is not necessarily a good predictor of its metabolic capacity.

While the PUL-PorB organization as observed in B. plebieus DSM 17135 was conserved in all the other strains we have identified, the gene organization of PUL-PorA was only conserved in 19/26 strains (Figure S5A). This genomic region was incompletely sequenced in B. plebeius strain AM09-36 and B. eggerthii strain AM42-16. In four other strains (B. stercoris strain AF05-4, B. sp. AF25-38AC B. stercoris strain AM32-16LB, B. xylanisolvens strain AF38-2) the homolog of the B. plebieus -agarase (GH86, Bacple_01694) appeared to have recombined with a DNA fragment composed of two genes predicted to encode metabolic enzymes: 2-dehydro-3-deoxygluconokinase and 2-dehydro-3-deoxyphosphogluconate aldolase (Figure S5B). For the same four strains, a deletion of the genes encoding the GH50 (Bacple_01683) and GH105|GH154 (Bacple_01684) located in PUL-PorC was also observed, highlighting the remodeling of the PUL by insertion and deletion of DNA segments (Figure S5B).

The Neighbor-Joining tree calculated with the six concatenated genes of the PUL-PorB (Figure 4A, based on 12,888 bp) revealed three clades with high bootstrap values, allowing to classify the strains in three groups accordingly to their PUL-PorB genetic diversity (Figure 4A).

In the first group (GI, 18 strains, including the B. plebeius DSM 17135 reference), all six genes presented 99-100% identity with the reference, except for B. plebeius strain AM49-7BH, which had only 97-98% identity for GH29 and GH16. In the second group (GII, six strains), the six genes were 99-100% identical to each other, but presented 94-99% identity with the B. plebeius reference. The genes of the B. plebeius strain AM49-7BH suggested that it resulted from a recombination event that occurred between SusD, SusC and GH2 genes from GI and Sulf and GH29 genes more related to GII. The recombination site was identified in the gene coding for GH16 (Figure S6). For this reason, we classified B. plebeius strain AM49-7BH as GIrec. The gene sequences of the B. plebeius strain AM09-36 and B. stercoris strain AF05-4 were more distantly related to the other strains and were classified as a third group (GIII). The genes of the B. stercoris strain AF05-4 presented 92-96% identity with the B. plebeius reference. On the other hand, the B. plebeius strain AM09-36 sequence seemed to result from a recombination event between GI and GIII (Figure 4A), since the SusD, SusC and GH2 genes presented 99-100% identity with GI genes, while the S1_11 sulfatase and GH29 genes were nearly identical to GIII genes. Again, the recombination site was identified in the gene coding for GH16 (Figure S6). For this reason, we classified B. plebeius strain AM09-36 as GIIIrec.

Distribution of PUL-PorB in assembled Human gut metagenomes

We explored the worldwide geographic distribution of the porphyran degradation system in diverse human gut microbiota by probing the six PUL-PorB genes against assembled human gut metagenomic datasets (Table S2). We scanned 19 projects (9 Asian, 7 North American, 1 South American and 2 African datasets), covering in total 2720 biosamples for which gut metagenomes were obtained and assembled. Only fully sequenced genes were retrieved. We found that at least one fully sequenced gene of PUL-PorB was present in the assembled gut metagenomes of 370 East Asians out of the 1144 individuals testedcorresponding to 36% of the Chinese individuals and 20 % of the Japanese individuals, while we found only ten positive individuals among the 1262 North American gut metagenomes (0.8% of the individuals) and none among the other geographic areas, including 117 individuals from South Asia (India, Bangladesh), 47 from South America (Peru) and 150 from Africa (Tanzania, Madagascar) (Table S2). Therefore, as observed in previous investigations, the occurrence of algal polysaccharides degradation systems seems to be mostly restricted to East Asian populations.

The gene and the corresponding individuals were characterized as belonging to GI/GIrec, GII, GIII and/or GIIIrec based on their gene sequence (Table S3). For most individuals (92.6%), the genes detected were attributed to a single group; however, few individuals (5.9%) presented PUL-PorB genes that belonged to two different groups. The individuals were then grouped by geographic location (Figure 4B). In order to obtain large sets of individuals, we distinguished populations living in North-East China (Hangzhou), South-East China (Shenzen) and Japan. In both set of Chinese individuals, the group GI/GIrec was present in more than half of the population (56% and 53%, respectively) followed by the group GII (33% and 39%, respectively). The group GIIIrec was only present in Chinese individuals, and more specifically, in North-East China were it was occurring in 15% of the individuals. In contrast, in Japan, GI and GII were found in similar proportion (42%), and the GIII group represented 20% of individuals. North-East China, South-East China and Japan populations were overall found to be significantly different (X-squared=49.78, p-value=2.5x10 -5 ). This difference was mostly driven by a significant difference between the Japanese population and each of the Chinese ones (q-value=0.0008 and 0.0323 for Japan versus Hangzhou and Shenzen, respectively, after FDR Bonferroni correction for multiple testing), while the two Chinese populations were not significantly different from one another (q-value=0.1280). We then tested separately which group was different in prevalence between populations and found that GIII and GIIIrec were significantly different (proportion test, q-value=0.0009 and 0.0406, respectively), while GI and GII were not (proportion test, q-value=1 and 0.8925, respectively).

Analysis of short read datasets

The analysis of the assembled metagenomes revealed that 92.6% of the individuals carried only one version of the PUL-PorB. This observation may however result from some bias inherent to the bioinformatics processes involved in the building up of the contigs.

Therefore, to analyze the data devoid of such biases, we have selected a 54 nucleotides sequences allowing to probe the different PUL-PorB groups directly on raw sequencing data (e.g. short reads). The probe was part of the Bacple_01703 (GH16) gene encoding the methylporphyranase, which presented characteristic mutations for the GI, GIrec, GII, GIII and GIIIrec groups (Figure 5). This method allowed us to explore new metagenomic datasets which were not assembled, thus expending the number of sampling sites, as well as the number individuals tested, with now a total of 2313 individuals (Table S2).

As observed previously, we confirmed by this independent approach that PUL-PorB is quite common in East Asian coastal populations (54% in China and 67% in Japan). We note that the prevalence of positive individuals is much higher with this approach than with the assembled datasets, with also observed a higher variance across populations within Japan (4-92%) than within China (20-69%). The histograms presented in Figure 5 show the number of short reads and their attribution to the GI, GIrec, GII, GIII and GIIIrec groups fort each individual. Among positive individuals, we found that most of them (80.0%) carried short reads attributed to a single group. Interestingly, when individuals presented two types of short reads (or more), one group was clearly predominant (the median ratio of the first highest number of reads to the second highest number of reads was 4). Consequently, we characterized individuals as belonging to a group if they carried only this group, or if this group was dominant, which we defined as having more than twice the amount of reads for this group compared to the second most prevalent group. Using this definition, we found that 93.8% of individuals carried short reads attributed to a single dominant group, a number very similar to what was observed with the assembled datasets.

Using these categories, analyses of the short reads confirmed that GI and GIrec mutations were predominant in Chinese populations (54.2% of individuals in total, versus 40.3% in Japan) while the GII ones were most common in Japanese population (51.0% of individuals, versus 42.4% in China). Individuals with the GIII mutation were found mainly in the Japanese population (13.8%) and in lower amount in Chinese population (2-6.7%). Finally, GIIIrec mutations were the most numerous in the Hangzhou population (9.3%), and in much less common elsewhere (0.3-5%).

Including only categories with large enough numbers (i.e., dominant GI, GIrec, GII, GIII, GIIIrec and co-dominant GI+GII groups), we found that the five populations were overall significantly different (X-squared= 93.2, p-value=2.0x10 -11 ) (Figure 6). This difference was mostly driven by a significant difference between the Japanese population and three of the Chinese ones (all except Beijing, q-value<0.004), as well as between Hongkong and Hangzhou (q-value=3.7x10 -7 ). We then tested separately which group was significantly different in prevalence between populations and found that all were significantly different (proportion test, q-value<0.015) except GIrec (proportion test, q-value=1), but GIII and GIIIrec presented the stronger difference (proportion test, q-value=2.2x10 -6 and 3.8x10 -8 , respectively), similarly to the MG results.

Discussion

Biochemical and crystallographic investigations showed that the porphyran PUL encodes tailored enzymes dedicated to the complete degradation of the polysaccharide, including the methylated fraction, giving galactopyranoside and 6-O-methyl-galactopyranoside as end-products. We showed that the 6-O-methyl-porphyranase (BpGH16C) presents an active site able to accommodate the methyl group, in contrast with the active sites of the previously investigated -porphyranases. Similarly, the exo-6O-L-galactose porphyran sulfatase Bp S1_11 can accommodate methylated oligosaccharides in its active site. The demethylation of 6-Omethyl-galactopyranoside was located in agarose PUL of several aerobic bacteria belonging to the genus Gammaproteobacteria, Bacteroidetes and Planctomycetes. The reaction was obtained with specific monooxygenase enzyme system, including ferredoxin, ferredoxin reductase and P450 monooxygenase [START_REF] Reisky | Oxidative demethylation of algal carbohydrates by cytochrome P450 monooxygenases[END_REF] . However, demethylation of galactose residue in anaerobic bacteria has not been elucidated yet and no encoding gene surrounding the porphyran utilization loci could be suspected to be involved in the demethylation pathway.

Because harvesting, drying and cooking Porphyra sp. were conserved practices in numerous populations living along the Pacific and because these ancestral traditions seem to follow the dispersion line of Homo sapiens from Asia to the South of the American continent, one could expect that the porphyran degradation system may be present in populations living along the Pacific coast. Analysis of metagenomics data of population in Peru or Fidji (Table S2) revealed no degradation system resembling the Asian porphyran PUL. In this context, analyses of the distribution and history of the porphyran PUL was restricted to the Chinese and Japanese populations. Homologous gene sequences of the PUL-PorB found in populations living by the sea were also present in continental population of China, showing that the PUL-PorB has spread among the Chinese population independently of the geographical location.

One hypothesize that the porphyran PUL was acquired, at first, by people living by the sea and then was spread among the continental population. The genes were conserved in the gut microbiota of the individual as long as the diet, including algae, maintained a selective pressure.

Investigations of the metagenomic data were conducted using assembled and raw (e.g. short reads) data. The analyses demonstrated that the investigated populations were structured with a limited number of genetic variations that were used to characterize 5 groups (GI, GIrec, GII, GIII and GIIIrec). Exploitation of the raw data allowed us to estimate a higher number of individuals carrying strains involved in the porphyran degradation, highlighting the likely better sensitivity of this approach. Also, analyses of the short reads revealed that one PUL variant was present or was dominating in individuals, giving an apparent haploid organization of the PUL.

The architecture of the microbiota at the bacterial level seems to drive the structuration at genetic and molecular level, or reciprocally.

The genetic organization and the dispersion of the various PUL-PorB suggested that they have diverged from an ancestral gene cluster. However, there is no evidence that the evolution from an ancestral PUL occurred in the human gut or that homologous PUL, which evolved independently were acquired at different periods. The low number of genetic variations of the PULs and their wide distribution suggests that acquisition by lateral transfer from marine bacteria living at the surface of red algae toward strains of the human gut remains a rare event.

Searching marine strain genomes and metagenomics data (e.g. Tara Ocean) failed to identify porphyran PUL similar or distantly related to those observed in the human gut. Therefore, there is yet no environmental data that could help to estimate the pace and number of lateral transfer events between marine bacterial strains and human gut microbiota.

Because the PUL-PorB GI and GII groups were common in the Chinese and Japanese populations and because they were identified in many different human gut strains, one can suppose that these clusters of genes were present in the early human gut microbiota, allowing their dispersion in the Asian populations and lateral exchanges between gut microbiota strains.

The occurrence of the GIIIrec group was found mostly in the Hangzhou population and in lower abundance in the Shenzen and Beijing populations, while almost absent in the Japanese population, suggesting a more recent appearance than the GI and GII groups. It is interesting to note that this recombinant strain is currently found in populations with the lowest prevalence of the parent strain GIII (i.e., China), but not in the only population with a non-negligible frequency of GIII (i.e., Japan). This might suggest that after having appeared, the GIIIrec might have locally replaced the GIII group, which bodes the question of the functional equivalence of these strains. In general, the rate of dispersion of the PUL variants inside a given population and between populations geographically distant may probably reflect exchanges between populations which involved displacement of individuals.

In conclusion, the architecture of the human gut microbiota across populations has mostly been studied and functionally interpreted at a genera or species level, while we see here with the example of the porphyran degradation system that the important functional information is carried at the strain, or even the genetic level. Indeed, the PUL system is carried by various species and strains of the Bacteroides genus, and the GI/GII/GIII groups are not associated with particular strains, so the inference of these groups based on taxonomy would not be possible.

Although it is still difficult to state which mutations were obtained by independent evolution inside the human gut population and which mutations were acquired by lateral transfer, in this study, we demonstrated that human populations geographically distant presented common (e.g. GI, GII) and specific genetic characteristics (e.g. GIrec, GIII and GIIIrec). The current structure of the investigated populations is likely the result of lateral transfer, recombination, as well as migration events, which reflect an evolution history of the gut microbiota and therefore, of its host. buffer pH 5.5 and 30% MPD; BpGH29 (Bacple_01702) was crystallized at 18 mg/ml in 20% Peg 8K and 0.1 M KH2PO4; BpGH16C (Bacple_01703) was crystallized at 34 mg/ml in 16% Peg 8K, 0.1 M Hepes pH 7.5 and 0.2M Calcium acetate. For diffraction experiments, the crystals were cryo protected in reservoir solution containing 20% MPD (BpS1_11), 25% glycerol (BpGH29), 30% ethylene glycol (BpGH16C) and flash frozen in liquid nitrogen. The diffraction data of all the crystals were collected at the 08ID beamline at the Canadian Light Source 34 .

The X-ray diffraction data were processed using XDS 35 . The structures were solved by molecular replacement using the program Phaser in Phenix package 36 . For BpS1_11, the structure solution was obtained by molecular replacement using the search model an arylsulfatase from Pseudomonas aeruginosa (PDB ID 1HDH) 37 . For BpGH29, the structure solution was obtained using the α-L-fucosidase from Fusarium graminearum as a search model (PDB ID 4NI3) 38 . BpGH16C was solved using porphyranase B from Zobellia galactanivorans (PDB ID 3JUU) 39 as a model. All the structures were refined with the Phenix software 40 and manual rebuilding and solvent placement was conducted with the COOT program 41 . The stereochemistry of all the models were validated with MolProbity 42 .

To obtain the complex structure, putative active site mutants of BpS1_11 (BpS1_11(H214N)) and BpGH16C (BpGH16C(E145L)) were made using the Quickchange site-directed mutagenesis protocol and using KOD polymerase. Briefly, the plasmid containing the gene of interest was amplified with the primer pairs carrying the mutation using KOD polymerase. After PCR, the template plasmid was digested with 1 µl of DpnI enzyme for an hour at 37°C. Five µl of PCR product was then transformed into 50µl of chemically competent E. coli DH5α cells. The clones carrying the desired mutation were confirmed by sequencing and proceeded for crystallization experiments.

Based on the structural comparison with other sulfatases, the His 214 residue in BpS1_11 was mutated to Asn (H214N). The mutant protein BpS1_11(H214N) was expressed and purified following the same protocol used for wild type enzyme. Purified Bacple_01701(H214N) was crystallized from solution containing 40% PEG 200, 0.1M sodium citrate buffer pH 5.5 and 30% MPD. The putative active site mutant of BpGH16C(E145L) was purified following the same protocol used for wild type enzyme. BpGH16C(E145L) was crystallized from solution containing 20% PEG 8K, 0.1 M Hepes pH 7.5 and 0.2 M calcium acetate. BpS1_11(H214N) and BpGH16C(E145L) crystals were soaked in the tetrasaccharide solution for an hour before the diffraction experiments. The soaked crystals were flash frozen in liquid nitrogen and diffraction data was collected at the 08ID beamline in the Canadian light source 34 . The diffraction data was processed with XDS 35 . One round of rigid body refinement was carried out using phenix refinement program (Adams et al., 2011). Manual rebuilding and substrate placement were done using coot (Emsley et al., 2010). The geometry was validated using MolProbity program 42 .

Human gut metagenomes analysis

More than 10,000 non redundant genomes of bacteria isolated from the human gut were were downloaded from various bacterial genome repositories including NCBI 43 , PATRIC 44 , references strains of the Human Metagenome Project 45 , the Human Gastrointestinal Bacteria Culture Collection (HBC) 46 , the Culturable Genome Reference (CGR) 47 and human Gut Microbial Biobank (hGMB) 48 . Bacterial strains harboring the porphyran PUL were identified after BLASTn of PUL-PorA (17,980 bp), PUL-PorB (12,888 bp) and PUL-PorC (19,198 bp) sequences against the set of genome. Genes of the identified porphyran PUL variant were retrieved from bacterial genome repositories. Genes and concatenated genes were aligned using Muscle and distant tree were built using MEGA6 allowing to distinguish GI, GIrec, GII, GIII and GIIIrec group of genes.

Assembled metagenomes data set (Table 2) were downloaded and analyzed locally.

Genes of the porphyran PUL were BLASTn and aligned with Muscle using MEGA6. Only coding genes with full length were kept and compared with those of the reference bacteria. The genes encoding the sulfatase were not used because of its low discriminative power. We did not discriminate between GI and GIrec because most of the time, there was not enough information.

To identify GIIIrec, we required to observe genes matching both GI and GIII on the same contig. When only part of the genes was available, we considered the major group (i.e., GI or GIII) as the most probable one, creating a slight bias downward for GIIIrec. When multiple contigs in the same individual were found to be positive and they corresponded to different groups, we considered the individual as multi-groups (e.g., GI+GII). When we could not assign a contig to one of the 4 groups defined earlier, we considered the individual to be unresolved.

Distant trees were built allowing to group the genes in GI, GIrec, GII, GIII or GIIIrec.

Analyses of short reads data set (Table S2) were conducted using 54 nt probes (Figure S7). The probes were used to search at Sequence Read Archive (SRA) data at NCBI. showing the active site region. The single mutation of Ser129/Thr137 (PorA/PorB) (sidechains shown as spheres) to Gly132 in BpGH16C provides space for accommodating a methyl group at C6 on L-galactose residue. E) A cros-section of the surface representation of the Gal sugar and the cavity in the binding site near Gly132. There is free space within the cavity sufficient to accommodate the methyl group of a methylated porphyran. S8 obtained from the analyses of short read sequencing of Chinese and Japanese metagenomics data. For 80% of the individuals, only one group of porphyran PUL was observed. For the other individuals, a dominant group was considered when the number short-reads detected was, at least, twice more abundant than the other. In other situation, the co-occurrence of groups was reported.

Captions for figures

Figure 1 :

 1 Figure 1: A) Chemical structure of the repetition moieties of porphyran. The sulfated disaccharides -porphyranobiosescan present methyl group at the position 6 of the Dgalactose giving methylated and none-methylated porphyran components. Agarobiose units, obtained by desulfation/cyclization of the L-galactose residue during the biosynthesis, are methylated accordingly. B) Organization of the B. plebieus porphyran PUL. The PUL is divided in three segments (Pul-PorA, -Por B and -PorC) based on previous transcriptomic analyses. The enzyme functions were determined in (1) Hehemann et al., 2012, (2) Giles et al., 2017, (3) Rodd et al., 2022 and (4) this study.

Figure 2 :

 2 Figure2: A) Size exclusion chromatography of the degradation products of porphyran incubated with the methyl-6O--poprhyranase BpGH16C (GH16_14). The degradation profile was compared with the two P. atlantica T6c porphyranases grouped in the GH16_12 (Patl_0824) and GH16_14 (Patl_0805) sub-families. B)1 H NMR of the disaccharides endproducts of BpGH16C.

Figure 3 :

 3 Figure3: A) Crystal structure of BpGH16C (Bacple_01703). Chain A of asymmetric unit is shown in cartoon representation and colored in spectrum from N to C terminal. The HEPES molecule in the active site is shown as stick representation. The N-and C-terminals are marked, and the strands are labeled in their order in the sequence. B) The substrate binding site in the mutant BpGH16C-E145L complexed with the D--Gal-(14)-6O-sulfate-L-Gal disaccharide end-product. The semitransparent surface of the protein is colored by the electrostatic potential. The residues forming the site are shown in stick representation. The sulfate group from L-galactose 6 sulfate is docked into a positively charged pocket. C) The hydrogen bonds between the disaccharide in substrate binding site residues. Several H-bonds are bridged by water molecules (W). Glu145 (in blue and thicker bonds) from the native structure is superimposed on the Leu145 in the mutants. D) Structure comparison of BpGH16C (green), PorA (PDB id-3ILF) (magenta) and PorB (PDB id-3JUU) (yellow) showing the active site region. The single mutation of Ser129/Thr137 (PorA/PorB) (sidechains shown as spheres) to Gly132 in BpGH16C provides space for accommodating a methyl group at C6 on L-galactose residue. E) A cros-section of the surface representation of the Gal sugar and the cavity in the binding site near Gly132. There is free space within the cavity sufficient to accommodate the methyl group of a methylated porphyran.

Figure 4 :

 4 Figure 4: A) Phylogenetic tree calculated with the concatenated gene sequences of PUL-PorB. Sequence homology observed allowed to identify recombination sites and to create groups of homologous porphyran PUL. B) Distribution of the different group of porphyran PUL in human gut metagenome assembly of South east and North East Chinese, and Japanese populations.

Figure 5 :

 5 Figure 5: The 54 nucleotides sequence characteristic of different group of porphyran PUL (top) were Blastn against the short read sequencing of metagnomic data recorded on Chinese and Japanese population. The number of detected short reads are indicated for each individuals and were grouped as a function of the geographic location of the sampling.

Figure 6 :

 6 Figure 6: Distribution of the various porphyran PUL among Chinese and Japanese population. The pie charts derive from the histogram presented in FigureS8obtained from the analyses of short read sequencing of Chinese and Japanese metagenomics data. For 80% of the individuals, only one group of porphyran PUL was observed. For the other individuals, a dominant group was considered when the number short-reads detected was, at least, twice more abundant than the other. In other situation, the co-occurrence of groups was reported.

Figure S1 :

 S1 Figure S1: Structural alignment of the B. plebieus DSM 17135 BpGH16C 6-O-methyl-porphyranase compared with the other characterized 6-O-methyl--porphyranases (P. atlantica T6c, Patl_0824; W. fucanilytica AXE80_06940) and characterized -porphyranases (B. uniformis_NP1, BuGH16A; P. atlantica T6c, Patl_0805; B. plebeius DSM 17135, BpGH16B). The blue circle indicates the serine/threonine amino acids observed in the active site of the -porphyranases which are replaced by a glycine amino acid in 6-O-methyl-porphyranases allowing to accommodate the methyl group.

Figure S2: 1 H

 1 Figure S2:1 H NMR recorded before and after incubation with the sulfatase BpS1_11 on purified disaccharides end-products of the BpGH16C. The occurrence of the methyl group located at the position 6 of the D-galactose didn't hinder the removal of the sulfate ester group of the L-galactose residue located at the non-reducing end.

Figure S3 :

 S3 Figure S3: A) Structural details of Bacple_01701. Schematic representation of Bacple_01701 shown in cartoon representation and colored from N to C terminus. Calcium ion is shown as green sphere. B) Conserved active site residues from S1 family sulfatases. The residues are shown in stick representation in the order Bacple_01701/ 5G2V/1HDH/6BIA. C) Surface representation of H214 mutant with the tetrasaccharide, showing pocket like architecture of active site. The 6S sulfate subsite is labelled as S. D) Complex structure of H214N mutant. The electron density (2fo-fc) map for the sugar contoured at 1.5σ. Sugars numbered from 0 to +3 and 6S sulfate subsite labelled as S. Residues at H bonding distance to the substrate are shown in stick representation. E) Structural comparison of H214 mutant with 1HDH. The residues involved in catalysis S83 (overlays with DDZ from 1HDH), H133, H214 (this is the mutated residue) are shown in stick representation.

Figure S4: 1 H

 1 Figure S4: 1 H NMR recorded after sequential incubation of the methylated tetrasaccharide end-products of the 6-O-methyl-porphyranase BpGH16C with the sulfatase BpS1_11 followed by the -L-galactosidase BpGH29 and the D-galactosidase BpGH2C.

Figure S5 :

 S5 Figure S5: Gene organization of the porphyran degradation system of B. plebieus compared with other identified human gut bacteria. A) Phylogenetic tree was calculated using concaneted gene sequences of PUL-PorB. PUL organizations of each strains were indicated based on the available sequencing data. B) Comparison of the two main gene organizations of the porphyran PUL observed in the set of bacteria investigated highlighting deletion and recombination events.

Figure S6 :

 S6 Figure S6: Recombination sites observed in the GH16 of PUL-PorB of the group GI/distant GII giving GIrec and the groups GI/GIII giving GIIIrec.

Figure S7 :

 S7 Figure S7: Phylogenetic trees calculated for all gene sequences of the members of the PUL-PorB.

Figure 1 :

 1 Figure 1: A) Chemical structure of the repetition moieties of porphyran. The sulfated disaccharidesporphyranobiosescan present methyl group at the position 6 of the D-galactose giving methylated and none-methylated porphyran components. Agarobiose units, obtained by desulfation/cyclization of the L-galactose residue during the biosynthesis, are methylated accordingly. B) Organization of the B. plebieus porphyran PUL. The PUL is divided in three segments (Pul-PorA, -Por B and -PorC) based on previous transcriptomic analyses. The enzyme functions were determined in (1) Hehemann et al., 2012, (2) Giles et al., 2017, (3) Rodd et al., 2022 and (4) this study.

Figure 2 :

 2 Figure 2: A) Size exclusion chromatography of the degradation products of porphyran incubated with the methyl-6O--poprhyranase BpGH16C (GH16_14). The degradation profile was compared with the two P. atlantica T6c porphyranases grouped in the GH16_12 (Patl_0824) and GH16_14 (Patl_0805) sub-families. B) 1 H NMR of the disaccharides end-products of BpGH16C.

Figure 3 :

 3 Figure 3: A) Crystal structure of BpGH16C (Bacple_01703). Chain A of asymmetric unit is shown in cartoon representation and colored in spectrum from N to C terminal. The HEPES molecule in the active site is shown as stick representation. The N-and C-terminals are marked, and the strands are labeled in their order in the sequence. B) The substrate binding site in the mutant BpGH16C-E145L complexed with the D--Gal-(14)-6O-sulfate-L-Gal disaccharide end-product. The semitransparent surface of the protein is colored by the electrostatic potential. The residues forming the site are shown in stick representation. The sulfate group from L-galactose 6 sulfate is docked into a positively charged pocket. C) The hydrogen bonds between the disaccharide in substrate binding site residues. Several Hbonds are bridged by water molecules (W). Glu145 (in blue and thicker bonds) from the native structure is superimposed on the Leu145 in the mutants. D) Structure comparison of BpGH16C (green), PorA (PDB id-3ILF) (magenta) and PorB (PDB id-3JUU) (yellow) showing the active site region. The single mutation of Ser129/Thr137 (PorA/PorB) (sidechains shown as spheres) to Gly132in BpGH16C provides space for accommodating a methyl group at C6 on L-galactose residue. E) A cros-section of the surface representation of the Gal sugar and the cavity in the binding site near Gly132. There is free space within the cavity sufficient to accommodate the methyl group of a methylated porphyran.

Figure 4 :

 4 Figure 4: A) Phylogenetic tree calculated with the concatenated gene sequences of PUL-PorB. Sequence homology observed allowed to identify recombination sites and to create groups of homologous porphyran PUL. For each strain, the identity of the PUL-PorB genes was compared with the corresponding gene of B. plebeius DSM 17135 B) Distribution of the different group of porphyran PUL in human gut metagenome assembly of South east and North East Chinese, and Japanese populations.

Figure 5 :

 5 Figure 5: The 54 nucleotides sequence characteristic of different group of porphyran PUL (top) were Blastn against the short read sequencing of metagnomic data recorded on Chinese and Japanese population. The number of detected short reads are indicated for each individuals and were grouped as a function of the geographic location of the sampling.

Figure 6 :

 6 Figure 6: Distribution of the various porphyran PUL among Chinese and Japanese population. The pie charts derive from the histogram presented in Figure S8 obtained from the analyses of short read sequencing of Chinese and Japanese metagenomics data. For 80% of the individuals, only one group of porphyran PUL was observed. For the other individuals, a dominant group was considered when the number short-reads detected was, at least, twice more abundant than the other. In other situation, the co-occurrence of groups was reported.

Figure S1 :

 S1 Figure S1: Structural alignment of the B. plebieus DSM 17135 BpGH16C 6-O-methyl-porphyranase compared with the other characterized 6-O-methyl--porphyranases (P. atlantica T6c, Patl_0824; W. fucanilytica AXE80_06940) and characterized -porphyranases (B. uniformis_NP1, BuGH16A; P. atlantica T6c, Patl_0805; B. plebeius DSM 17135, BpGH16B). The blue circle indicates the serine/threonine amino acids observed in the active site of the -porphyranases which are replaced by a glycine amino acid in 6-O-methyl--porphyranases allowing to accommodate the methyl group.

Figure S2: 1 H

 1 FigureS2:1 H NMR recorded before and after incubation with the sulfatase BpS1_11 on purified disaccharides end-products of the BpGH16C. The occurrence of the methyl group located at the position 6 of the D-galactose didn't hinder the removal of the sulfate ester group of the L-galactose residue located at the non-reducing end.

Figure

  Figure S3: A) Structural details of Bacple_01701. Schematic representation of Bacple_01701 shown in cartoon representation and colored from N to C terminus. Calcium ion is shown as green sphere. B) Conserved active site residues from S1 family sulfatases. The residues are shown in stick representation in the order Bacple_01701/ 5G2V/1HDH/6BIA. C) Surface representation of H214 mutant with the tetrasaccharide, showing pocket like architecture of active site. The 6S sulfate subsite is labelled as S. D) Complex structure of H214N mutant. The electron density (2fo-fc) map for the sugar contoured at 1.5σ. Sugars numbered from 0 to +3 and 6S sulfate subsite labelled as S. Residues at H bonding distance to the substrate are shown in stick representation. E) Structural comparison of H214 mutant with 1HDH. The residues involved in catalysis S83 (overlays with DDZ from 1HDH), H214 (this is the mutated residue) are shown in stick representation.

Figure S4: 1 H

 1 Figure S4: 1 H NMR recorded after sequential incubation of the methylated tetrasaccharide endproducts of the 6-O-methyl-porphyranase BpGH16C with the sulfatase BpS1_11 followed by the -Lgalactosidase BpGH29 and the D-galactosidase BpGH2C.

Figure S5 :

 S5 Figure S5: Gene organization of the porphyran degradation system of B. plebieus compared with other identified human gut bacteria. A) Phylogenetic tree was calculated using concaneted gene sequences of PUL-PorB. PUL organizations of each strains were indicated based on the available sequencing data. B) Comparison of the two main gene organizations of the porphyran PUL observed in the set of bacteria investigated highlighting deletion and recombination events.

Figure S6 :

 S6 Figure S6: Recombination sites observed in the GH16 of PUL-PorB of the group GI/distant GII giving GIrec and the groups GI/GIII giving GIIIrec.

Figure S7 :

 S7 Figure S7: Phylogenetic trees calculated for all gene sequences of the members of the PUL-PorB.

Material and methods

Purification of porphyran

6 g of dried Porphyra columbina were suspended in 120 ml de distilled H2O and autoclaved for 30 min at 120°C. After 14 h at room temperature, the suspension was centrifuged at 8000 rpm during 30 min at 4°C. The supernatant was added to 120 ml of pure ethanol (50% v/v EtOH final concentration) and the solution was maintained at 4°C for 2h. After centrifugation (8000 rpm, 30 min, 4°C), the porphyran, present in the supernatant, was precipitated by addition of 130 ml of pure ethanol (67% v/v EtOH final concentration). After centrifugation the porphyran pellet was dissolved in distilled water, and dialysis for 3 days against distilled water using dialysis membrane with a cut-off 1000 Da (pre-wetted Spectra/Por® 6 dialysis tubing -Sprectrumlabs). The polysaccharide was lyophilized and the purity was verified by 1 H-NMR. The yield of purification was about 15-20 % w/w dried algae.

Heterologous expression of Bacteroides plebeius DSM 17135 PUL-PorB enzymes

The genes encoding the predicted glycoside hydrolases (BpGH29, Bacple_01702; BpGH16C, Bacple_01703; BpGH2C, Bacple_01706) and sulfatase (BpS1_11, Bacple_01701) from B. plebeius DSM 17135 were cloned using genomic DNA as template in the pET-28 or pFO4 expression plasmid 32 without their signal peptides identified with SignalP 33 . The expression strains harboring the recombinant expression plasmids were grown in Luria Bertani (LB) medium supplemented with 50 µg/ml kanamycin (pET28a plasmid) or 100 µg/ml ampicillin (pFO4 plasmid) until the OD600nm reached 0.6 in a shaking incubator working at 180 rpm and 37°C. After the addition of isopropyl--D-thiogalactopyranoside (IPTG), the temperature was cooled down at 20°C and maintained overnight.

Cultures (200 mL) were centrifuged at 6000 g for 15 min and the bacterial pellet was suspended in 10 ml of buffer A (20 mM Tris pH 8, 500 mM NaCl, 20 mM Imidazole). The cells were lysed using a cell disrupter (Constant system Ltd). Insoluble fractions were removed by centrifugation at 20000 g during 30 min at 4°C and the supernatant was loaded on a 1 mL HisTrap TM HP column (GE Healthcare) connected to a NGC chromatography system (BioRad).

The proteins were eluted with a imidazole gradient from 20 to 300 mM gradient of imidazole.

Pure enzymes were obtained after a size exclusion chromatography using a HiLoad 16/600 Superdex 75 pg column in buffer B (10 mM Tris-HCl pH 8, 50 mM NaCl).

Enzymatic assays

(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Enzymatic degradations were monitored by analytical gel permeation chromatography using Superdex S200 10/300 and Superdex peptide 10/300 (GE Healthcare) columns mounted in series and connected to a high-performance liquid chromatography (HPLC) Ultimate 3000 system (Thermo Fisher). The injection volume was 20 µL and the elution was performed at 0.4 mL.min-1 in 0.1 M NaCl. Oligosaccharides were detected by differential refractometry (Iota 2 differential refractive index detector, Precision Instruments).

The oligosaccharide products were purified by semi-preparative gel permeation chromatography using three HiLoad® 26/600 Superdex® 30 pg (GE Healthcare) columns mounted in series and connected to a semi-preparative size-exclusion chromatography system which consisted of a Knauer pump (pump model 100), a refractive detector (iota2 Precision instrument) and a fraction collector (Foxy R1) mounted in series. The elution was conducted at a flow rate of 1.2 mL.min -1 at room temperature using 100 mM (NH4)2CO3 as eluent. The collected fractions were freeze-dried prior NMR and mass spectrometry analyses.

Samples were exchanged twice with D2O and were transferred to a 5 mm NMR tube. 1 H NMR spectra were recorded at 323 K using an Advance III 400 MHz spectrometer (Bruker).

Chemical shifts are expressed in ppm in reference to water. The HOD signal was not suppressed. The homogenous fractions of proteins obtained from gel permeation chromatography were concentrated and crystallization experiments were attempted. Initial crystals were obtained by screening against wide range of commercial screens and the hits were optimized by hanging drop diffusion method. The drop containing 1 l of protein and 1 l of reservoir solution was incubated over 1 ml of reservoir solution and crystal growth was monitored regularly. The conditions of the best diffracting crystals of the three proteins were listed: BpS1_11 (Bacple_01701) was crystallized at 23 mg/ml in 40% Peg200, 0.1 M Sodium Citrate (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Table S1: List of the bacterial strains which assembled genome contained encoding genes of the porphyran degradation system. Bacteroides and Phocaeicola are basonyms. The name of the strains were reported as reported in databank.

GenBank assembly accession

Assembly level

Biosamples Strains

GCA_003466465.


Mohyeddine Omrane

Thomas J Melia

Abdou Rachid Thiam

LC3 conjugation to lipid droplets

Keywords: Lipid droplets, noncanonical autophagy, LC3B, ATG3, prolonged starvation, organelle biology, membrane contact sites LD: lipid droplet, ER: endoplasmic reticulum, ATG: autophagy-related protein, LC3B: microtubule-associated protein 1 light chain 3 beta, LIR motif: LC3B-interacting region, P62: Sequestosome-1, PLIN1: Perilipin-1, ATGL: Adipose Triglyceride Lipase, VSP13: Variable small protein 13, DFCP1: Double FYVE-containing protein 1, PE: Phosphatidylethanolamine, ATP: Adenosine triphosphate, E2 enzyme: Ubiquitin-conjugating enzyme

Autophagy and lipid droplet (LD) biology are intricately linked, with autophagosomes breaking down LDs in response to different signals. LDs play crucial roles in forming autophagosomes possibly by providing essential lipids and serving as a supportive autophagosome assembly platform at the ER-LD interface. LDs and autophagosomes share common autophagic proteins, such as VPS13, Atg2, DFCP1, and Atg14, but their dual functions remain poorly understood. In our recent study, we found that prolonged starvation leads to Atg3 localizing to large LDs and lipidating LC3B, revealing a non-canonical autophagic role on LDs. In vitro, Atg3 associated with purified and artificial LDs, and conjugated Atg8 proteins. In long-term starved cells, only LC3B was found on the specific large LDs, positioned near LC3B-positive membranes that undergo lysosome-mediated acidification. This implies that LD-lipidated LC3B acts as a tethering factor, connecting autophagosomes to LDs and promoting degradation. Our data also support the notion that certain LD surfaces may function as lipidation stations for LC3B, which may move to nearby sites of autophagosome formation. Overall, our study unveils an unknown non-canonical implication of LDs in autophagy processes.

Autophagosomes are double-membrane vesicles that enclose and deliver micrometer-range cargos to lysosomes for degradation. A crucial reaction in the formation of autophagosomes is the conjugation of Atg8 proteins (LC3s or GABARAPs) to PE phospholipids of the membrane. This process is mediated by the E2 enzyme Atg3, which is recruited to the membrane by the ATG5-ATG12-ATG16L1 complex. Atg8 is the major ligand on autophagosomes that binds to cargo receptors having an LIR motif.

Autophagosomes degrade lipid droplets (LDs), the major cellular lipid storage organelles, through a process known as lipophagy. Lipophagy-targeted LDs are typically in the µm size range or smaller. Recently, we investigated the fate of larger LDs during starvation 1 . We utilized 3T3-L1 differentiated adipocytes and Huh7 liver cells, which can contain large LDs, e.g. >> 1µm in size. To allow sufficient time for LD consumption and explore various regimes of autophagy, the cells were starved for extended periods, up to 72hrs.

After approximately 6hrs of starvation, only LC3B and no other Atg8 proteins localized to a few large LDs. LD proteins like Plin1 and ATGL, which have LIR motifs, were not responsible for LC3B localization to the LDs. Similarly, P62, a major adaptor of LC3B that could bind to ubiquitinated LD proteins, was not responsible for the LD localization of LC3B. LC3B presence on LDs depended on Beclin1 and Atg5. Based on our overall observations, we developed two hypotheses: LC3B on LDs was from closely apposed autophagic membranes, or conjugation to the LD phospholipid monolayer. We ruled out the first possibility e.g. by employing a hypotonic medium-based organelle swelling approach, which demonstrated that the LC3B around LDs was not associated with a bilayer. Hence, the reliance of this phenomenon on Atg5/Beclin1 could be attributed to 1) the involvement of the proteins in the direct recruitment of conjugation machinery to LDs or 2) an interplay with classical autophagy to develop appropriate LD substrates. The latter model would involve the action of autophagy during the initial stages of starvation to release cues that trigger other autophagic mechanisms upon prolonged fasting. Hence, inhibiting the autophagic machinery would impede the release of these cues and prevent LC3B lipidation to LDs.

We investigated the possibility of LC3B being conjugated to LDs (Figure 1a). During prolonged starvation, we found that Atg3 was associated with specific larger LDs, supporting that it acted on LDs. When we deleted Atg3 or introduced a mutation (V15K) to reduce its binding efficiency, the localization of LC3B to LDs was hindered. Conversely, a mutant (K11W) that enhanced binding significantly increased LC3B localization to LDs. These results suggested that Atg3 may bind and conjugate LC3B to the LD phospholipids. To test this model, we adopted in vitro experiments using purified or artificial LDs. Atg3 bound and lipidated LC3B to these LDs; it was also capable of lipidating GABARAPLs. Lipidation always occurred, independent of artificial LD size, as long as PE was present. Hence, the preference for larger LDs was likely due to the interfacial phospholipid and proteome composition rather than the size of the LDs.

Atg3 demonstrated binding affinity to model LDs on its own, whereas it did not bind to giant unilamellar vesicles of similar lipid composition. This suggests that the physical recruitment of ATG3 to LDs in cells might be independent or differentially reliant on the Atg5-12-16 complex compared to bilayer membranes. The preference for the LD surface likely stems from the greater occurrence of phospholipid packing defects on the LD monolayer in comparison to the bilayer.

The LDs that were lipidated with LC3B were consistently in close proximity to LC3Bmembrane structures, likely autophagosomes, which appeared as more intense regions in Figure 1b. Through live imaging, we observed that LC3B-positive LDs were specifically recognized by LC3B-membrane structures. To demonstrate the affinity of this interaction, we conducted experiments using purified LC3B-positive vesicles, which were brought into contact with artificial LDs that either were conjugated or not to LC3B. Only LC3B-lipidated LDs developed contacts with the LC3B-vesicle.

The docked autophagosome-like structure became positive for lysosomal markers, indicating the degradation of compounds on or near the LDs, Figure 1c. Accordingly, the few LC3B-positive LDs often lacked PLIN1, suggesting that either PLIN1 was removed prior to LC3B lipidation or it was degraded following lipidation. Further investigations are needed to fully understand the fate of these LDs and PLIN1 in this process. Nonetheless, our findings support the idea that LC3B on the surface of LDs acts as a cofactor facilitating the recruitment of autophagosomes. Given that LC3B can undergo trans homodimerization, it is possible that LC3B-LC3B interactions mediated the interaction between LDs and autophagosomes.

Our findings also agree with the model that LD-lipidated LC3B is transferred to autophagosomes in close proximity, Figure 1c. This transfer could happen directly from the LDs to the autophagosomes via possible physical connections between LDs and autophagosomes, or via the ER. Notably, the trans-homodimerization of LC3B induces hemifusion of bilayers, meaning that, in principle, a contiguous LD-autophagosome leaflet could form. Such a pathway would be particularly activated by prolonged starvation, as a mechanism to alleviate the classical autophagosome biogenesis steps. Intriguingly here, Atg3 might associate with LDs alone and mediate LC3B lipidation, without the physical involvement of the upstream autophagic components. In this viewpoint, a constant supply of PE to the LDs would be necessary for LC3B continuous lipidation. PE supply for LD lipidation could occur through PE diffusion from the ER to LDs or transfer by lipid transfer proteins at ER-LD contact sites. Further investigations are required to explore these aspects and shed light on this possible pathway.

In summary, the localization of Atg3 and the subsequent lipidation of LC3B to LDs introduces a new perspective on the dynamic interplay between autophagy and LDs. By expanding the list of proteins common between autophagosome and LDs, VPS13, Atg2, DFCP1, Atg14, and now Atg3 and Atg8, we emphasize the intricate relationship between these organelles, which is important for their dynamics. It would be unsurprising to discover additional autophagy-related proteins localizing to LDs, depending on specific cues.

  1. Omrane M, M'Barek KB, Santinho A, Nguyen N, Nag S, Melia TJ, Thiam AR. LC3B is lipidated to large lipid droplets during prolonged starvation for noncanonical autophagy. Developmental Cell 2023; (FRM: EQU202103012564) and the Programme Emergence de la Ville de Paris (DDEEES 165). Figure legend: Figure 1 a) During prolonged starvation, ATG3 binds alone to LDs and lapidates LC3B, delivered by ATG7. Plin1 level on the LD is decreased on the LD that is lipidated. b) Confocal imaging of EGFP-LC3B (green) and LDs (LipidTox, magenta) in differentiated adipocytes (3T3-L1) expressing EGFP-LC3B and incubated in EBSS for 48 hours. The inset displays a 3D reconstructed image from z-stacks of the imaged LD. An LC3Bpositive membrane, possibly autophagosomes, is in close apposition to the LC3Blipidated c) Possible mechanisms. Lipidated LC3B on LD surfaces serves as a tethering factor for autophagosome recruitment. The recruited autophagosomes fuse with lysosomes and undergo acidification, without the LD being fully engulfment. LDs may function as a platform for LC3B lipidation, with lipidated LC3B potentially being transferred to recruited/nascent autophagosomes close to the LD surface, directly or via the ER membrane.

Acknowledgment. This work was supported by the Fondation pour la Recherche Médicale


H Chepfer

S Bony

D Winker

M Chiriaco

J L Dufresne

G Sèze

Use of CALIPSO lidar observations to evaluate the cloudiness simulated by a climate model

New spaceborne active sensors make it possible to observe the threedimensional structure of clouds. Here we use CALIPSO lidar observations together with a lidar simulator to evaluate the cloudiness simulated by a climate model: modeled atmospheric profiles are converted to an ensemble of subgridscale attenuated backscatter lidar signals from which a cloud fraction is derived. Except in regions of persistent thick upperlevel clouds, the cloud fraction diagnosed through this procedure is close to that actually predicted by the model. A fractional cloudiness is diagnosed consistently from CALIPSO data at a spatiotemporal resolution comparable to that of the model. The comparison of the model's cloudiness with CALIPSO data reveals discrepancies more pronounced than in previous model evaluations based on passive observations. This suggests that spaceborne lidar observations constitute a powerful tool for the evaluation of clouds in largescale models, including marine boundarylayer clouds.

Introduction

Clouds are the primary modulators of the Earth's radiation budget and still constitute the main source of uncertainty in model estimates of climate sensitivity [START_REF] Randall | Climate Models and Their Evaluation[END_REF]. Clouds simulated by climate models have long been evaluated using passive remotesensing [e.g. [START_REF] Zhang | Comparing Clouds And Their Seasonal Variations in 10 Atmospheric General Circulation Models With Satellite Measurements[END_REF], making the vertical structure of clouds difficult to assess. With the new generation of satellites carrying cloud profiling radar and lidar instruments such as CloudSat [START_REF] Stephens | The CloudSat mission and the Atrain[END_REF], ICESat/GLAS [START_REF] Spinhirne | Cloud and aerosol measurements from the GLAS: Overview and initial results[END_REF] and CALIOP/CALIPSO [START_REF] Winker | Initial Performance Assessment of CALIOP[END_REF]] a nearglobal view of the threedimensional distribution of clouds becomes possible. However, to make the comparison between modeled clouds and these observations meaningful, it is necessary to take into account the effects of cloud overlap, spatial resolution, and active signal attenuation. For this purpose, we propose a methodology of comparison of CALIOP/CALIPSO lidar data with modeled clouds that is intermediate between "modeltosatellite" [Morcrette 1991[START_REF] Doutriauxboucher | Simulation of Satellite Lidar and Radiometer Retrievals of a GCM ThreeDimensional Cloud Dataset[END_REF], Chiriaco et al. 2005, Chepfer et al. 2007] and "satellitetomodel" approaches: gridscale model outputs are converted to an ensemble of subgridscale lidar signals while highresolution CALIPSO lidar signals are averaged to the model vertical resolution , and then a cloud fraction is diagnosed from these signals similarly for model and for observations. The lidar simulator is described in Section 2, the processing of CALIPSO data is presented in Section 3 and the cloud detection method in Section 4. The methodology of modeldata comparison is then applied to the general circulation model (GCM) LMDZ4 [START_REF] Hourdin | The LMDZ general circulation model: climate performance and sensitivity to parameterized physics with emphasis on tropical convection[END_REF], and the results are presented in Section 5.

Simulation of lidar profiles from climate model outputs

Lidar equation

Our aim is to diagnose from climate model outputs the lidar profiles that would be observed by CALIPSO if the satellite was flying above an atmosphere similar to that predicted by the model. At the wavelength of the CALIOP/CALIPSO lidar (λ = 532 nm), atmospheric cloud particles and gas molecules contribute to scattering but not to absorption. In these conditions, the lidar ATtenuated Backscattered (ATB) signal corrected of geometrical effects and normalized to the molecular signal is given by :

∫ - z . . 2 TOA z dz (z)) + (z) e (z)) + (z) ( = ATB(z) mol sca, part sca, ( mol sca, part sca, α η (1)
where β sca,part , β sca,mol are lidar backscatter coefficients [m 1 sr 1 ] and α sca,part , α sca,mol attenuation coefficients [m 1 ] for particles and molecules, and η is a multiple scattering coefficient [Chiriaco et al. 2005]. The molecular properties (β sca,mol and α sca,mol ) are computed from [START_REF] Collis | Laser Monitoring of the Atmosphere[END_REF] as a function of temperature and pressure. The particle backscattered coefficient β sca,part is linked to the attenuation coefficient α sca,part and the backscattertoextinction ratio k part [sr 1 ] through :

(z) (z) k = (z) part sca, part part sca,
.

(2)

∫ z)dr Q(r)n(r, r π = (z) part sca, 2 (3)
and k part = P(π)/4π, where r is the particle radius, n(r,z) the particle size distribution, P(π) the backscattering phase function and Q(r) the scattering efficiency. To be consistent with the GCM, we assume that the cloud particles are spherical. Therefore, P(π) is parameterized as a function of the effective radius using Mie theory and, as most of the cloud particles are larger than λ = 532nm, Q(r) is set to 2. The multiple scattering coefficient (η) is variable in space and time as a function of the lidar footprint diameter, and the size, shape and density of the particles in the atmosphere. Theoretically it ranges between 0 and 1, but for CALIOP it is about 0.7 in ice clouds and less in water clouds [START_REF] Winker | Accounting for Multiple Scattering in Retrievals from Space Lidar[END_REF]]. Here we use η=0.7 and check the sensitivity of the results to this value in section 5.

To highlight the contribution of particles to the lidar signal, it is common to consider the scattering ratio SR(z), defined as the ratio of the total ATB signal (Equation 1) over the molecular ATB signal (which is the ATB signal in the absence of particules, i.e. when α sca,part (z)=0). SR(z) is equal to 1 in the absence of particles (gas molecules only) and, when the lidar signal is not attenuated, is greater than unity in the presence of particles.

Application to GCM outputs

The LMDZ4 GCM has an horizontal resolution of 2.5° in latitude, 3.75° in longitude, and is discretized into 19 vertical levels [START_REF] Hourdin | The LMDZ general circulation model: climate performance and sensitivity to parameterized physics with emphasis on tropical convection[END_REF]]. In each gridbox, the model calculates the vertical profiles of temperature, pressure, cloud fraction (assuming a maximumrandom overlap), cloud condensate and effective radius of cloud droplets and ice crystals. These profiles are converted to an ensemble of subgridscale profiles by dividing each gridbox into 50 subcolumns generated stochastically using the Subgrid Cloud Overlap Profile Sampler [START_REF] Klein | Validation and sensitivities of frontal clouds simulated by the ECMWF model[END_REF][START_REF] Webb | Combining ERBE and ISCCP data to assess clouds in three climate models[END_REF] : in each subcolumn, the cloud fraction is assigned to be 0 or 1 at every model level, with the constraint that the cloud condensate and cloud fraction averaged over all subcolumns is consistent with the gridaveraged model diagnostics and the cloud overlap model assumption.

In each subcolumn, the particle backscatter coefficient β sca,part and the attenuation coefficient α sca,part are the sums of liquid and ice cloud particles coefficients (aerosols are not considered in this study). They are computed from equations ( 2) and (3) after the model mixing ratios of incloud liquid and ice water content have been converted

to particle concentrations. This is done by assuming that the cloud particles are spherical with a radius equal to the effective radius predicted by the model. Finally, equation ( 1) is used to compute in each subcolumn of each gridbox the molecular and total ATB lidar profiles, and hence the SR profile.

CALIPSO observations

The CALIOP lidar, that operates at 532 and 1064 nm, was launched in April 2006 onboard the heliosynchro neous CALIPSO platform [START_REF] Winker | Initial Performance Assessment of CALIOP[END_REF]]. This study uses the total ATB measurements derived from the Version 1.10 Level 1 CALIOP dataset collected at 532 nm during night time from January to March 2007. The horizontal and vertical resolutions of these measurements are 333 m (333 m along track and 75 m cross track) and 30 m (respectively) below 8 km, and 1 km and 60 m at higher altitudes.

At each time and location, a molecular ATB profile is computed by using GMAO (Global Modeling and Assimil ation Office) temperature and pressure profiles [START_REF] Bey | Global modeling of tropospheric chemistry with assimilated meteorology: Model description and evaluation[END_REF]]. This profile is normalized (scaled) so as to match the Level 1 ATB profile measured between 27 and 29 km, where the atmosphere is free of particles, and then the scattering ratio (SR) is inferred by dividing the measured ATB profile by this molecular profile. Then, we reduce the vertical resolution of the CALIOP ATB profiles to that of the GCM (19 pressure levels) by vertical averaging. This makes the vertical resolution of CALIPSO profiles consistent with that of model compu tations and strongly increases the signaltonoise ratio of the lidar signals (section 4 and supplementary material).

Cloud diagnostics

The presence of cloud layers is then diagnosed consistently for the model and for observations from the vertical profile of lidar SR. For each CALIOP/CALIPSO profile or model subcolumn, we use different SR thresholds to label each atmospheric layer as "clear" (if 0.01 < SR ≤ 1.2), "cloudy" (if SR ≥ 3 of if the layer is the first fully attenuated layer encountered from top), "unclassified" (if 1.2 < SR < 3) or "fully attenuated" (if SR ≤ 0.01). The latest flag is set when optically thick clouds occurring in upper layers obscure the layers below and fully attenuate the laser beam. The conversion of SR in terms of optical depth depends on many factors. As an informal guesstimate, an SR value around 3 may be obtained with a liquid water cloud of 250500m depth, having a cloud optical depth of 0.015 0.03 and a liquid water path of 0.10.2 g/m 2 . The main cloud layers (which are associated with the largest SR values in the full resolution CALIOP profile) are detected by applying these thresholds to the vertically averaged CALIOP profile, but the optically thin cloud layers such as subvisible cirrus clouds and potentially some semitransparent clouds can be labelled as "unclassified" (Supplementary Material). This illustrates that the detection of a cloud layer obviously depends on the criteria used for this detection. Advanced algorithms of cloud fraction retrieval such as those used for level 2 CALIOP products would provide a better detection of the different cloud layers (especially of the optically thinnest and subvisible ones) and thus would likely report larger cloud fractions than our simple detection method. Although this is not critical for this modeldata comparison study because the cloud layers are diagnosed similarly in the model and in observations, this implies that our study ignores the optically thinnest cloud layers which may not be detected (in simulations like in observations) by our simple cloud detection criteria. Then, a "lidarderived cloud fraction" is computed for each vertical level at the resolution of the model gridbox as the ratio of the number of "cloudy" layers (or profiles) encountered within this gridbox over the total number of layers that are not "fully attenuated". Layered cloud fractions are also computed for three atmospheric layers : upper levels (between the 50 and 440 hPa pressure levels), middle levels (between 440 and 680 hPa) and low levels (altitudes below the 680 hPa level).

Results

5.1. Model cloud fractions diagnosed from the lidar simulator.

The lidar simulator described in section 2 is applied to daily atmospheric profiles derived from a LMDZ4 simulation forced by observed sea surface temperatures. The daily cloud fractions diagnosed from the simulator are then averaged in time to produce a monthlymean CALIPSOlike cloud fraction.

Sensitivity calculations indicate that the frequency of the model inputs and simulator calculations (once a day, every 3 hour, 1.5 hour or 30 min), as well as the value of the multiple scattering coefficient assumed in the simulator (e.g. using η = 0.3 instead of η = 0.7) have negligible impacts (less than 1 %) on the monthly mean CALIPSOlike cloud fraction (not shown). Using a different threshold (e.g. SR≥ 5 instead of SR≥ 3) does not change drastically the results neither, which suggests that the model does not simulate many optically thin clouds.

The layered cloud fractions derived from the lidar simulator and the simple cloud detection method are first compared to the cloud fractions originally predicted by the GCM. The CALIPSOlike cloud fractions (Figures 1b) are generally close to their original counterpart (Figures 1a),. Since the attenuation of the backscattered lidar signal at a given height depends on the presence of overlapping clouds at higher altitudes, such a result is expected at upper levels but not necessarily at lower levels. As shown by zonal means (Figures 2 ab), the original and CALIPSOlike cloud fractions at low and middlelevels differ mainly in regions of persistent convective activity (such as in the intertropical and southernpacific convergence zones and over the warm pool) and in extratropical regions associated with frontal cloudiness. However, over a large fraction of the Earth, the simulated upperlevel clouds are either not optically thick enough (the laser is typically attenuated when the optical depth is larger than 3) or not frequent or not persistent enough to obscure the view of low atmospheric layers from above. Despite the lidar attenuation effects, the CALIPSO simulator thus constitutes a powerful tool to diagnose and then to evaluate the different cloud types predicted by the GCM, including marine boundary layer clouds.

Before the era of spaceborne active remote sensing, GCM clouds were evaluated mostly against passive observations using the International Cloud Climatology Project (ISCCP) simulator [START_REF] Klein | Validation and sensitivities of frontal clouds simulated by the ECMWF model[END_REF][START_REF] Webb | Combining ERBE and ISCCP data to assess clouds in three climate models[END_REF][START_REF] Zhang | Comparing Clouds And Their Seasonal Variations in 10 Atmospheric General Circulation Models With Satellite Measurements[END_REF]. When applied to our GCM, the ISCCP and CALIPSO simulators both diagnose an upperlevel cloud fraction that resembles that originally predicted by the GCM (not shown). However, the lowlevel cloud fraction diagnosed from the CALIPSO simulator (Figure 1b) is much closer to that predicted by the GCM (Figure 1a) than that derived from the ISCCP simulator (Figures 1c).

Cloud fractions derived from CALIPSO observations.

The layered and zonal mean vertical distributions of the cloud cover diagnosed from CALIPSO data using the methodology discussed in sections 3 and 4 are shown in Figures 3cd,1d, and 2c respectively. As expected for the JanuaryFebruaryMarch season, the upper level cloud fraction exhibits maxima in deep convective regions located over tropical land areas, over the maritime continent and the pacific warm pool and over the southern pacific convergence zone. Secondary maxima occur at extratropical latitudes. Middlelevel clouds occur in middlelatitudes and in regions of deep convection. CALIPSO data suggest that lowlevel clouds occur nearly everywhere over ocean. Maximum lowlevel cloud fractions are found at subtropical latitudes, at the eastern side of the ocean basins which are known to be covered by stratus and stratocumulus clouds, but cloud fractions of 40 to 60 % are also reported over the tropical trade winds and at middle latitudes. Note that owing to the pronounced diurnal cycle of oceanic tropical low clouds [START_REF] Rozendaal | An observational study of diurnal variations of marine stratiform cloud[END_REF], these nightime values might be higher than daytime or diurnally averaged values.

Sensitivity tests using the CALIPSO lidar color ratio to detect the presence of aerosols (defined as cases where the ratio between the ATB measured at 1064 and 532 nm exceeds 0.6, [START_REF] Liu | CALIOP Algorithm Theoretical Basis Document, Part 3: Scene Classification Algorithms[END_REF]) show that the cloud detection criteria defined in section 4 do not misinterpret the presence of aerosols as lowlevel clouds : atmospheric layers with large aerosol burden, such as over Africa or over ocean regions west of the Sahara, are labelled as "unclassified". On the other hand, other tests show that the cloud fraction reported in the trades critically depends on the SR threshold used for cloud detection (e.g. it is substantially reduced when the threshold is changed from 3 to 5), while it is not the case at the eastern side of the ocean basins or in the extratropics. This suggests that the large cloud fractions reported in the trades are related to the ability of lidar measurements to detect optically thin or small broken clouds owing both to their high horizontal resolution and to their high sensitivity to the presence of cloud particles., This is corroborated by the cloud fractions derived from ICESat/GLAS lidar measurements (level 2 dataset, GLA09, release 26, [START_REF] Palm | GLAS atmospheric data products. GLAS algorithms theoretical basis document[END_REF] on October and November 2003 which are of same order of magnitude as those derived from CALIPSO at the same season (not shown). The difficulty of detecting thin or/and broken lowlevel clouds with passive remote sensing [START_REF] Turner | Thin liquid water clouds: their importance and our challenge[END_REF] presumably explains why data from the International Cloud Climatology Project [ISCCP, [START_REF] Rossow | Advances in understanding clouds from the ISCCP[END_REF] report substantially smaller lowlevel cloud fractions (Figure 1e).

Evaluation of the GCM cloudiness against CALIPSO observations.

The comparison of the cloud fraction derived from the lidar simulator with that derived from CALIPSO observations reveals large biases in the cloudiness simulated by the GCM (Figures 1, 2 and 3). The upperlevel cloud fraction is underestimated by several tens of % in the tropical deep convective regions and the middlelevel cloudiness is systematically underestimated, especially in the intertropical convergence zone (ITCZ) and in middlelatitude regions associated with frontal clouds. But the largest biases appear for the lowlevel cloud fraction, which is strongly overestimated (by about 25 %) over midlatitude oceans and is strongly underestimated (by up to a factor of 4) in the subtropics, both over the cold waters at the eastern side of the Pacific, Atlantic and Indian oceans, and offshore over the subtropical oceans. Except in midlatitudes and in the ITCZ, these biases are unlikely to result from the attenuation of the lidar signal by overlapping optically thick clouds (section 5.1). Both in the tropics and in middle latitudes, the modeldata discrepancies revealed by this comparison are more pronounced than in previous comparisons using ISCCP observations and the ISCCP simulator (Figure 1) This suggests that using spaceborne lidar observations to evaluate largescale models will help to point out systematic biases in the cloudiness predicted by the models, especially in the simulation of marine boundarylayer clouds.

Conclusion

This paper proposes a methodology to evaluate the cloud fraction simulated by a GCM using space borne lidar observations : on the one hand, the atmospheric profiles predicted by the GCM are converted (through a lidar simulator) to an ensemble of subgridscale lidar signals similar to those observed from space; on the other hand, the high vertical resolution of the lidar profiles measured by CALIOP/CALIPSO is degraded down to the vertical resolution of the GCM; then, a fractional cloudiness is diagnosed from the simulated or observed profiles of lidar scattering ratio using simple and consistent criteria. When applied to the LMDZ4 GCM, the CALIPSOlike cloud fraction derived from this methodology reproduces well the cloudiness originally predicted by the model, and this despite the attenuation effect of the lidar signal and whatever the assumed value of the multiple scattering coefficient. When applied to CALIPSO data, this simple methodology efficiently detects the main cloud layers, does not misinterpret them with the presence of aerosols, and provides cloud fraction estimates that are of same order of magnitude than those derived from the IceSAT/GLAS level 2 dataset. The comparison of the cloud fraction derived from the lidar simulator with that derived from CALIPSO data reveals large biases in the GCM cloudiness, which are more dramatic than those pointed out previously using passive observations. In particular, striking differences between GCM and CALIPSO observations (found for all seasons, not shown) occur for the marine boundarylayer clouds observed over most of the subtropical oceans.

As the response of these clouds to climate change constitutes a key uncertainty for GCM estimates of climate sensitivity [START_REF] Bony | Marine boundary layer clouds at the heart of tropical cloud feedback uncertainties in climate models[END_REF], it is of great importance to improve their representation in climate models. Efforts to improve the parameterization of shallow clouds in LMDZ4 are under way [START_REF] Rio | A Thermal Plume Model for the Convective Boundary Layer: Representation of Cumulus Clouds[END_REF], and the methodology presented in this paper will allow us to quantify the extent to which the simulated cloud fraction is improved at the large scale. This study thus shows the large potential of CALIPSO observations and of the lidar simulator for the evaluation of clouds simulated by climate models. Within the framework of the Cloud Feedback Model Intercomparison Project (CFMIP, http://www.cfmip.net), the lidar simulator will be applied to several other climate models to examine intermodel differences and systematic biases in the simulation of clouds by GCMs. For this purpose, the GCM oriented CALIOP dataset processed and used in this study will be made available to the modelling community. A comparison between this dataset and the CALIOP Level 2 product is underway (both for nightime and daytime data) and will be reported in a future paper.

Figure 1 :

 1 Figure 1: Lowlevel cloud fractions averaged for the JanuaryFebruaryMarch season. (a) The cloud fractions actually predicted by the GCM, (b) the cloud fractions derived from the lidar simulator and (d) from CALIOP/CALIPSO observations using similar consistent criteria. Also reported at the bottom is (c) the lowlevel cloud fraction derived from ISCCPD2 data and (e) that diagnosed from the ISCCP simulator.

Figure 2 :

 2 Figure 2: Vertical distribution of the zonally averaged cloud fraction for JanuaryFebruaryMarch: a) Original cloud fraction predicted by the GCM; b) GCM cloud fraction diagnosed from the lidar simulator; c) Cloud Fraction derived from CALIOP/CALIPSO data.

Figure 3 :Figure 1 :

 31 Figure 3: (Top) Upperlevel and (bottom) middlelevel layered cloud fractions averaged for the January FebruaryMarch season. The cloud fractions derived from the lidar simulator (left panels) and from CALIOP/CALIPSO observations (right panels) are diagnosed consistently by using similar criteria.

Acknowledgments.

The authors thank Vincent Noel for fruitful discussions and for his help in the processing of CALIOP level 1 data, and to Ionela Musat for her help in GCM data processing. NASA, CNES, ICARE and Climserv are acknowledged for giving access to the CALIOP/CALIPSO data. Constructive comments from two anonymous reviewers greatly helped to improve the manuscript.