Kristyna Brejchova

Veronika Paluchova

Marie Brezinova

Tomas Cajka

Laurence Balas

Thierry Durand

Marcela Krizova

Zbynek Stranak

Ondrej Kuda
email: ondrej.kuda@fgu.cas.cz

Triacylglycerols containing

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Nourishment plays a key role in the early stage of life. Breast milk provides the ideal cocktail of nutrients for infants in this period of rapid development and growth and breastfeeding have a long-term impact on health [START_REF] Floris | Human milk fatty acid profile across lactational stages after term and preterm delivery: A pooled data analysis[END_REF]. Although total contents of milk macronutrients are relatively stable, the levels of distinct components can be altered by several factors such as maternal dietary habits, mode of delivery, time of delivery, and/or stage of lactation [START_REF] Ballard | Human Milk Composition: Nutrients and Bioactive Factors[END_REF][START_REF] Koletzko | Human Milk Lipids[END_REF][START_REF] Hobbs | The impact of caesarean section on breastfeeding initiation, duration and difficulties in the first four months postpartum[END_REF][START_REF] Lee | Compositional Dynamics of the Milk Fat Globule and Its Role in Infant Development[END_REF].

Milk-producing mammary gland is a dynamic organ composed by the ductal tree system, which is embedded in mammary fat pad. During pregnancies, repeatable lobulo-alveolar expansions and adipocyte re-differentiations occur to support sufficient milk production during lactation periods [START_REF] Watson | Mammary development in the embryo and adult: new insights into the journey of morphogenesis and commitment[END_REF][START_REF] Wang | Reversible Dedifferentiation of Mature White Adipocytes into Preadipocyte-like Precursors during Lactation[END_REF][START_REF] Morroni | Reversible transdifferentiation of secretory epithelial cells into adipocytes in the mammary gland[END_REF]. The first fluid produced after delivery -colostrum -is rich in proteins, immunologic components, and developmental factors and it plays a primary role in immunomodulation and promotion of intestinal development [START_REF] Ballard | Human Milk Composition: Nutrients and Bioactive Factors[END_REF]. In few days postpartum (3-6 days), transitional milk replaces colostrum and fully mature milk is produced after two weeks post-delivery.

Human milk lipids provide a major energy supply for infants and the main lipid species are triacylglycerols (TAGs) which form the core of the milk fat globule covered by phospholipid tri-layer membrane [START_REF] Lee | Compositional Dynamics of the Milk Fat Globule and Its Role in Infant Development[END_REF]. The membrane comprises predominantly amphipathic lipids that originates from mammary gland epithelia (glycerophospholipids, sphingolipids, cholesterol) [START_REF] Koletzko | Human Milk Lipids[END_REF][START_REF] Lee | Compositional Dynamics of the Milk Fat Globule and Its Role in Infant Development[END_REF]. Besides the structural and nutritional lipids, lipid-based bioactive components such as specialized proresolving mediators and fatty acid esters of hydroxy fatty acids (FAHFA) were identified in the milk [START_REF] Arnardottir | Human milk proresolving mediators stimulate resolution of acute inflammation[END_REF][START_REF] Brezinova | Levels of palmitic acid ester of hydroxystearic acid (PAHSA) are reduced in the breast milk of obese mothers[END_REF].

FAHFAs belong to a family of branched lipids, some of which have anti-inflammatory and antidiabetic properties [START_REF] Brejchova | Understanding FAHFAs: From structure to metabolic regulation[END_REF][START_REF] Aryal | Distinct Biological Activities of Isomers from Different Families of Branched Fatty Acid Esters of Hydroxy Fatty Acids (FAHFAs)[END_REF]. One specific FAHFA family called palmitic acid ester of hydroxystearic acid (PAHSAs) could stimulate glucagon-like peptide 1 (GLP-1) secretion in the gut and improve systemic glucose homeostasis [START_REF] Syed | Palmitic Acid Hydroxystearic Acids Activate GPR40, Which Is Involved in Their Beneficial Effects on Glucose Homeostasis[END_REF]. Oral administration of 5-PAHSA isomer prevented gut mucosal damage and colitis in mice [START_REF] Lee | Branched Fatty Acid Esters of Hydroxy Fatty Acids (FAHFAs) Protect against Colitis by Regulating Gut Innate and Adaptive Immune Responses[END_REF]. Therefore, PAHSAs from food can modulate metabolic homeostasis and immune response in the intestine [START_REF] Brejchova | Understanding FAHFAs: From structure to metabolic regulation[END_REF][START_REF] Riecan | Branched and linear fatty acid esters of hydroxy fatty acids (FAHFA) relevant to human health[END_REF].

FAHFAs are either de novo synthesized in cells capable of (de novo) lipogenesis or supplied by dietary sources [START_REF] Brezinova | Levels of palmitic acid ester of hydroxystearic acid (PAHSA) are reduced in the breast milk of obese mothers[END_REF][START_REF] Paluchova | Lipokine 5-PAHSA Is Regulated by Adipose Triglyceride Lipase and Primes Adipocytes for De Novo Lipogenesis in Mice[END_REF][START_REF] Kolar | Linoleic acid esters of hydroxy linoleic acids are anti-inflammatory lipids found in plants and mammals[END_REF][START_REF] Kuda | Nrf2-Mediated Antioxidant Defense and Peroxiredoxin 6 Are Linked to Biosynthesis of Palmitic Acid Ester of 9-Hydroxystearic Acid[END_REF][START_REF] Yore | Discovery of a class of endogenous mammalian lipids with anti-diabetic and anti-inflammatory effects[END_REF]. Free FAHFAs can be esterified to glycerol giving rise to FAHFA acylglycerols, which in combination with other two fatty acids (FAs) result in the formation of TAG estolides (TAG EST) [START_REF] Paluchova | Lipokine 5-PAHSA Is Regulated by Adipose Triglyceride Lipase and Primes Adipocytes for De Novo Lipogenesis in Mice[END_REF][START_REF] Brejchova | Distinct roles of adipose triglyceride lipase and hormone-sensitive lipase in the catabolism of triacylglycerol estolides[END_REF][START_REF] Brezinova | Exercise training induces insulin-sensitizing PAHSAs in adipose tissue of elderly women[END_REF][START_REF] Tan | Discovery of FAHFA-Containing Triacylglycerols and Their Metabolic Regulation[END_REF]. FAHFAs can be released from TAG estolide storage pool by adipose triglyceride lipase (ATGL), which is also free FAHFA hydrolase and the enzyme responsible for acylremodelling reaction among acylglycerols [START_REF] Brejchova | Distinct roles of adipose triglyceride lipase and hormone-sensitive lipase in the catabolism of triacylglycerol estolides[END_REF][START_REF] Tan | Discovery of FAHFA-Containing Triacylglycerols and Their Metabolic Regulation[END_REF]. Free FAHFAs' estolide bond can be hydrolyzed by carboxyl ester lipase (CEL) [START_REF] Kolar | Branched Fatty Acid Esters of Hydroxy Fatty Acids Are Preferred Substrates of the MODY8 Protein Carboxyl Ester Lipase[END_REF], androgen-induced gene 1, androgen-dependent TFPI-regulating protein [START_REF] Ertunc | AIG1 and ADTRP are endogenous hydrolases of fatty acid esters of hydroxy fatty acids (FAHFAs) in mice[END_REF][START_REF] Parsons | AIG1 and ADTRP are atypical integral membrane hydrolases that degrade bioactive FAHFAs[END_REF], and hormone sensitive lipase [START_REF] Brejchova | Distinct roles of adipose triglyceride lipase and hormone-sensitive lipase in the catabolism of triacylglycerol estolides[END_REF]. Importantly, CEL belongs to abundant proteins in human milk and it contributes to TAG digestion in the immature gut of infants [START_REF] Olivecrona | Human milk lipases and their possible role in fat digestion[END_REF][START_REF] Miller | Carboxyl ester lipase from either mother's milk or the pancreas is required for efficient dietary triglyceride digestion in suckling mice[END_REF].

We have previously identified 5-PAHSA as a signalling lipid, which was reduced in the colostrum of obese mothers. In the present study we hypothesized that 1) timing, mode of delivery, and stage of lactation has an impact on levels of bioactive FAHFAs; and 2) TAG estolides are present in human milk and could be released by the CEL in milk.

Methods

Chemicals

LC-MS-grade solvents and mobile phase modifiers were obtained from VWR International, (Czech Republic). FAHFA standards were from Cayman Europe (Tallinn, Estonia), and 5-PAHSA was synthesized as previously described [START_REF] Balas | Regiocontrolled syntheses of FAHFAs and LC-MS/MS differentiation of regioisomers[END_REF]. The TAG estolides were synthesized at the Institute of Biomolecules Max Mousseron (Montpellier) [START_REF] Brejchova | Distinct roles of adipose triglyceride lipase and hormone-sensitive lipase in the catabolism of triacylglycerol estolides[END_REF]. All other chemicals were obtained from Merck (Czech Republic) unless stated otherwise.

Sample collection

Milk samples from healthy lactating mothers who delivered spontaneously in term (VB, vaginal birth), or by caesarean section (CS), or preterm (PB, preterm birth) were collected at the Institute for Mother and Child Care in Prague, Czech Republic. Mothers, who delivered spontaneously in term or underwent caesarean section, gave birth between 38-41 weeks of gestation (Table 1). Their newborns had a chronological age 72 hours and 28 days at the time of milk collection, respectively. Preterm infants were delivered before 32nd week of gestation (Table 1, Figure 1) and samples of milk were collected after 72 hours and 28 days. Moreover, another milk sample was acquired at the time when the infants would have reached the 36th week of gestation. We considered the preterm delivery condition more relevant for data interpretation than the mode of preterm delivery (Table 1). Therefore, the PB group contains both vaginal birth and caesarean section cases. Milk samples were frozen in -20 °C within 30 minutes after collection and transported to -80 °C freezer in 1-3 days. All donors signed informed consent before inclusion in the study. The study protocols were approved by the Research Ethics Board of the Institute (Approval number: 2016-12-19/ 3).

Sample extraction

Before extraction, frozen milk samples were warmed up to 37 °C in a water bath and homogenized within 3 minutes using an ultrasonic bath to get samples of consistent quality. Total lipids were extracted using biphasic solvent system of cold methanol, methyl tert-butyl ether (MTBE), and water according to published method [START_REF] Brezinova | Levels of palmitic acid ester of hydroxystearic acid (PAHSA) are reduced in the breast milk of obese mothers[END_REF][START_REF] Brezinova | Exercise training induces insulin-sensitizing PAHSAs in adipose tissue of elderly women[END_REF]. In brief, 300 µl of milk was mixed with pre-chilled 660 µl of methanol, 2,400 µl of MTBE and 10% methanol containing internal standards. Samples were centrifuged at 5,200 x g, 10 minutes at 4 °C and aliquots of the organic phase for lipidomics (100 µL) and TAG estolide/FAHFA (2 mL) analysis were collected. Samples of water phase (70 µL) were used for metabolomics and protein pellets for proteomic analysis. TAG estolide and FAHFA extraction was performed using SPE columns (HyperSep Silica 500 mg/10 mL, 40-60 µm, 60 Å, ThermoFisher Scientific, Czech Republic) similarly as was described before [START_REF] Brezinova | Levels of palmitic acid ester of hydroxystearic acid (PAHSA) are reduced in the breast milk of obese mothers[END_REF].

CEL activity assay

To determine CEL hydrolase affinity to TAG estolide substrates the previously published CEL activity assay method was modified and performed [START_REF] Kolar | Branched Fatty Acid Esters of Hydroxy Fatty Acids Are Preferred Substrates of the MODY8 Protein Carboxyl Ester Lipase[END_REF][START_REF] Xiao | Pancreatic lipase-related protein-2 (PLRP2) can contribute to dietary fat digestion in human newborns[END_REF]. Positional isomers of TAG EST 70:4 containing either 5-or 9-PAHSA esterified at sn-1/sn-3 (e.g. 5-PAHSA/18:2/18:2) or sn-2 positions of glycerol backbone (e.g. 18:2/5-PAHSA/18:2) as well as trilinolein (18:2/18:2/18:2) were used as CEL substrates. In brief, 3.8 µg of purified CEL (BioVendor R&D, Czech Republic) was used per reaction in the assay buffer containing 150 mM NaCl, 1 mM Tris-HCl, 2 mM CaCl 2 , pH 8. TAG estolide substrate (250 µM) was emulsified in the assay buffer with 10 mM sodium cholate. To prevent product inhibition, fatty acid-free bovine serum albumin (MP Biomedicals, France) was added to the reaction at a final concentration of 1%. Volume 100 μL of the enzyme solution and 100 μL of the substrate were mixed and incubated at 37 °C for 60 min. The reaction was terminated in ice-cold water bath. MTBE extraction was performed as before [START_REF] Brejchova | Distinct roles of adipose triglyceride lipase and hormone-sensitive lipase in the catabolism of triacylglycerol estolides[END_REF]. Control reactions were performed under identical conditions without the enzyme. Reaction products were analyzed using lipidomics platforms.

Liquid chromatography and mass spectrometry (LC-MS) analyses.

The LC-MS systems consisted of a Vanquish UHPLC System (Thermo Fisher Scientific, Bremen, Germany) coupled to a Q Exactive Plus mass spectrometer (Thermo Fisher Scientific, Bremen, Germany) and UltiMate 3000 RSLC UHPLC system coupled to a QTRAP 5500/SelexION mass spectrometer (SCIEX, Darmstadt, Germany). LIMeX (LIpids, Metabolites, and eXposome compounds) LC-MS workflow was used for the untargeted analysis of complex lipids and polar metabolites as before [START_REF] Lopes | Metabolomics atlas of oral 13C-glucose tolerance test in mice[END_REF]. Targeted milk FAHFA analysis was performed as before [START_REF] Brezinova | Levels of palmitic acid ester of hydroxystearic acid (PAHSA) are reduced in the breast milk of obese mothers[END_REF].

Milk metabolome and lipidome dataset analysis

Metabolomics data were analyzed using statistical enrichment approach based on chemical similarity [START_REF] Barupal | Chemical Similarity Enrichment Analysis (ChemRICH) as alternative to biochemical pathway mapping for metabolomic datasets[END_REF]. Customized ChemRICH software was used to generate metabolite clusters annotated according to Medical Subject Headings classes. Kolmogorov-Smirnov test was used to calculate adjusted cluster p value and z-score for each cluster. The sign of z-score defines whether the clustered metabolites are downregulated (z-score < 0) or upregulated (z-score > 0). Clusters with |z-score| > 2 and -log 10 (adjusted cluster p value) > 1 were considered significantly altered. ChemRICH uses classic 881 bit molecular fingerprints defined by PubChem which is not optimal for lipid-rich datasets. Therefore, the algorithm generated one big cluster of TAGs (n = 140). To overcome this limitation, we performed the enrichment analysis of the TAG cluster using 1024 bit 'hybridization' type of fingerprint which explored 'depth' of 7 atoms and considered the presence of double bonds [START_REF] Guha | Chemical Informatics Functionality in R[END_REF]. The final 8 TAG sub-clusters were annotated using Lipid Mini-On and the TAG double bond number was considered as the significant sub-cluster discriminator (q values < 0.05) [START_REF] Clair | Lipid Mini-On: mining and ontology tool for enrichment analysis of lipidomic data[END_REF].

Statistics

GraphPad Prism 9.2.0 statistical tools were used for bar graphs. [START_REF] Faul | G*Power 3: a flexible statistical power analysis program for the social, behavioral, and biomedical sciences[END_REF].

Results

Impact of lactation stage, gestational age, and mode of delivery on metabolite composition of breast milk

We analysed metabolite composition of breast milk by LC-MS in the three groups of mothers at different time points of lactation and selected VB samples at 72h after childbirth as a reference group. Annotated data were sorted into 55 clusters using chemical similarity approach and up-or down-regulated clusters visualized (Figure 1A, Table S1). Enrichment analysis showed that CS colostrum had significantly lower levels of TAGs and higher levels of several phosphatidylinositol (PI) and phosphatidylserine (PS) classes than VB. In contrast, PB colostrum was highly enriched in several phospholipid and sphingomyelin classes as compared to VB group. (Figure 1B, Table S2).

The largest cluster of TAGs was further sub-clustered based on acyl chain composition to gain a better insight into the changes of milk lipid composition in time (Figure 2). Heatmap representation of the clusters showed that metabolic profiles in all three groups followed similar patternsdecrease of phospholipids and increase of unsaturated acylglycerols. Composition of PB mature milk caught up with the VB and CS groups after 28 days and remained stable until gestational week 36.

We focused only on FAHFA-related lipids further on. 1.

Milk FAHFA isomer profiles change in time

We have previously shown that human colostrum contains several FAHFA species and that 5-PAHSA levels were lower in obese compared to lean mothers [START_REF] Brezinova | Levels of palmitic acid ester of hydroxystearic acid (PAHSA) are reduced in the breast milk of obese mothers[END_REF]. Here we explored also other FAHFA with documented biological activities [START_REF] Aryal | Distinct Biological Activities of Isomers from Different Families of Branched Fatty Acid Esters of Hydroxy Fatty Acids (FAHFAs)[END_REF][START_REF] Kolar | Linoleic acid esters of hydroxy linoleic acids are anti-inflammatory lipids found in plants and mammals[END_REF][START_REF] Yore | Discovery of a class of endogenous mammalian lipids with anti-diabetic and anti-inflammatory effects[END_REF][START_REF] Benlebna | Long-term high intake of 9-PAHPA or 9-OAHPA increases basal metabolism and insulin sensitivity but disrupts liver homeostasis in healthy mice[END_REF].

First, we tested whether the term or mode of delivery affects FAHFA concentration. We detected significantly lower levels of both 5and 9-PAHSA in colostrum samples of mothers who delivered preterm compared to those who gave birth spontaneously in term. Significantly lower concentration of 5-PAHSA was also found in milk samples of mothers in CS group (Figure 3A). The levels of both PAHSAs normalized in mature milk (28 days). The only difference was detected in PB group, in which the concentration of 5-PAHSA was still significantly reduced in comparison with the VB group. None of the other bioactive FAHFAs in colostrum and mature milk revealed significant difference among groups of delivery (Figure 3B). 1).

Second, we looked into the change of FAHFAs levels during different stages of lactation. In the VB group, the only significant change observed was a decrease in 5-PAHSA and a rise in 13-LAHLA levels in mature milk compared to colostrum (Figure 3B). Significantly increased concentrations of seven FAHFAs (5-PAHSA, 10-PAHSA, 9-PAHPA, 9-, 10-, 12/13-OAHSA and 13-LAHLA) and three FAHFAs (5-PAHSA, 10-and 12/13-OAHSA) were found in mature milk from CS and PB groups, respectively (Figure 3B).

TAG estolides in the milk

Milk fat globule TAG core is synthesized from cytoplasmic lipid droplets by mammary gland alveolar epithelial cells. Metabolism of FAHFAs and their storage form, TAG estolides, is closely associated with ATGL that acts on lipid droplet surface. Therefore, we tested whether the milk contains TAG estolides present in human adipose tissue [START_REF] Brezinova | Exercise training induces insulin-sensitizing PAHSAs in adipose tissue of elderly women[END_REF]. We confirmed the presence of several TAG estolides in the milk: TAG EST 68:1, 68:2, 68:3, 68:4, 70:2, 70:3, 70:4, and 70:5. However, their levels were very variable and low or undetectable in some samples and no significant difference was found among the groups. The most abundant TAG estolides represented a mixture of isomers containing PAHSAs, palmitic acid esters of hydroxylinoleic acid, oleic acid esters of hydroxystearic acids, and oleic acid esters of hydroxyoleic acids in agreement with free FAHFA species (Figure 4). 1).

CEL hydrolyzes glycerol ester bonds and FAHFA estolide bond in TAG estolides

We hypothesized that the TAG estolides in the milk could be degraded by the CEL, because PAHSAs are known CEL substrates [START_REF] Kolar | Branched Fatty Acid Esters of Hydroxy Fatty Acids Are Preferred Substrates of the MODY8 Protein Carboxyl Ester Lipase[END_REF]. We tested the affinity of purified CEL to TAG estolides with a small combinatorial library of one PAHSA and two linoleic acids (LA) bound to the glycerol backbone and compared it with a common CEL substrate, glyceryl trilinoleate. The assay was focused on three main TAG estolide structural features (Figure 5). First, to evaluate substrate preference to different TAG estolide/FAHFA regioisomers, we used TAG estolides containing either 5-PAHSA or 9-PAHSA. Second, to evaluate sn-position preference, we used two positional isomers, in which PAHSA was esterified either at the primary (sn-1/3) alcohol position or at secondary alcohol position (sn-2) of glycerol. Third, to test whether CEL hydrolyzes FAHFA estolide bond when bound to TAG estolide, we monitored hydroxylated acylglycerol intermediates (Figure 5). CEL completely degraded all substrates giving rise to several intermediates and products of the cleavage (Figure 5). CEL can cleave LA-glycerol bond, release LA and produce an intermediatediacylglycerol estolide (DAG EST) (Figure 5, green cloud). This could also serve as a CEL substrate generating another molecule of LA and a monoacylglycerol estolide (MAG EST). Our results showed that CEL effectively released LA from both TAG EST regioisomers and their positional isomers as well as from trilinolein. Corresponding reaction products such as DAG EST and MAG EST were also found in all TAG EST substrate samples (Figure 6A).

CEL can also cleave the PAHSA-glycerol bond resulting in the production of a PAHSA and a DAG (Figure 5, orange cloud). Presence of free PAHSA found in all reactions confirms the ability of the lipase to release 5-and 9-PAHSA from both TAG estolide positional isomers, although CEL was more effective when PAHSA was bound at sn-2 position. This preference is evident not only from significantly higher PAHSAs, but also from lower levels of MAG ESTs in samples of these positional isomers (Figure 6A,B). Moreover, considerably lower values of MAG EST 9P than MAG EST 5P of both positional isomers (p < 0.0001 for both positional isomers, compare absolute values between panels) indicate higher affinity of CEL to 9-PAHSA-than 5-PAHSA-glycerol ester bond (Figure 6A).

The reaction scheme in the blue cloud (Figure 5, blue cloud) depicts possible cleavage of PAHSA estolide bond within TAG estolide structure resulting in the release of free palmitic acid (PA) and hydroxylated TAG. LA could be subsequently hydrolyzed from this intermediate and hydroxylated DAGs and MAGs may arise (Figure 5, blue cloud). However, CEL could also react with a pool of intermediate products (DAG EST and MAG EST) in the same way (Figure 5, green cloud) and generate hydroxylated DAGs and MAGs. Presence of all of these products in samples with TAG estolide containing 5-PAHSA provided an evidence that CEL has such ability and prefers TAG estolide with 5-PAHSA bound at the primary (sn-1,3) alcohol position in the glycerol (Figure 6C).

Since hydrolysis of 9-PAHSA TAG estolide substrates gave rise to negligible levels of hydroxylated DAG, it is plausible that CEL rather cleaves 9-PAHSA-glycerol ester bond than FAHFA estolide bond (Figure 6C). It has been shown that free 9-PAHSA is more effectively hydrolyzed by CEL than 5-PAHSA [START_REF] Kolar | Branched Fatty Acid Esters of Hydroxy Fatty Acids Are Preferred Substrates of the MODY8 Protein Carboxyl Ester Lipase[END_REF]. This rule also apply to TAG estolide substrates, as significantly higher levels of the final products (9-HSA and PA, p < 0.0001), as well as lower levels of free 9-PAHSA (p < 0.0001) were found compared to reactions with TAG EST containing 5-PAHSA substrate (Figure 6B). Regioisomer comparison between panels is described in the text for clarity. No suitable internal standards were available for quantitation of hydroxylated acylglycerols.

Discussion

Our data suggest that milk from VB group has the optimal balanced combination of free FAHFAs and TAG estolides. Premature birth and/or caesarean section represent a situation in which the mammary gland is not ideally stimulated or ready for milk production, but the milk lipid composition normalizes in time.

The caesarean section rate is expected to increase above the medically necessary threshold in future neglecting the potential negative effects for the mother and the child [START_REF] Betran | Trends and projections of caesarean section rates: global and regional estimates[END_REF]. The caesarean section has a negative impact on breastfeeding initiation, duration, and difficulties (4, 37). Prior et al. systematically reviewed association between CS and breastfeeding and concluded that rates of early breastfeeding initiation were lower in CS compared with VB group [START_REF] Prior | Breastfeeding after cesarean delivery: a systematic review and meta-analysis of world literature[END_REF]. Delayed initiation of milk production could affect lipid composition of the colostrum acquired at 72 h postpartum (Figure 1A), resulting in higher phospholipid levels and lower TAG levels in CS group. Higher levels of phospholipid and sphingomyelin classes, that form the milk fat globule tri-layer, were obvious also in PB colostrum in agreement with previous observations [START_REF] Pérez-Gálvez | Effect of gestational age (preterm or full term) on lipid composition of the milk fat globule and its membrane in human colostrum[END_REF]. Differences in TAG and phospholipid levels or in their ratio could be linked to milk fat globule size and potential digestibility [START_REF] Argov | Size-Dependent Lipid Content in Human Milk Fat Globules[END_REF][START_REF] Garcia | The size and interfacial composition of milk fat globules are key factors controlling triglycerides bioavailability in simulated human gastro-duodenal digestion[END_REF]. In comparison to colostrum samples, mature milk metabolomic and lipidomic profiles were affected by multiple factors (diet, sample acquisition, etc.), but the differences vanished with advancing lactation.

The logical source of both free FAHFAs and TAG estolides in the milk is the breast white adipose tissue that transdifferentiated into milk-secreting glands during pregnancy [START_REF] Morroni | Reversible transdifferentiation of secretory epithelial cells into adipocytes in the mammary gland[END_REF]. FAHFA levels produced into the milk are comparable to serum/plasma FAHFA concentrations [START_REF] Riecan | Branched and linear fatty acid esters of hydroxy fatty acids (FAHFA) relevant to human health[END_REF] and the structure of milk fat globules [START_REF] Lee | Compositional Dynamics of the Milk Fat Globule and Its Role in Infant Development[END_REF] suggests that the TAG lipid core could store TAG estolides.

FAHFA production is dependent on de novo lipogenesis and/or dietary lipids [START_REF] Brejchova | Understanding FAHFAs: From structure to metabolic regulation[END_REF][START_REF] Yore | Discovery of a class of endogenous mammalian lipids with anti-diabetic and anti-inflammatory effects[END_REF]. However, the de novo lipogenesis in a mammary gland is different from white adipose tissue and the pathway produces mainly saturated fatty acids with 6-14 carbons [START_REF] Libertini | Purification and properties of a thioesterase from lactating rat mammary gland which modifies the product specificity of fatty acid synthetase[END_REF][START_REF] Hachey | Human lactation. II: Endogenous fatty acid synthesis by the mammary gland[END_REF][START_REF] Mohammad | De novo synthesis of milk triglycerides in humans[END_REF]. This medium-chain fatty acid synthesis in the human mammary gland starts with parturition, irrespective of length of pregnancy [START_REF] Spear | Human mammary gland function at the onset of lactation: medium-chain fatty acid synthesis[END_REF]. Therefore, fatty acids for FAHFAs containing longer acyl chains (e.g. C16:0 and C18:0 for PAHSAs) are acquired from endogenous sources (liver and adipose tissue DNL) or diet and probably only the estolide bond formation takes place in the mammary gland. Preferential handling of fatty acids for TAG or FAHFA synthesis remains enigmatic. Although colostrum contains comparable levels of e.g. TAG-bound essential LA to mature milk, LA-derived polyunsaturated 13-LAHLA appeared only in mature milk. This might be linked to the unique TAG stereo-specific composition of human milk TAGs which contain conserved proportion of saturated FAs (mainly palmitic acid) at the sn-2 position [START_REF] Innis | Dietary triacylglycerol structure and its role in infant nutrition[END_REF].

TAG cores of milk fat globules are synthesized from cytoplasmic lipid droplets budding from the endoplasmic reticulum [START_REF] Lee | Compositional Dynamics of the Milk Fat Globule and Its Role in Infant Development[END_REF][START_REF] Monks | Organellar Contacts of Milk Lipid Droplets[END_REF]. This pathway could produce also TAG estolides that have been detected only in adipose tissue (presumably in lipid droplets) and not in circulation [START_REF] Paluchova | Lipokine 5-PAHSA Is Regulated by Adipose Triglyceride Lipase and Primes Adipocytes for De Novo Lipogenesis in Mice[END_REF][START_REF] Brezinova | Exercise training induces insulin-sensitizing PAHSAs in adipose tissue of elderly women[END_REF][START_REF] Tan | Discovery of FAHFA-Containing Triacylglycerols and Their Metabolic Regulation[END_REF]. ATGL is one of several known FAHFA hydrolases, but this enzyme can also remodel acyl chains on TAGs and synthesize FAHFA-TAG combinations [START_REF] Brejchova | Distinct roles of adipose triglyceride lipase and hormone-sensitive lipase in the catabolism of triacylglycerol estolides[END_REF]. ATGL is highly expressed in goat mammary gland, induced during lactation, and regulates TAG formation [START_REF] Li | Adipose triglyceride lipase regulates lipid metabolism in dairy goat mammary epithelial cells[END_REF]. Therefore, both FAHFAs and TAG estolides can be present within the lipid droplets and their levels can be controlled via ATGL within the epithelial cells or via CEL in the milk. However, it is currently unknown why are FAHFAs and TAG estolides exported to the milk together with their hydrolase CEL.

CEL activity could modulate free FAHFA levels, because TAG estolide concentration is an order of magnitude higher than free FAHFA levels. CEL concentrations are higher in mature milk than in colostrum, but no difference was observed between full term and preterm colostrum [START_REF] Sha | The level of bile salt-stimulated lipase in the milk of Chinese women and its association with maternal BMI[END_REF]. We randomly selected 4 and 4 colostrum samples from the VB and PB group and performed proteomic analysis. Relative protein levels of CEL in colostrum were similar in VB and PB groups (Figure S1) in agreement with published data [START_REF] Sha | The level of bile salt-stimulated lipase in the milk of Chinese women and its association with maternal BMI[END_REF]. However, the CEL activity was reported higher in preterm colostrum, resulting in a higher fat digesting potential during the early postnatal phase [START_REF] Pamblanco | Bile salt-stimulated lipase activity in human colostrum from mothers of infants of different gestational age and birthweight[END_REF]. Therefore, preterm colostrum might be optimized to help with the fat digestion at the expense of TAG estolides and FAHFAs.

The CEL in vitro assay showed that the enzyme hydrolyzes the glycerol ester and the estolide bonds at all positions of TAG estolide substrates. In agreement with data on CEL substrate preferences [START_REF] Kolar | Branched Fatty Acid Esters of Hydroxy Fatty Acids Are Preferred Substrates of the MODY8 Protein Carboxyl Ester Lipase[END_REF], we observed lower affinity to 5-PAHSA than 9-PAHSA isomer. This preference explains high levels of 5-PAHSA in VB colostrum samples and also high 5-PAHSA levels in colostrum from an independent study [START_REF] Brezinova | Levels of palmitic acid ester of hydroxystearic acid (PAHSA) are reduced in the breast milk of obese mothers[END_REF]. In mouse experiments, 5-PAHSA successfully passed the gastrointestinal tract and reached circulation [START_REF] Brezinova | Levels of palmitic acid ester of hydroxystearic acid (PAHSA) are reduced in the breast milk of obese mothers[END_REF][START_REF] Paluchova | Lipokine 5-PAHSA Is Regulated by Adipose Triglyceride Lipase and Primes Adipocytes for De Novo Lipogenesis in Mice[END_REF]. The privileged presence of 5-PAHSA in the colostrum could be linked to its anti-inflammatory effects and protection against ulcerative colitis [START_REF] Lee | Branched Fatty Acid Esters of Hydroxy Fatty Acids (FAHFAs) Protect against Colitis by Regulating Gut Innate and Adaptive Immune Responses[END_REF]. A recent screening study concluded that FAHFAs with lower branching position are more likely to be anti-inflammatory [START_REF] Aryal | Distinct Biological Activities of Isomers from Different Families of Branched Fatty Acid Esters of Hydroxy Fatty Acids (FAHFAs)[END_REF], thus could help to prevent necrotizing enterocolitis in a newborn. Future studies will be needed to define whether FAHFA-rich milk from milk banks could be the preferred substitute for preterm newborns at risk for conditions such as necrotizing colitis [START_REF] Lee | Branched Fatty Acid Esters of Hydroxy Fatty Acids (FAHFAs) Protect against Colitis by Regulating Gut Innate and Adaptive Immune Responses[END_REF]. The degradation products of PAHSAcontaining TAG estolides could also serve as precursors for the odour-active molecules in the milk that attract a newborn [START_REF] Loos | The odour of human milk: Its chemical variability and detection by newborns[END_REF]. Milk TAGs containing 5-hydroxy FA could produce stable δ-lactones, e.g. δ-octadecalactone as semivolatile flavor component in the cream fat [START_REF] Schlutt | Sensory-directed identification of creaminess-enhancing volatiles and semivolatiles in full-fat cream[END_REF][START_REF] Yoshinaga | Analysis of hydroxy triacylglycerol as a lactone precursor in milk fat using liquid chromatography electrospray ionization tandem mass spectrometry[END_REF].

Conclusion

Here we identified TAG estolides as a new component of human milk and broadened substrate specificity of CEL that can modulate milk FAHFA composition. Mothers and health care professional should be aware of the indication that CS and associated delayed breastfeeding initiation could negatively affect levels of anti-inflammatory lipid mediators in the milk.

Limitations of the study.

Qualified personnel processed the colostrum samples at the maternity hospital, but mature milk samples were collected by the nursing mothers at home. Therefore, higher variability in milk composition has to be anticipated (mode and time of sample collection, diet, etc.). Our analytical approach is unable to separate TAG estolide isomers and we report the acyl composition of major species based on MS/MS data. The CEL lipolytic activities were measured under optimal in vitro assay conditions. However, the expected physical form of the TAG estolide substrate is a TAG lipid droplet core of a milk fat globule surrounded by a tri-layer phospholipid membrane. Digestion of the substrate in this unique microenvironment is a complex process and cannot be simulated in vitro.

Clinical perspectives

 Breast milk provides the ideal nutrition for babies, but its composition is very complex and not completely characterized at specific conditions such as preterm birth or caesarean section. Especially bioactive lipids represent an extensive class of poorly explored molecules.  Colostrum levels of 5-PAHSA, a lipid mediator with anti-inflammatory properties preventing colitis, were negatively affected by PB and CS. Furthermore, new FAHFA species and their TAG estolide storage form were identified in the milk. CEL substrate specificity assay towards TAG estolide isomers confirmed a special role of 5-PAHSA in human milk.  Mothers delivering term babies vaginally produce colostrum rich in 5-PAHSA. CS and/or PB delay initiation of de novo lipogenesis, milk production, and FAHFA synthesis. Free 5-PAHSA, resistant to CEL hydrolysis, could contribute to prevention of intestinal inflammation in newborns. Our results could indicate a new approach in the use of donor milk in preterm infants.

Figure 3 :

 3 Figure 3: PAHSAs concentration in human breast milk differ among groups of delivery. (A) The levels of 5-and 9-PAHSA were analysed and compared in milk samples obtained 72 hours (72h) and 28 days (28d) after delivery from mothers who delivered vaginally in term (VB), by caesarean section in term (CS), or preterm (PB). (B) Bioactive FAHFA levels change during lactation period.In PB group, FAHFA levels were measured also at the time when the infants would have reached the 36th week of gestation age (36w). Concentrations of each FAHFA in milk from each time point were compared within each group. Moreover, FAHFA levels in milk samples from PB group at 36w were compared with mature milk from VB group. OAHSA, oleic acid ester of hydroxystearic acids; PAHPA, palmitic acid ester of hydroxypalmitic acid; LAHLA, linoleic acid ester of hydroxylinoleic acid; ⁎, statistically significant at p < 0.05. Data are means ± SEM, n = 27-10 (see Table1).

Figure 4 :

 4 Figure4: TAG estolides are present in milk of all groups and persist during lactation period. The levels of TAG estolides were analysed in milk samples obtained 72 hours (A) and 28 days (B) after delivery from mothers who delivered vaginally in term (VB), by caesarean section (CS) in term, and from mothers who delivered preterm (PB). Data are means ± SEM, n = 26-14 (see Table1).

Figure 5 :

 5 Figure 5: Scheme of TAG estolide substrate and theoretical reactions catalysed by CEL. DAG, diacylglycerol; FA, fatty acid; HSA, hydroxystearic acid; MAG, monoacylglycerol; TAG, triacylglycerol; X:Y stands for number of carbon atoms:number of double bonds.

Figure 6 :

 6 Figure 6: CEL degradation products of TAG EST 70:4 containing either 5-PAHSA (5P) or 9-PAHSA (9P). TAG 5P #1 columns indicate samples, in which TAG EST 70:4 containing 5-PAHSA bound at sn-1/3 alcohol position of glycerol was used as a substrate, TAG 5P #2 represents samples with TAG EST 70:4 and 5-PAHSA at sn-2 position. TAG 9P #1 and #2 mark positional isomers of TAG EST containing 9-PAHSA. (A) Reaction products of LA-glycerol bond cleavage. (B) Reaction products of PAHSAglycerol and estolide bond hydrolysis. (C) Products of estolide bond disruption within the TAG EST molecule. TAG 54:6 stands for trilinolein; DAG EST, diacylglycerol estolide; MAG EST, monoacylglycerol estolide; HSA, hydroxystearic acid; TAG (OH), hydroxylated triacylglycerol; DAG (OH), hydroxylated diacylglycerol; and MAG (OH), hydroxylated monoacylglycerol. Data are means ± SEM, n = 3. T-test was used to compare sn positions (* marks) or regioisomers (# marks), p < 0.05. Regioisomer comparison between panels is described in the text for clarity. No suitable internal standards were available for quantitation of hydroxylated acylglycerols.

Table 1 :

 1 Maternal and infant anthropometrics.

			Term delivery	Preterm delivery
		Vaginal birth (VB) Caesarean section (CS)	Preterm birth (PB)
	Maternal characteristics
	Sample size (n)	27	26	20
	Age at birth (years)	32.4 ± 0.8	32.6 ± 0.8	31.0 ± 0.9
	BMI (kg/m 2 )	22.2 ± 0.5	23.0 ± 0.5	21.7 ± 0.6
	Gestational weight gain (kg)	12.2 ± 1.0	14.4 ± 1.3	8.9 ± 0.9
	Gravidity	2.4 ± 0.2	2.2 ± 0.2	1.5 ± 0.1
	Parity	1.7 ± 0.1	1.7 ± 0.1	1.3 ± 0.1
	Delivery method
	Vaginal	27	0	6
	Caesarean section	0	26	14
	Milk sample
	72 hours	27	26	20
	28 days	18	21	15
	36 th gestational week	-	-	10
	Newborn characteristics
	Gestational age at delivery (wk+d)	40+3 ± 1d	39+2 ± 1d	29+4 ± 1d
	Sex (M/F)	11 / 16	17 / 9	8 / 12
	Birth weight (g)	3613 ± 67	3323 ± 74	1255 ± 65
	Birth height (cm)	50.4 ± 0.3	50.0 ± 0.3	-

  Table 1 values and bar graph data are means ± SEM. Results of One-way ANOVA and Students t-test were considered significant at p < 0.05. Based on published differences in serum 5-PAHSA concentrations in humans (19), the minimum sample size of 19 per group was calculated using G*Power software (t-test, two groups, one tail, power 0.95, α = 0.05)

Acknowledgments:

The authors would like to acknowledge the Metabolomics Core Facility at the Institute of Physiology of the Czech Academy of Sciences for metabolomics and lipidomics profiling and Proteomics service laboratory at the Institute of Physiology and the Institute of Molecular Genetics of the Czech Academy of Sciences for proteomics profiling. The authors wish to thank Jana Kollarova for help with milk sampling and are grateful to all the mothers who took part in the study.

Funding

This work was supported by a grant from the Czech Science Foundation (20-00317S) and by the Czech Academy of Sciences (Lumina quaeruntur LQ200111901).

CRediT

Conceptualization: O.K.; Data curation: K.B., V.P. and T.C.; Formal analysis: K.B., V.P., T.C., M.B. and O.K.; Funding acquisition: O.K.; Investigation: K.B., V.P. and O.K.; Methodology: K.B., V.P. and T.C.; Project administration: M.C., Z.S. and O.K.; Resources: L.B., T.D., M.C. and O.K.; Supervision: M.C. and Z.S.; Validation: T.C. and O.K.; Visualization: K.B. and O.K.; Writing -original draft: K.B.; Writingreview & editing: K.B., V.P., T.C., L.B., T.D., M.C., Z.S. and O.K.;

Data Availability

All raw LC-MS profiling data have been uploaded to the Zenodo repository and can be accessed under accession number: ZENODO #...

Competing Interests

The authors declare that there are no competing interests associated with the manuscript.

Supplementary Information

Supplementary methods

Milk protein mass spectrometry analysis

Protein-containing pellets from the MTBE extraction were processed according to Hartmannova et al. [START_REF] Hartmannova | Acadian variant of Fanconi syndrome is caused by mitochondrial respiratory chain complex I deficiency due to a non-coding mutation in complex I assembly factor NDUFAF6[END_REF]. Briefly, pellets corresponding to 100 µg of protein were solubilized using sodium deoxycholate and digested sequentially with Lys-C and trypsin. LC-MS data were acquired using Exploris 480 (Thermo Scientific) mass spectrometer. Resulting raw files were processed in MaxQuant. Relative amounts of proteins (label-free quantification) were analyzed in Metaboanalyst 5.0 [START_REF] Pang | MetaboAnalyst 5.0: narrowing the gap between raw spectra and functional insights[END_REF]. Samples were log-transformed and auto-scaled.

Supplementary figures

Figure S1: Proteomic profiling of the milk. A) Volcano plot of the identified proteins. False discovery rate set to p-value 0.1 with equal group variance. Volcano plot showed only one significantly altered protein (q-value 0.037) 'Ganglioside GM2 activator' which was higher in PB group. B) Heatmap representation of the 25 top-ranked proteins according to t-test score, distance measure Euclidean, clustering algorithm Ward. Heatmap visualization sorted according to the lowest p-value highlighted a panel of proteins with a higher relative concentration in PB group. Most of the proteins belong to leukocyte mediated immunity, vesicle-mediated transport, and secretion pathways. C) Milk levels of CEL, data are means ± S.D., n = 4. VB, vaginal birth; PB, preterm birth.


Cesar Pairetti

Raphaël Villiers

Stéphane Zaleski
email: stephane.zaleski@sorbonne-universite.fr

On shear layer atomization within closed channels: numerical simulations of a cough-replicating experiment

Aerosol generation during coughing and sneezing has gained major relevance due to the current COVID pandemic. The atomization involved in this process takes place in the complex context of the respiratory system and develops very rapidly. In order to get further insights on the early spray generation, we introduce a simplified model of physiological coughing or sneezing, in the form of a thin liquid layer subject to a rapid (30 m/s) air stream. The setup is simulated using the Volume-Of-Fluid method with octree mesh adaptation, the latter allowing grid sizes small enough to capture the Kolmogorov length scale. The results confirm the trend to an intermediate distribution between a Log-Normal and a Pareto distribution P(d) ∝ d -3.3 for the distribution of droplet sizes in agreement with a previous re-analysis of experimental results by one of the authors. The mechanism of atomization does not differ qualitatively from the multiphase mixing layer experiments and simulations. No mechanism for a bimodal distribution, also sometimes observed, is evidenced in these simulations.

Introduction

Coughing and sneezing are two processes by which a large number of droplets of muco-salivary fluid are exhaled and subsequently travel large distances in the environment [START_REF] Bourouiba | Violent expiratory events: on coughing and sneezing[END_REF][START_REF] Bourouiba | Turbulent gas clouds and respiratory pathogen emissions: potential implications for reducing transmission of covid-19[END_REF]. These phenomena have acquired an acute interest in the context of the so-called aerosol transmission route of the Covid-19 pandemic [START_REF] Balachandar | Host-to-host airborne transmission as a multiphase flow problem for science-based social distance guidelines[END_REF], but have been studied for near a century in the context of respiratory diseases in general [START_REF] Bourouiba | The fluid dynamics of disease transmission[END_REF]. Early investigations by Duguid [START_REF] Duguid | The size and the duration of air-carriage of respiratory droplets and droplet-nuclei[END_REF] and Loudon & Roberts [START_REF] Loudon | Relation between the airborne diameters of respiratory droplets and the diameter of the stains left after recovery[END_REF] of the number and size of droplets emitted in coughing and sneezing events have yielded rich data, incorporating very large numbers of droplets. More recently, Xie et. al. [START_REF] Xie | Exhaled droplets due to talking and coughing[END_REF] proposed corrected droplet size distribution by estimating the in-flight evaporation process. In this context, a comprehensive analysis of the characteristics of these droplets, both in size and velocity, would be of immense interest, as it is a prerequisite for modeling the droplet cloud propagating further downstream from the mouth [START_REF] Chong | Extended lifetime of respiratory droplets in a turbulent vapour puff and its implications on airborne disease transmission[END_REF][START_REF] Dbouk | On coughing and airborne droplet transmission to humans[END_REF].

Rapid photographic imaging [START_REF] Scharfman | Visualization of sneeze ejecta: steps of fluid fragmentation leading to respiratory droplets[END_REF] has revealed features similar to those observed in other atomization processes, including thin liquid sheets, ligaments and droplets. In this context, the statistical distribution of droplet sizes is of particular interest. Although the log-normal distribution has been frequently mentioned in connection with exhalations [START_REF] Nicas | Toward Understanding the Risk of Secondary Airborne Infection: Emission of Respirable Pathogens[END_REF][START_REF] Wells | Airborne Contagion and Air Hygiene. An Ecological Study of Droplet Infections[END_REF] as for other atomizing flows [START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF] many other distributions N(d) of the diameter d have been put forward such as compound gamma distributions [START_REF] Villermaux | Drop fragmentation on impact[END_REF] and many others. A re-analysis of the data of Duguid and Loudon & Roberts for sneezes has however revealed [START_REF] Balachandar | Host-to-host airborne transmission as a multiphase flow problem for science-based social distance guidelines[END_REF] a N(d) ∼ d -2 scaling over an impressive three orders of magnitude. This scaling allows determining the proportion of millimeter-sized droplets that travel short distances and the proportion of much smaller droplets that can be incorporated in a turbulent puff or particle-laden jet and travel long distances as discussed in ref. [START_REF] Bourouiba | Violent expiratory events: on coughing and sneezing[END_REF][START_REF] Balachandar | Host-to-host airborne transmission as a multiphase flow problem for science-based social distance guidelines[END_REF]. The fraction of the exhaled muco-salivary fluid in each class of droplet sizes may also be determined in this way and the probability of having viral loads in each class can be inferred under adequate hypotheses, such as a homogeneous distribution of the virus in volume or surface.

In order to better understand the fluid mechanics of exhalations, King, Brock & Lundell [START_REF] King | Clearance of mucus by simulated cough, physiology[END_REF] have designed a physical model of violent exhalation that may be nicknamed a "cough machine". Air is flowing at high speed (from 10 to 30 m/s) in a rectangular-section duct, with a flow rate analog of the observed human cough. A thin layer of muco-salivary fluid is deposited at the bottom of the duct. While King, Brock & Lundell use the cough machine to observe non-atomizing waves on the thin layer, it can also be used to simulate droplet formation at higher speeds and/or lower velocities. The device would then achieve atomization through a process similar to that of shear flow atomization of planar sheets [START_REF] Gorokhovski | Modeling primary atomization[END_REF][START_REF] Villermaux | Drop fragmentation on impact[END_REF][START_REF] Agbaglah | Numerical simulation of the initial destabilization of an air-blasted liquid layer[END_REF][START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF][START_REF] Zandian | Length-scale cascade and spread rate of atomizing planar liquid jets[END_REF][START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF]. (See also the recent review of numerical approaches in [START_REF] Constante-Amores | Direct numerical simulations of transient turbulent jets: vortexinterface interactions[END_REF].) However there are important differences since the liquid layer is initially at rest and the airflow is impulsive. These differences could lead to different distributions of droplet sizes and velocities, which are the focus of the current investigation.

Numerical method and setup

We model both the muco-salivary fluid and the surrounding gas as incompressible Newtonian fluids. This hypothesis is valid considering that non-newtonian effects on the liquid are only relevant for very small scales. In a context of uniform temperature and high relative humidity, such as a human trachea or a long channel experiment, heat transfer and evaporation will be negligible. Moreover, during the short duration of the atomization process, heat and mass transfer will not have time to develop significantly. Thus the flow is described by the Navier-Stokes equations:

∇ • u = 0 (1) ∂ρ u ∂t + ∇ • (ρ u u) = -∇p + ∇ • (2µD) + σκ n s δ s (2)
where u( x, t) is the velocity field and p( x, t) is the pressure field. Tensor D is equal to 1 2 ∇ u + ∇ u T . ρ and µ are the flow density and viscosity respectively. The last term on the right-hand side on the Navier-Stokes Equation 2 represents the surface tension force, where the coefficient σ is considered constant. The force only acts at the free surface, hence the Dirac function δ s , and also depends on the interface shape, particularly on its curvature κ and normal n s .

We solve this set of equations by a simple finite volume discretization in the onefluid numerical approach for two-phase flow, using a Volume-Of-Fluid (VOF) [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF].

The VOF function at a grid cell is f = cdV, where c indicates the presence of the liquid (c( x, t) = 1) or gas (c( x, t) = 0) phase. The fluid properties at the cell are defined by arithmetic averages:

ρ = f ρ l + (1 -f )ρ g µ = f µ l + (1 -f )µ g (3)
The discrete equations for this model can then be expressed as in [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]:

f n+ 1 2 -f n-1 2 ∆t + ∇ • ( u n f n ) = c c ∇ • u n (4) ρ u * -ρ u n ∆t + ∇ • ρ n+ 1 2 u n u n = ∇ • [µ n+ 1 2 (D n + D * )] + (σκδ s n s ) n+ 1 2 (5) ∇ •       ∆t ρ n+ 1 2 ∇p n+ 1 2       = ∇ • u * (6)
u n+1 = u * - ∆t ρ n+ 1 2 ∇p n+ 1 2 (7)
where c c is the contraction function used in the split volume fraction advection scheme from [START_REF] Weymouth | Conservative volume-of-fluid method for freesurface simulations on cartesian-grids[END_REF] that give consistency to the mass balance. The advective fluxes on Equation ( 5) are integrated explicitly based on a momentum-conserving Bell-Collela-Glaz scheme [START_REF] Pairetti | Mesh resolution effects on primary atomization simulations[END_REF] and the viscous term is computed in a semi-implicit manner. The surface tension force term is based on a discrete-balanced formulation reducing spurious currents [START_REF] Popinet | Numerical models of surface tension[END_REF]. We compute the interface curvature using second-order height functions defined by PLIC plane position based on the analytical formulation from [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF]. Equations ( 6) and ( 7) are the projection steps that will ensure mass conservation for the velocity field at the step (n + 1).

Simulation setup

The physical domain subject to the simulation is a cube of dimensions L 3 , as shown in Figure 1. The flow is channeled through a tube of length x and rectangular cross section with height y and width z . The bottom plate of the tube (z = L/2z /2) is covered with a liquid phase of thickness h, density ρ l and viscosity µ l , depicted in the subfigure (c). The gas phase, modeling the exhaled air, has density ρ a , and viscosity µ a . It enters the flow from the left wall at x = 0 with uniform and constant velocity U through an inflow boundary condition in the tube inner section, |y -L/2| < y /2 and hz /2 < z -L/2 < z /2. This condition ensures the inflow of air blows just above the initial position of the liquid layer. The surface tension of the liquid is noted σ. Gravity is neglected as it is small compared to inertial effects g z /U 2 1. The main dimensionless parameters are the Reynolds number of the air based on the channel height Re a = ρ a U y /µ a = 18, 000, the Weber number of the air based on the channel height We a = ρ a U 2 y /σ = 360, and the Reynolds number of the liquid based on the height of the liquid Re l = ρ l Uh/µ l = 30, 000. The values of these parameters are quite high which makes it surprising that the simulations, performed with a Navier-Stokes code, are still classified as Direct Numerical Simulations (DNS).

Indeed the assessment of the DNS character of the simulation may be inferred from the value of the dissipation in the similar setup of refs. [START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF][START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF]. There the kinetic energy dissipation per unit volume [START_REF] Pope | Turbulent flows[END_REF] was measured and it was observed that the maximum value of /(ρ a U 3 / z ) was about 0.01 which yields an estimate of . The Kolmogorov length scale is then η = (ρ a ν 3 a / ) 1/4 With the value of estimated above, we have η 1 24.2 microns while the size of the smallest grid cell with the thirteen levels of refinement in the simulation reported here is ∆ = 2 -13 L 18.3 microns. According to the DNS resolution criterion given by Pope [START_REF] Pope | Turbulent flows[END_REF], the smallest turbulent scales will be well resolved if ∆/η ≤ 2.1 while we have ∆/η 1 0.75 Even if one uses the more conservative estimate η 2 = y Re -3/4 a as in [START_REF] Constante-Amores | Direct numerical simulations of transient turbulent jets: vortexinterface interactions[END_REF] one gets ∆/η 2 2.8. The perhaps surprising result that the simulations may be qualified as DNS can be explained by the fact that while the less extreme simulations of [START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF][START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF] were limited to z /∆ = 256 grid points in the gas jet thickness, our simulations using octree refinement go up to the equivalent of z /∆ = 1092.

Regarding Adaptive Mesh Refinement (AMR), we use a wavelet-based criterion [START_REF] Van Hooft | Towards adaptive grids for atmospheric boundarylayer simulations[END_REF] to limit local f and u estimated numerical error. The finest simulation reported employs a maximum refinement level of thirteen (L13), resulting on a grid size of 18.3 µm. A cartesian grid with this accuracy would have ≈ 5.5 × 10 11 grid cells while the corresponding octree-grid simulation only scale up to ≈ 3.2×10 9 cells for the spray stage. Simulating 9 ms of physical time take 82 880, 210 446, and 599 040 CPUhours when using 11, 12 and 13 levels respectively. We run all the simulations at Irene machine (TGCC).

Results

Simulations are initialized with zero velocity in the air and liquid although the incompressibility condition results in a non-zero velocity field everywhere in the gas immediately after time zero. The simulation is continued for a total time T=9 ms. The liquid surface is quickly significantly perturbed with waves present over the entire length x of the tube, with a much larger wave near the inlet and some secondary waves near the outlet (Figure 2). The wave stretches into a thin liquid layer that fills with air under the combined effect of pressure and vorticity. Eventually, two mechanisms lead to the formation of droplets: the formation of fingers or ligaments at the end of the sheets and the puncturing of the sheets. The latter causes the formation of holes in the sheets and the subsequent expansion of those holes, leading to the formation of ligaments that eventually pinch and break into droplets. This process has been described in other experimental [START_REF] Scharfman | Visualization of sneeze ejecta: steps of fluid fragmentation leading to respiratory droplets[END_REF] and numerical [START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF] investigations. In the current "closed channel" configuration similar mechanisms of droplet formation are observed. It is seen on Figure 2 that the small dark red droplets near the channel exit have moved quickly since the drag law scaling for a sphere F d ∼ ρ a d 2 U 2 /4 (where d is the droplet radius) implies that the small droplets accelerate much faster than the large ones. A clear phenomenon that couples with droplet formation is growth of a large sheet or wave near the inlet and its progression downstream. The velocity of such waves is typically the Dimotakis velocity [START_REF] Dimotakis | Entrainment and growth of a fully developed, two-dimensional shear layer[END_REF] and it is expressed as

U D = U √ ρ a /( √ ρ a + √ ρ l ).
A similar kind of solitary wave progressing at the Dimotakis velocity has been observed in the simpler setting of an infinite vortex sheet between air and liquid [START_REF] Hoepffner | Self-similar wave produced by local perturbation of the Kelvin-Helmholtz shear-layer instability[END_REF]. The liquid surface is excited by a local perturbation of the flat vortex sheet, and in this case the inhomogeneity of the flow at the entrance plays the role of the localised perturbation, while less localised waves are seen further downstream. In our case the Dimotakis velocity is U D ∼ 1.004 m/s which may be compared to a rough measurement from the simulations (Fig. 3) of U D,num 1.25 m/s.

The identification of droplets or connected fluid components, and the computation of their volumes V d allows to define an equivalent droplet diameter d = (6V d /π) 1/3 . Figure 4 shows the droplet size distribution for the three simulations performed. We compare only the numerical results, for resolved spherical drops, with the experimental data from Xie et. al. [START_REF] Xie | Exhaled droplets due to talking and coughing[END_REF]; in this article, the Probability Density Function of droplet size is assembled with data from several cough samples, considering a correction for evaporation to better represent the size distribution at the mouth. Because of the presence of some non-spherical droplets in the simulation, we filtered the counts, incorporating only near-spherical droplets using a sphericity index. We find that this filtering removes only large droplets, with d > 2mm;there are less than ten of such ligaments. The normalized frequency is shown in Figure 4 together with the statistical error bars, the error being defined as one standard deviation of the binomial. It is seen that the distribution is close to the power law N(d) ∼ d -3.32 for small sizes and then inches down. On Fig. 4 we do not show droplet counts for d < 100 microns and d > 400 microns, since the small sizes are plagued by grid resolution errors and the larger sizes by statistical errors. Regarding the comparison with experiments, the simulations predict a smaller amount of large (d > 200 µm) drops. Nevertheless, both distributions show an acceptable agreement with the Pareto fit proposed.

Discussion

An often considered number frequency (NF) model is the log-normal, which reads N(d) = (B/d) exp -(ln dμ) 2 /(2 σ2 ) where B is a normalization constant, μ is the ex-pected value of ln d, also called the geometric mean, and σ is the standard deviation of ln d, also called the geometric standard deviation (GSD, see [START_REF] Nicas | Toward Understanding the Risk of Secondary Airborne Infection: Emission of Respirable Pathogens[END_REF]). If we plot y = dP(d) versus x = ln d the Log-Normal frequency distribution appears as a parabola. This is done in Fig. 5. That figure shows a fit with a GSD σ = 0.98. The dimensionless GSD is of the same order of magnitude as the GSD ( σ ∼ 1) in similar experiments and simulations [START_REF] Herrmann | On Simulating Primary Atomization Using the Refined Level Set Grid Method[END_REF][START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF][START_REF] Marty | Contribution à l'étude de l'atomisation assistée d'un liquide : instabilité de cisaillement et génération du spray[END_REF], but in the authors' opinion, this may not capture universal underlying physics but the similar range of scales that is within numerical or experimental reach in the literature cited. Comparing the NF at various resolution in [START_REF] Herrmann | On Simulating Primary Atomization Using the Refined Level Set Grid Method[END_REF][START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF] it is seen that as the grid size is reduced the NF shifts to the left, with its geometric mean decreasing and its GSD increasing. This also seen partially in Fig. 4 where the L11 resolution peaks "before" the two others. It is thus possible that as resolution is increased the curved NF seen on Figs 4 and 5 would progressively asymptote to a power law, that is a Pareto, instead of a converged Log-Normal. Such a Pareto NF necessarily has a lower bound, if only the molecular size. This lower bound is not attainable numerically or perhaps even experimentally. A tempting hypothesis is to associate it to the thickness of the sheets as they break or perforate, which involves mechanisms that are not modelled in this study.

Conclusion

We have numerically analyzed a simple physical analog of the physiological mechanism of coughing, based on a previously studied experiment, that is strongly reminiscent of planar sheet atomization processes and presents a reasonable agreement with experimental data [START_REF] Xie | Exhaled droplets due to talking and coughing[END_REF]. A Pareto d -2 distribution was not found but is not excluded at very small diameters. Perspectives include higher resolutions simulations and laboratory experiments in the regime of this numerical experiment, and using the characteristics of the numerically estimated droplet sizes and velocity to predict the further evolution of the droplet cloud using Lagrangian particle methods such as those of Chong et al. [START_REF] Chong | Extended lifetime of respiratory droplets in a turbulent vapour puff and its implications on airborne disease transmission[END_REF].

Figure 1 :

 1 Figure 1: Simulation domain (a), rectangular channel detail (b) and liquid film (c) dimensions.

Figure 2 :

 2 Figure 2: Isometric view of the air-liquid interface at t=9 ms. The interface is colored by the U x velocity component.The droplets by the exit of the tube are close to the stream velocity, as shown in the detail to the right. The left detail employs a different color scale to depict the formation and acceleration of sheets and ligaments in the main wave.

Figure 3 :

 3 Figure 3: In the middle cross-section of the channel, the liquid phase and the solid are shown in black and the gas phase in white. Three different snapshots at regularly spaced time intervals are shown. The oblique line connects the positions, thus giving a graphical display of the wave velocity U D , agreeing with the Dimotakis velocity discussed in the text.

Figure 4 :

 4 Figure 4: Droplet counts in each bin, proportional to the Probability Distribution Function of droplet diameters from Xie et. al. [7] and for the simulations at t = 9 ms and for three grid resolutions L13: ∆x = 18µm, L12: ∆x = 36µm and L11: ∆x = 72µm.

Figure 5 :

 5 Figure 5: Droplet counts in "log-normal coordinates" in which the log-normal NF appears as a parabola, with labels L11-13 as in Fig. 4.

Table 1 :

 1 Dimensional values of the fluid and geometrical parameters (SI units)

Acknowledgements

We acknowledge the funding of the CONICET, the Bec.ar programme from the Ministerio de Educación Argentino, the ERC-ADV grant TRUFLOW, the ANR NAN-ODROP grant funded by the Fondation de France and the PRACE Covid grant of computer time. We thank Lydia Bourouiba, Pallav Kant, Detlef Lohse and Stéphane Popinet for stimulating and fruitful discussions, and for sharing useful material, visualisations and codes.