Annabelle Dupont

Antoine Rauch

PhD Senna Staessens

Mouhamed Moussa

Mickael Rosa

Delphine Corseaux

Emmanuelle Jeanpierre

Julien Goutay

Morgan Caplan

Pauline Varlet

PhD Antoine Pharmd

Rauch

Msc ; Mouhamed Moussa

Guillaume Lefevre

PhD Fanny Lassalle

Anne Bauters

Karine Faure

PhD Marc Lambert

PhD Alain Duhamel

PhD Julien Labreuche

Delphine Garrigue

Simon F De Meyer

Bart Staels

Phd ; Flavien

Natacha Rousse

Eric Kipnis

Peter Lenting

PhD Julien Poissy

Sophie Susen
email: sophiesusen@aol.com

Vascular Endothelial Damage in the Pathogenesis of Organ Injury in Severe COVID-19

Keywords: SARS-CoV-2, endothelium, respiratory failure, renal failure, liver dysfunction, critically ill, death vascular disease, thrombosis, biomarkers, inflammation translational study TOC subcategory: vascular biology neutrophil extracellular trap PAI-1: plasminogen activator inhibitor-1 TFPI: tissue factor pathway inhibitor VV-ECMO: veno-venous extracorporeal membrane oxygenation VWF:Ag: von Willebrand factor antigen

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

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) that causes coronavirus disease 2019 (COVID- 19), preferentially infects respiratory epithelial cells, but can also result in extra-pulmonary manifestations such as myocardial, hepatocellular or acute kidney injury. [START_REF] Zhou | Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study[END_REF][START_REF] Gupta | Extrapulmonary manifestations of COVID-19[END_REF] COVID-19 is also associated with a pro-thrombotic phenotype, elevated laboratory markers of coagulopathy and an increased risk of thrombosis first identified through manifestations of large vessels thrombosis. [START_REF] Connors | COVID-19 and its implications for thrombosis and anticoagulation[END_REF][START_REF] Helms | High risk of thrombosis in patients with severe SARS-CoV-2 infection: a multicenter prospective cohort study[END_REF][START_REF] Klok | Incidence of thrombotic complications in critically ill ICU patients with COVID-19[END_REF][START_REF] Poissy | Pulmonary Embolism in COVID-19 Patients: Awareness of an Increased Prevalence[END_REF][START_REF] Middeldorp | Incidence of venous thromboembolism in hospitalized patients with COVID-19[END_REF][START_REF] Zhang | Deep Vein Thrombosis in Hospitalized Patients With COVID-19 in Wuhan, China: Prevalence, Risk Factors, and Outcome[END_REF] There is also growing evidence suggesting that microvascular thrombosis is likely a major pathophysiologic event in COVID-19. [START_REF] Pons | The vascular endothelium: the cornerstone of organ dysfunction in severe SARS-CoV-2 infection[END_REF][START_REF] Rovas | Microvascular dysfunction in COVID-19: the MYSTIC study[END_REF] The angiotensin-converting enzyme 2 (ACE2) receptor, essential for SARS-CoV-2 uptake by host cells, is highly expressed on endothelial cells. Post-mortem lung histological analyses have reported the presence of vascular inflammation and severe endothelial injury as a direct consequence of SARS-CoV-2 intracellular infection. [START_REF] Varga | Endothelial cell infection and endotheliitis in COVID-19[END_REF][START_REF] Ackermann | Pulmonary Vascular Endothelialitis, Thrombosis, and Angiogenesis in Covid-19[END_REF] Although direct viral tissue damage is a plausible mechanism of endothelial injury, a dysregulated immune response could also be involved. [START_REF] Gupta | Extrapulmonary manifestations of COVID-19[END_REF] Using cultured human pulmonary microvascular endothelial cells, we demonstrated that plasma collected from critically-ill COVID-19 patients triggers endothelial damage in vitro. [START_REF] Rauch | Lille Covid Research Network(LICORNE). Endotheliopathy is induced by plasma from critically-ill patients and associated with organ failure in severe COVID-19[END_REF] Interestingly, this plasma-induced endothelial cytotoxicity was associated in patients with a more pronounced hypoxemia and higher circulating levels of biomarkers of endothelial damage and a dysregulated immune response. This in vitro response could be a marker of in vivo microvascular damage. Once injured, endothelial cells can activate coagulation leading to widespread microvascular thrombi and promote multi-organ failure. To date, reports assessing the presence of endothelial damage and microvascular thrombosis in COVID-19 patients focused on respiratory failure or death as main outcomes. [START_REF] Smadja | Angiopoietin-2 as a marker of endothelial activation is a good predictor factor for intensive care unit admission of COVID-19 patients[END_REF][START_REF] Goshua | Endotheliopathy in COVID-19-associated coagulopathy: evidence from a single-centre, cross-sectional study[END_REF][START_REF] Masi | Systemic Inflammatory Response Syndrome is a Major Contributor to COVID-19-Associated Coagulopathy: Insights from a Prospective SingleCenter Cohort Study[END_REF][START_REF] Liu | The Science Underlying COVID-19: Implications for the Cardiovascular System[END_REF][START_REF] Yeh | Hypercoagulability and COVID-19 associated hypoxemic respiratory failure: Mechanisms and emerging management paradigms[END_REF] To address the question whether endothelial damage and microvascular thrombosis could occur in other organs, we assessed whether the degree of endotheliopathy at intensive care unit (ICU)-admission is an early indicator of organ failure (hepatic, renal, hemodynamic or multi-organ failure) in addition to respiratory failure by investigating a large panel of endothelial damage markers (von Willebrand factor [VWF], plasminogen activator inhibitor-1 [PAI-1], tissue factor pathway inhibitor [TFPI], soluble thrombomodulin [sThrombomodulin]) and endothelial-glycocalyx disruption [syndecan-1]). As a dysregulated immune response in COVID-19 has also been proposed to mediate organ dysfunction through several pathways inducing endothelial damage and microvascular thrombosis, we evaluated if endotheliopathy and organ failure are associated with a dysregulated immune response. [START_REF] Libby | COVID-19 is, in the end, an endothelial disease[END_REF][START_REF] Perico | Immunity, endothelial injury and complement-induced coagulopathy in COVID-19[END_REF] Finally, we investigated the composition of thrombi collected in living critically-ill COVID-19 patients.

Material and methods

Patients

Consecutive patients admitted for COVID-19 were recruited from the intensive-care-unit (ICU) of Lille-University-Hospital between March 21 th and April 16 th , 2020. Inclusion criteria were: individuals aged 18 years or older with a positive SARS-CoV-2 RT-PCR. Patients received enoxaparin or unfractionated heparin according to their renal status and the need for invasive procedures. In overweight and obese patients, the dosing regimen was adapted according to the European Society of Cardiology proposals. [START_REF] Rocca | Antithrombotic therapy and body mass: an expert position paper of the ESC Working Group on Thrombosis[END_REF] From April 6 th onwards, patients received thromboprophylaxis according to the GFHT/GIHP proposals. [START_REF] Susen | Prevention of thrombotic risk in hospitalized patients with COVID-19 and hemostasis monitoring[END_REF] This study was approved by the French institutional authority for personal data protection (Commission Nationale de l'Informatique et des Libertés (CNIL), registration number DEC20-086), and ethics committee (ID-CRB-2020-A00763-36) and informed consent was obtained from all participants. The protocol was registered as a clinical trial (registration number NCT-04327180).

Outcomes

The outcomes were the occurrence of any organ injury, respiratory failure and/or development of renal, hemodynamic, liver injury or multi-organ injury within 14 days of ICUadmission or all-cause 28-day mortality. Secondary respiratory failure was defined by the need for mechanical ventilation. Kidney, liver and shock components of the SOFA score were used to assess the presence of renal, liver or hemodynamic injuries. Renal-failure was defined as presence or any aggravation of renal SOFA and/or KDIGO score from ICU admission; liver-injury as presence or any aggravation of hepatic SOFA score from ICU admission or alanine aminotransferase (ALT) or aspartate transaminase (AST) values more than thrice the upper limit of normal; Hemodynamic dysfunction as a need for epinephrine or norepinephrine support at a dose > 0.1 µg/kg/min. SOFA scores were evaluated at ICUadmission and on ICU days 4, 7 and 14. If the patient was discharged from ICU before day 14, the last SOFA score recorded during ICU stay was used. Multi-organ injury was defined by the presence of more than one dysfunction. Follow-up period was limited to 14 days to minimize the contribution of long-term ICU complications unspecific to COVID-19 as cofounding event in outcome assessment.

Laboratory testing

Blood was collected at ICU-admission for each patient on 0.109M trisodium citrate tube (BD Vacutainer®, BD Diagnostics, Sparks, MD, USA). All analysis were performed on plateletpoor plasma obtained after a double centrifugation at 2500g for 15 min at room temperature. ELISA assays were performed to measure sThrombomodulin, syndecan and nucleosomes using Human Thrombomodulin/BDCA-3 Quantikine® ELISA kit (R&D Systems, Minneapolis, MN, USA), Human sCD138 ELISA kit (Diaclone, Besançon, France) and Cell Death Detection ELISAPLUS kit (Roche Diagnostics, Bâle, Suisse), respectively. TFPI was measured with Asserachrom total TFPI (Diagnostica Stago, Asnières sur Seine, France) and PAI-1 with Asserachrom PAI-1 (Diagnostica Stago). C5a and C5b-9 were measured using standard ELISA kits (Quidel, San Diego, CA, USA). Prothrombin time (PT) was measured on a STA R Max analyser (Diagnostica Stago) using STA Neoplastin R (Diagnostica Stago). Ddimers levels were measured in µg/ml Fibrinogen Equivalent Units (FEU) using an immunoturbidimetric latex-particle assay (Liatest DDI-Plus, Diagnostica Stago) on the STA R Max analyser (Diagnostica Stago). VWF:Ag was measured using an immunoturbidimetric assay (LIAPHEN vWF:Ag, HYPHEN BioMed, Andresy, France), both on a CS 2400 analyser (Sysmex, Kobe, Japan). ADAMTS-13 activity was measured with the Technozym® ADAMTS13 activity ELISA kit (Technoclone, Vienna, Austria). Cell-free DNA was measured by spectrofluorimetry at 520 nm after excitation at 480 nm with Quant-iT PicoGreen dsDNA assay (Thermo fischer Scientific, Waltham, MA, USA) on a SAFAS spectrophotometer (Monaco, France). For MPO-DNA complexes quantification, we used the Cell Death Detection Elisa kit (Roche, Basel, Switzerland). All samples were analysed in duplicate. Interleukin-6 (IL-6), TNFα and soluble interleukin-2 receptor (IL-2R) concentrations were assessed using the Ella Automated Immunoassay System (ProteinSimple, San Jose, CA, USA). Other laboratory blood dosages including a complete blood count, C-reactive protein (CRP), liver transaminases, bilirubin and creatinine by standard methods as part of patient's care in the Biology and Pathology Center of Lille University Hospital.

Thrombi analysis

Histological analysis was performed on thrombi collected immediately after veno-venous extracorporeal membrane oxygenation (VV-ECMO) removal from 8 COVID-19 patients with successful weaning of ECMO and 3 non-COVID-19 patient requiring VV-ECMO(Supplementary 1 for description) . Thrombi were washed in saline, incubated in 4% paraformaldehyde for 24 hours at room temperature and analyzed as previously described. [START_REF] Staessens | Structural analysis of ischemic stroke thrombi: histological indications for therapy resistance[END_REF][START_REF] Laridan | Neutrophil extracellular traps in ischemic stroke thrombi[END_REF] Statistical analysis Quantitative variables were expressed as means (±standard deviation (SD)) or medians (interquartile range, IQR) and categorical variables as frequencies (percentage). We assessed the association of the endothelial markers and nucleosomes with the occurrence of organ injury (censored at 14 days) and death during ICU-stay using Student test and ANCOVA for adjusted analysis on age, sex and body-mass-index (BMI). A log-transformation of original markers data was performed if required. For each outcome, p-values were adjusted for multiple comparisons using false discovery rate (FDR) procedure. We computed the adjusted effect sizes of the 6 markers using the Cohen's d formula for ANCOVA with the 95% confidence intervals (CI). Spearman correlation coefficient was used to test the association between variables and Mann-Whitney or Kruskal-Wallis test for univariate analysis of continuous variables. Statistical analyses were performed using the SAS software package version 9.4 (SAS Institute, Cary, NC, USA).

Results

Patients characteristics

A total of 82 COVID-19 patients were admitted in ICU (Table 1). At admission, median PaO2/FiO2 was 146 mmHg and median SOFA score was 5.9 [2-9]. All patients had increased oxygen requirements with 48 patients (59%) requiring mechanical ventilation, 19 (23%) high-flow oxygen therapy and 15 (18%) low-flow oxygen therapy. Respiratory failure requiring mechanical ventilation was observed in 12 more patients in the 14 days after ICUadmission. Eleven patients (13%) presented with hemodynamic failure while hepatic or renal dysfunction unrelated to pre-existing liver or renal disease were observed in 15 (18%) and 34 (41%) patients respectively. Worsening of the hemodynamic, hepatic or renal status was observed in 12 (15%), 16 (20%) and 6 (7%) patients, respectively. Hepatic dysfunction and hemodynamic failure occurred at a significantly later stage than renal or respiratory failure (p<0.001). Ten patients required a VV-ECMO support and their characteristics are presented in Supplementary Table 1. Eleven patients were receiving heparin (unfractionated heparin, n=8; low molecular weight heparin, n=3).

Endotheliopathy and dysregulated immune response are associated with altered baseline markers of respiratory and liver function

Baseline levels of VWF:Ag (p=0.001), syndecan-1 (p=0.0005), and TFPI (p=0.03) were significantly higher in patients requiring high-flow oxygen therapy or mechanical ventilation at admission than in patients with lower oxygen requirements (Supplementary Figure 1). Bilirubin, AST and ALT correlated with most of the endothelial damage biomarkers (VWF:Ag, syndecan-1, TFPI and PAI-1) (Table 2). We also observed a correlation between either PaO2/FiO2 ratio, ALT or bilirubin and markers of complement activation (C5a) and NETosis. Bilirubin was also correlated with IL-6. AST was correlated with C5a, cytokines and markers of NETosis (cell-free DNA and nucleosomes). Among cytokines, NETs and complement activation markers, nucleosomes were the immune variable that was the most consistently and significantly associated with impaired liver and respiratory functions. No correlation was observed between endothelial markers, cytokines, complement or NETosis markers with creatinine levels. Furthermore, a severe imbalance between VWF and ADAMTS13 concentrations (ADAMTS13/VWF:Ag ratio = 0.17 [0.12-0.26]) was present at admission suggesting an impairment in hemostasis promoting microvascular thrombosis. Overall, when assessing the pattern of organ injury at ICU-admission, the association between immune response and biomarkers of organ dysfunction was mostly apparent for the lungs and the liver.

Occurrence of respiratory failure, liver injury, multi-organ injury and death are associated with endotheliopathy

As shown in Figure 1 and supplementary Table 2, most endothelial damage markers (VWF:Ag, syndecan-1, TFPI and PAI-1) were significantly associated with respiratory failure at admission or occurring during follow-up. After adjustment for age, sex and BMI, these associations remained significant with the strongest effect size for syndecan-1. Overall, this association of endothelial damage at admission with the occurrence of respiratory failures independently of age, sex and BMI, suggests its direct involvement in COVID-19 pneumonia pathogenesis. We further assessed whether the degree of endothelial damage was also associated with the occurrence of specific organ dysfunctions. Occurrence of liver injury was strongly associated with baseline levels of VWF:Ag, syndecan-1, TFPI and PAI-1. After adjustment for age, sex and BMI, theses associations remained significant with the highest effect sizes for TFPI and syndecan-1. As thrombomodulin levels are dependent on renal function in patients without endothelial damage [START_REF] Hergesell | Plasma thrombomodulin levels are dependent on renal function[END_REF] , we evaluated the association of renal failure and sThrombomodulin/creatinine ratio. Only high TFPI and sThrombomodulin/creatinine ratio baseline levels were significantly associated with the occurrence of renal failure. This association remained significant after adjustment. No association was observed between the other endothelial damage markers at admission and the occurrence of hemodynamic failure. Most endothelial damage markers (VWF:Ag, syndecan-1, TFPI and PAI-1) at ICU-admission were significantly associated with the occurrence of multi-organ injury. After adjustment for age, sex and BMI, this association remained significant for all parameters (see Supplementary Table 2 for unadjusted and adjusted for age, sex, and BMI relations). It was also noticed that VWF:Ag, syndecan-1, TFPI and sThrombomodulin were increased according to the number of organ-injuries (supplementary figure 2). VWF:Ag, syndecan-1 and PAI-1 levels at ICU-admission were significantly associated with all-cause-of-death but this association was no longer significant after adjustment.

Interestingly, the 10 patients who required a VV-ECMO support for acute respiratory distress syndrome (ARDS) presented at ICU-admission higher levels of endothelial damage markers (VWF, syndecan-1 and TFPI) and of IL-6 and a tendency for higher nucleosomes levels (p=0.09) (Supplementary Figure 3). Overall, we demonstrated a strong association between endothelial damage markers at ICUadmission and the occurrence of liver injury and multi-organ failure in critically-ill COVID-19 patients.

Endotheliopathy is associated with dysregulated immune response

We further investigated whether cytokine levels, complement activation and NETosis were associated with the intensity of endotheliopathy (Figure 2). The relation between baseline values of endotheliopathy biomarkers and dysregulated immune response at baseline is presented in Supplementary Table 3. All cytokines correlated with VWF:Ag levels. IL-2R and TNFα also correlated with syndecan-1, PAI-1 levels and sThrombomodulin/creatinine levels. All NETs markers correlated with VWF:Ag while cell-free DNA and nucleosomes both correlated with syndecan-1 and TFPI levels. Cell-free DNA was also correlated with PAI-1 levels. When assessing complement activation, only VWF:Ag correlated with C5b-9 levels. We also observed a correlation between most of the immune markers and D-dimer levels (IL-2R: r=0.29, p=0.008; IL-6: r=0.26, p=0.02; cell-free DNA: r=0.27, p=0.02; nucleosomes: r=0.25, p=0.03; and C5b-9: r=0.35, p=0.001), a coagulopathy marker consistently associated with worse outcomes across COVID-19 studies. Overall, dysregulated immune response including excessive NETosis was associated with endothelial damage markers highlighting an interaction between inflammatory-driven processes and endothelial damage leading to organ failure in critically-ill COVID-19 patients.

Immune response and endotheliopathy promote immunothrombosis

In inflammatory disorders, microvascular thrombosis as a consequence of endothelial damage is involved in organ failure. We aimed to investigate the main drivers of thrombogenicity in COVID-19 through the combined analysis of circulating blood biomarkers reflecting endothelial damage and by investigating the composition of in vivo generated thrombi collected in alive patients. We investigated whether the composition of thrombi collected in critically-ill living COVID-19 patients reflected the involvement of both circulating markers of endothelial damage, measured by VWF staining and immune dysregulation, measured by neutrophils and NETs staining. We performed a histological analysis of thrombi collected from VV-ECMO tubings from 8 COVID-19 patients and 3 non-COVID-19 patients presenting with ARDS (for main characteristics of patients included in this histological analysis see Supplementary Table 1). As expected, Martius Scarlet Blue (MSB) and haematoxylin and eosin (H&E) stainings revealed a strong presence of fibrin and VWF. No prominent differences in red blood cells, fibrin, platelet or VWF content between COVID-19 and the control thrombi were observed (Figure 3). In contrast, a striking accumulation of nucleated cellular infiltration was found in thrombi from COVID-19 patients compared to control. Neutrophil elastase staining revealed that leukocyte infiltrates in COVID-19 thrombi consisted predominantly of neutrophils, with few macrophages detected. Quantitative analysis revealed a significant higher area of H3Cit staining in COVID-19 patients compared to non-COVID-19 patients (Supplementary Figure 4). In plasma, neutrophil count correlated with markers of NETosis (cell-free DNA: r=0.51, p<0.0001; nucleosomes: r=0.45, p=0.0005), endothelial damage (VWF: r=0.46, p=0.0001; TFPI: r=0.55, p<0.0001; syndecan-1: r=0.43, p=0.0003; PAI-1: r=0.26, p=0.04) and cytokine (IL-2R: r=0.29, p=0.02). COVID-19 patients displayed higher levels of CRP and fibrinogen then the 3 non COVID-19 patients (295 vs 50 mg/L, p=0.01; 8.62 vs 6.1g/L, p=0.05 respectively). Compared to control thrombi, COVID-19 thrombi contained significant amounts of intracellular and extracellular citrullinated histones H3 (H3Cit), suggestive of NETosis. Overall, we provide evidence for immunothrombosis through endothelial damage and dysregulated immune responses, in particular NETs formation, in the pathogenesis of severe COVID-19. VWF released from damaged endothelial cells interacts with NETs released from neutrophils to provide a scaffold for platelet adhesion and thrombus formation. [START_REF] Kolaczkowska | Molecular mechanisms of NET formation and degradation revealed by intravital imaging in the liver vasculature[END_REF]

Discussion

Using endotheliopathy circulating biomarkers, we demonstrate that endothelial damage is associated with respiratory failure, liver injury, multi-organ failure and death in COVID-19. Most of these associations were still significant after adjustment for age, sex and BMI further underlining the direct involvement of endotheliopathy in organ failure. We also observed a relationship between endothelial damage and the extent of immune inflammatory responses highlighting that inflammatory-driven processes are likely primary drivers of endothelial damage leading to microvascular thrombosis and organ failure in COVID-19. These findings are consistent with our previous finding on the direct in vitro acute cytotoxic effect on vascular endothelial cells of plasma collected from critically-ill COVID-19. [START_REF] Rauch | Lille Covid Research Network(LICORNE). Endotheliopathy is induced by plasma from critically-ill patients and associated with organ failure in severe COVID-19[END_REF] Finally we provide evidence for a major role of neutrophils, NETosis and VWF as a scaffold for thrombus formation and new associative data to support the concept that endothelial damage and dysregulated immune response are involved in thrombus formation and organ failure.

The pathophysiology of multi-organ dysfunction in COVID-19 patients is still poorly understood. COVID-19 patients seem to have a specific respiratory failure pattern characterized by an initial good compliance despite severe hypoxemia. [START_REF] Grasselli | COVID-19 Lombardy ICU Network. Baseline Characteristics and Outcomes of 1591 Patients Infected With SARS-CoV-2 Admitted to ICUs of the Lombardy Region, Italy[END_REF] This discrepancy could be linked to ventilation-perfusion mismatch and the loss of hypoxic pulmonary vasoconstriction, because of endothelial damage with subsequent activation of coagulation and widespread microvascular thrombi formation in the lung and other organs, as well as in the same time hypoperfusion of well ventilated area (REF : Grasselli G et al. Lancet Respir Med 2020) . In the lung, the progressive injury increases vascular permeability then causes the alveolar changes, ground-glass appearance in the parenchymal lung injury, and clearly ARDS. [START_REF] Marini | Management of COVID-19 Respiratory Distress[END_REF][START_REF] Leisman | Facing COVID-19 in the ICU: vascular dysfunction, thrombosis, and dysregulated inflammation[END_REF] Supporting this concept, we observed that endothelial damage markers are significantly associated with the development of respiratory failure in ICU upon admission to ICU. [START_REF] Rauch | Coagulation biomarkers are independent predictors of increased oxygen requirements in COVID-19[END_REF] Endothelial glycocalyx damage was evidenced by increased plasma levels of syndecan-1 and TFPI. As TFPI is mainly associated with endothelial cells through binding to heparan sulfate proteoglycans, a main component of the glycocalyx at the surface of endothelial cells, this result further strengthens the hypothesis that glycocalyx disruption is a key-feature of critically-ill COVID-19 patients. [START_REF] Kato | Regulation of functions of vascular wall cells by tissue factor pathway inhibitor: basic and clinical aspects[END_REF] Endothelial damage was also evidenced by unusually high VWF concentrations indicative of an acute release of endothelial Weibel-Palade bodies content in COVID-19 patients, as previously reported for angiopoietin also stored in the same organelles. [START_REF] Smadja | Angiopoietin-2 as a marker of endothelial activation is a good predictor factor for intensive care unit admission of COVID-19 patients[END_REF] However, it should be mentioned that the presence of endothelial injury is not specific of COVID-19 since the prognostic value of endothelial biomarkers levels was already reported in patients with septic shock or acute lung injury. [START_REF] Hyseni | Active von Willebrand factor predicts 28-day mortality in patients with systemic inflammatory response syndrome[END_REF] Prior reports have also extensively link aberrant NET formation to pulmonary diseases, particularly ARDS. Neutrophils from patients with pneumonia-associated ARDS appear "primed" to form NETs, and both the extent of priming and the level of NETs in blood correlate with disease severity and mortality. [START_REF] Adrover | Programmed 'disarming' of the neutrophil proteome reduces the magnitude of inflammation[END_REF][START_REF] Bendib | Understanding the role of neutrophil extracellular traps in patients with severe pneumonia and ARDS[END_REF] It is therefore likely that NETs, as a source of extracellular histones, contribute to ARDS and sepsis. [START_REF] Chaput | Sepsis: the dark side of histones[END_REF][START_REF] Lefrançais | Neutralizing Extracellular Histones in Acute Respiratory Distress Syndrome. A new role for an endogenous pathway[END_REF] There are increasing reports addressing the pathophysiology of ARDS related to COVID-19 using seasonal influenza patients as controls. To date, most of the data were derived from autopsy studies allowing comparisons with historical controls such as patients died from seasonal influenza. Ackermann et al showed distinctive vascular features between patients died from COVID-19 and influenza with lung autopsy of COVID-19 specimens showing more severe endothelial injury with disrupted endothelial cell membranes and alveolar capillary microthrombi. [START_REF] Ackermann | Pulmonary Vascular Endothelialitis, Thrombosis, and Angiogenesis in Covid-19[END_REF] Likewise, Nicolai et al reported a significantly higher percentage of vessels with feature of immunothrombosis in the lung vasculature in COVID-19 patients than in non-COVID pneumonia ones. [START_REF] Nicolai | Immunothrombotic dysregulation in COVID-19 pneumonia is associated with respiratory failure and coagulopathy[END_REF][START_REF] Nicolai | Vascular neutrophilic inflammation and immunothrombosis distinguish severe COVID-19 from influenza pneumonia[END_REF] While being pivotal in the understanding of the pathophysiology of respiratory failure in COVID-19, the major limitation of such studies is inherent to their post-mortem nature with potential confoundings related to long-term ICU complications unspecific to COVID-19. Our study was designed to overcome such limitation using a dual approach including the dosing of an extensive set of biomarkers related to endothelial damage and NETosis at ICU admission and a thrombus analysis performed on thrombi retrieved from alive COVID-19 patients at ECMO weaning. Our study gives further evidence supporting a critical role for endotheliopathy and immunothrombosis in COVID-19 ARDS pathophysiology since we observed that these biomarkers upon admission were significantly associated with the development not only of respiratory failure but also hepatic and multiorgan failure in ICU. Beyond respiratory failure, kidney and liver injuries were frequently observed in our COVID-19 cohort as previously reported. [START_REF] Yeh | Hypercoagulability and COVID-19 associated hypoxemic respiratory failure: Mechanisms and emerging management paradigms[END_REF][START_REF] Guan | China Medical Treatment Expert Group for Covid-19. Clinical Characteristics of Coronavirus Disease 2019 in China[END_REF][START_REF] Zhang | Liver injury in COVID-19: management and challenges[END_REF][START_REF] Batlle | COVID-19 and ACE2 in Cardiovascular, Lung, and Kidney Working Group. Acute Kidney Injury in COVID-19: Emerging Evidence of a Distinct Pathophysiology[END_REF] The direct relation between biomarkers of endothelial damage, NETosis with different organ dysfunctions has not been established so far during COVID-19 apart from autopsy findings. [START_REF] Varga | Endothelial cell infection and endotheliitis in COVID-19[END_REF][START_REF] Ackermann | Pulmonary Vascular Endothelialitis, Thrombosis, and Angiogenesis in Covid-19[END_REF] It is noteworthy to report that the extent of endotheliopathy and microvascular damage at ICU admission in COVID-19 patients is associated not only with respiratory failure but also with the occurrence of multi-organ failure. Indeed, high baseline levels of most endothelial biomarkers were associated with the development of liver injury while this association was not found for hemodynamic failure and only observed with TFPI levels for renal failure suggesting that different mechanisms are contributing to liver and renal injuries in COVID-19 patients. Moreover, liver injury occurred at a later stage than renal dysfunction. Thrombomodulin is increased in patients with renal dysfunction, a situation frequently encountered in critically-ill COVID-19 patients. [START_REF] Hergesell | Plasma thrombomodulin levels are dependent on renal function[END_REF][START_REF] Batlle | COVID-19 and ACE2 in Cardiovascular, Lung, and Kidney Working Group. Acute Kidney Injury in COVID-19: Emerging Evidence of a Distinct Pathophysiology[END_REF] As such, we only found a week association between the sThrombomodulin over creatinine ratio and the development of renal injury. A renal tropism of SARS-CoV-2 with a direct cytopathic effect on podocytes and renal tubular epithelial cells could explain our results [START_REF] Su | Renal histopathological analysis of 26 postmortem findings of patients with COVID-19 in China[END_REF] , although ultrastructure evidence of direct viral replication in the kidney is lacking.

Liver injury occurs in 14 to 53% of COVID-19 patients according to the definition criteria used across studies and patient's illness severity. [START_REF] Chen | Clinical features and risk factors of COVID-19-associated liver injury and function: A retrospective analysis of 830 cases[END_REF][START_REF] Kullar | Hepatic Injury in Patients With COVID-19[END_REF] Raised transaminase levels are also frequently reported despite the fact that ACE2 is poorly expressed on hepatocytes and that SARS-CoV-2 RNAs were not detected in hepatocytes from COVID-19 patients. [START_REF] Yang | Binding of SARS coronavirus to its receptor damages islets and causes acute diabetes[END_REF][START_REF] Wang | Histopathological features of multiorgan percutaneous tissue core biopsy in patients with COVID-19[END_REF] In contrast, given that ACE2 is highly expressed in the endothelial layer of small hepatic blood vessels, hepatic dysfunction could be related to endotheliopathy rather than to a cytopathic viral effect. The significant association between the degree of endotheliopathy and the development of liver injury is in line with this hypothesis. Liver injury at admission was most frequently observed in patients with elevated cytokines, complement activation and NETosis suggesting that inflammatory-driven processes could promote hepatocyte or liver sinusoidal cell damage in COVID-19 patients. COVID-19 patients with non-alcoholic steatohepatitis (NASH) have a higher likelihood of developing a severe illness, irrespective of obesity. [START_REF] Targher | Risk of severe illness from COVID-19 in patients with metabolic dysfunction-associated fatty liver disease and increased fibrosis scores[END_REF] Interestingly, altered hepatic immune function has been also recently reported to mediate the progression of fatty liver to NASH through pronounced changes in blood immune cell populations. [START_REF] Haas | Transcriptional Network Analysis Implicates Altered Hepatic Immune Function in NASH development and resolution[END_REF] Whether NASH presence exacerbates the SARS-CoV-2-induced inflammatory process thereby contributing mechanistically to severe COVID-19 deserves further research.

We also clearly demonstrated that thrombi are infiltrated with neutrophils and contain markers suggestive of NETs, confirming the role of neutrophils and NETs in the immunothrombotic complications that characterize COVID-19 pathology. [START_REF] Nicolai | Immunothrombotic dysregulation in COVID-19 pneumonia is associated with respiratory failure and coagulopathy[END_REF][START_REF] Nicolai | Vascular neutrophilic inflammation and immunothrombosis distinguish severe COVID-19 from influenza pneumonia[END_REF][START_REF] Skendros | Complement and tissue factor-enriched neutrophil extracellular traps are key drivers in COVID-19 immunothrombosis[END_REF][START_REF] Middleton | Neutrophil extracellular traps contribute to immunothrombosis in COVID-19 acute respiratory distress syndrome[END_REF] We also provide new evidence for a major role of VWF, a binding partner of NETs, as a scaffold for thrombus formation in critically-ill COVID-19 patients. [START_REF] Grässle | von Willebrand factor directly interacts with DNA from neutrophil extracellular traps[END_REF] Albeit being observational, our results strengthen the concept that mechanical vessel obstruction by aggregated NETs is a major pathogenic determinant of severe COVID-19 leading to vascular damage, thrombosis and organ dysfunction. In vitro, it has been shown that sera from COVID-19 patients trigger NET release from control neutrophils. [START_REF] Veras | SARS-CoV-2-triggered neutrophil extracellular traps mediate COVID-19 pathology[END_REF] There is also some evidence that SARS-CoV-2 infection promotes an immunological response unseen with seasonal influenza. This difference is well-illustrated in our ECMO cohort with COVID patients displaying significantly higher levels of CRP and fibrinogen than the 3 patients with non-COVID ARDS). This is in line with the recently published case series also reporting a significantly higher content of NETs in the coronary thrombi of (COVID-19 patients with STEMI compared with historical non-COVID-19 STEMI. [START_REF] Blasco | Assessment of Neutrophil Extracellular Traps in Coronary Thrombus of a Case Series of Patients With COVID-19 and Myocardial Infarction[END_REF] Whether this increased NETosis is due to complement activation or to specific interactions between monocytes and neutrophils or to other mechanisms is currently under investigation. Our results clearly suggest a specific role of neutrophils and NETs in the thrombotic process of circulating blood. A marked imbalance between VWF and ADAMTS13 was present at admission indicative of a reduced capacity of ADAMTS13 to process ultra-large VWF multimers, a finding associated with microvascular thrombosis and organ failure in sepsis. [START_REF] Fukushima | Ratio of von Willebrand factor propeptide to ADAMTS13 is associated with severity of sepsis[END_REF] Given that the pro-thrombotic role of VWF fibers included in NETs is dependent on their multimeric size, potent inhibitors of ADAMTS13 activity such as thrombospondin-1 and αdefensins released locally by activated platelets and neutrophils, respectively could potentiate thrombus growth by VWF fibers. [START_REF] Sorvillo | Plasma Peptidylarginine Deiminase IV Promotes VWF-Platelet String Formation and Accelerates Thrombosis After Vessel Injury[END_REF][START_REF] Bonnefoy | Thrombospondin-1 controls vascular platelet recruitment and thrombus adherence in mice by protecting (sub)endothelial VWF from cleavage by ADAMTS13[END_REF][START_REF] Pillai | Human neutrophil peptides inhibit cleavage of von Willebrand factor by ADAMTS13: a potential link of inflammation to TTP[END_REF] Targeting dysregulated immune responses including neutrophils and NETosis or endotheliopathy might be a valuable and feasible therapeutic approach to prevent thrombosis and improve clinical outcomes of COVID-19. Different targets are currently under investigation. [START_REF] Jin | Endothelial activation and dysfunction in COVID-19: from basic mechanisms to potential therapeutic approaches[END_REF] Drugs that target NETs include molecules that block neutrophils function such as anticytokine therapy (anti-interleukin 6, anti-interleukin 1), that prevent NETosis such as neutrophil elastase, peptidyl arginine deiminase type 4 (PAD4) and that dismantle NETs with DNAses I. [START_REF] Jin | Endothelial activation and dysfunction in COVID-19: from basic mechanisms to potential therapeutic approaches[END_REF][START_REF] Thierry | Neutrophil Extracellular Traps and By-Products Play a Key Role in COVID-19: Pathogenesis, Risk Factors, and Therapy[END_REF] Improving NETs clearance through exogenous DNase treatment with recombinant human DNase I is currently under investigation for safety and efficacy in nine clinical trials in COVID-19. Targeting the VWF-ADAMTS13 axis is also an attractive hypothesis. N-acetylcysteine that dissolves ultra-large VWF multimers under different conditions in both animal models and humans is also currently under investigation. [START_REF] Martinez De Lizarrondo | Abbreviations: CI: confidence interval; ES: effect size; ICU: intensive care unit; PAI-1: plasminogen activator inhibitor-1; s: soluble; TFPI: tissue factor pathway inhibitor, VWF:Ag: von Willebrand factor antigen. Abbreviations: AU: arbitrary units, C5b-9: complement membrane attack complex, IL-2R: soluble interleukin-2 receptor, TFPI: tissue factor pathway inhibitor, TNFα: Tumor Necrosis Factor alpha, VWF:Ag: von Willebrand factor antigen. Correlations were evaluated with Spearman's rank-correlation statistical tool[END_REF] The main strengths of our study are: i) our outcomes were not restricted to respiratory failure or death and other organ dysfunctions were investigated, ii) we assessed the occurrence of organ injury in the first 14 days only after ICU-admission and the composition of extrapulmonary thrombi collected at the acute phase of the disease. Thus our results are less likely to be affected by confounding related to long-term complications unspecific to COVID-19, iii) we combined a large panel of markers and adjusted our results to the main COVID-19 comorbidities to disentangle the role of SARS-CoV-2 on top of comorbidities. Assessing whether the biomarkers associated with multi-organ dysfunction in our study are specific to SARS-CoV-2 or could be translated to other viruses such as SARS, MERS or influenza is a major limitation.

In conclusion our results suggest that organ injury due to endothelial and microvascular damage is associated to immunopathology and may occur in parallel to SARS-CoV-2 intracellular infection.

FiO 2 :Figure1:

 2 Abbreviations: ALT: alanine aminotransferase, AST: aspartate aminotransferase, FiO 2 : fractional inspired oxygen, Il: interleukin, PAI-1: plasminogen activator inhibitor-1, PaO2: arterial oxygen partial pressure, sTM: soluble thrombomodulin, TFPI: tissue factor pathway inhibitor, TNF: tumor necrosis factor, VWF:Ag: von Willebrand factor antigen, p: p-value. * Squared Spearman correlation coefficient

Figure 2 Figure 2 :Figure 3 Figure 3 :

 2233 Figure 2

Table 2 : Correlations between baseline values of all biomarkers and respiratory and liver function criteria

 2

			VWF:Ag	Syndecan-1	TFPI	PAI-1	sTM	C5a	C5b-9	TNFα	IL-6	IL-2R	Neutrophil s	Cell free DNA	MPO-DNA	Nucleosom es
		R*	-0.28	-0.43	-0.33	-0.13	-0.10	-0.25	-0.032	0.006	-0.14	-0.14	-0.36	-0.20	0.039	-0.31
	PaO2
	/FiO2	p	0.01	<0.0001	0.005	0.26	0.38	0.02	0.77	0.96	0.21	0.21	0.004	0.08	0.73	0.008
		R*	0.35	0.34	0.28	0.09	0.046	0.19	0.25	-0.02	-0.009	0.03	0.26	0.17	0.088	0.25
	ALT	p	0.001	0.002	0.02	0.4	0.69	0.09	0.02	0.83	0.93	0.76	0.04	0.15	0.44	0.03
		R*	0.41	0.38	0.37	0.35	0.28	0.24	0,04	0.31	0.21	0.31	0.27	0.38	-0.014	0.27
	AST	p	0.0001	0.0004	0.001	0.001	0.014	0.03	0.7	0.005	0.05	0.004	0.03	0.001	0.9	0.02
	Bilirubine	R*	0.24	0.24	0.27	0.26	0.09	0.37	0.40	0.12	0.30	0.14	0.32	0.44	0.31	0.32
		p	0.03	0.03	0.02	0.017	0.44	0.001	0.0002	0.26	0.006	0.21	0.009	0.0001	0.006	0.005

Abbreviations: ALT: alanine aminotransferase, AST: aspartate aminotransferase,

Acknowledgments

The authors would like to thank all physicians and medical staff involved in patient care. Special thanks are addressed to Eric Boulleaux, Laureline Bourgeois, Vincent Dalibard, Aurélie Jospin, Catherine Marichez, Bénédicte Pradines, Sandrine Vanderziepe, all the biologists and technicians of the Hemostasis Department for their support during the COVID-19 pandemic.

Abbreviations

BMI: body-mass-index ICU: intensive care unit

Sources of funding

This study was supported by the French government through the Programme Investissement d'Avenir (I-SITE ULNE / ANR-16-IDEX-0004 ULNE) and with the Hospital-University Research in Health program (Recherche Hospitalo-Universitaire, WILL-ASSIST HEART ANR-17-RHUS-0011) managed by the Agence Nationale de la Recherche and by KU Leuven through IOF C3/DOA and BOF ISP funding.

Disclosures: none conflict of interest reported

Supplemental Materials

Supplementary Table 1

Values are number (%) otherwise stated: mean (standard deviation) or median [interquartile]. Secondary respiratory failure was defined by the need for mechanical ventilation. Kidney, liver and shock components of the SOFA score were used to assess the presence of renal, liver or hemodynamic dysfunctions. Renal failure was defined as presence or any aggravation of renal SOFA and/or KDIGO score from ICU admission. Liver injury was defined as presence or any aggravation of hepatic SOFA score from ICU admission or ALT or AST values more than thrice the upper limit of normal in accordance with the modified version of the US National Cancer Institute's (NCI's) Common Terminology Criteria for Adverse Events (CTCAE) (2017). Hemodynamic dysfunction was defined as a need for epinephrine or norepinephrine support at a dose > 0.1 µg/kg/min. SOFA scores were evaluated for all ICU-patients at ICU-admission and on ICU days 4, 7 and 14. If the patient was discharged from ICU before day 14, the last SOFA score recorded during ICU stay was used. Multi-organ dysfunction was defined by the presence of more than one dysfunction.


Dissimilarity factor (and co) Algorithm factor Number of clusters factor User factor

Experiments To go further

The variable e↵ect

Medicine 1 : diseases may be classified by etiology (cause), pathogenesis (mechanism by which the disease is caused), or by symptom(s). Alternatively, diseases may be classified according to the organ system involved, though this is often complicated since many diseases a↵ect more than one organ.

And so on. . .

z ẑ intermediate rand(z, ẑ) G 1 = {a, b, c} Ĝ1 = {a, b} A = 2, B = 7 0.6 G 2 = {d, e, f } Ĝ2 = {c, d, e} C = 4, D = 2 Ĝ3 = {f }
Caution: use the adjusted rand index [Hubert and Arabie 1985] to compare rand(z, ẑ) and rand(z, z) when K 6 = K

Stability of a clustering result

Do not forget that ẑ is just an estimate of (a hypothetical true) z

Statistical properties of this estimate should be addressed, as it stability (variance) A simple (but computational demanding) attempt: E↵ect of the metric M (1/5)

X = R 2 , M = ✓ a 0 0 1 ◆ , x 1 = ✓ 0 0 ◆ , x 2 = ✓ 3 0 ◆ , x 3 = ✓ 0 1 ◆ M (x 1 , x 2 ) 2 = (x 1 x 2 ) 0 ✓ a 0 0 1 ◆ (x 1 x 2 ) = a(x 21 x 11 ) 2 = 9a M (x 1 , x 3 ) 2 = (x 1 x 3 ) 0 ✓ a 0 0 1 ◆ (x 1 x 3 ) = (x 32 x 12 ) 2 = 1
E↵ect of the metric M (2/5)

M (x 1 , x 2 ) 2  M (x 1 , x 3 ) 2 , a  1 9
The distance is impacted by the metric, thus the clustering could be also Somewhere the metric is also related to variable selection (try a = 0. . . ) What is "and co"?

E↵ect

Notice also obviously that:

Kernel clustering result depends on the kernel choice Spectral clustering result depends on the Laplacian choice . . .

A meaningful way to be less metric (and co) dependent: idea

Clustering interesting if separated clusters

If separated clusters, partition less metric dependent Thus, the problem is partially reported on choosing K (see later in this lesson)

It will be an interesting element to be used in next lesson also This idea is also applicable for hierarchy, kernel, spectral clustering. . .

A meaningful way to be less metric (and co) dependent: illustration A meaningful way to be less metric (and co) dependent: limit However, it is not always su cient. . . Wk not enough

Spectral clustering K -means clustering

Outline

W K = arg min z with K clusters W (z) W K +1  W K Wn = 0 Figures/XKMEANS-eps-converted-to.pdf

Typology of methods for choosing the number of clusters 5

There exists many other empirical criteria (ex: cross-validation)

There exists clustering methods including automatic choice of K (ex: DBSCAN)

Number of clusters for large data sets

When n increases, K is expected to do so See next lesson for formalizing that point

Gold rule

Retain a useful nb of clusters Some previous criteria are just here for guiding among a set of candidate K values Elbows are interesting for this task

Dendrogram and large data set

Little readable. . . 48 percentages that a given word appears in an e-mail ("make", "you'. . . ) 6 percentages that a given char appears in an e-mail (";", "$". . . )

Transformation of continuous descriptors into binary descriptors

1

  Nosologie méthodique, dans laquelle les maladies sont rangées par classes, suivant le système de Sydenham, & l'ordre des botanistes, par François Boissier de Sauvages de Lacroix. Paris, Hérissant le fils, and co) Algorithm factor Number of clusters factor User factor Experiments To go further Need to compare partitions: empirical error rate Two partitions z and ẑ ⌧ : all permutations on {1, . . . , {d, e, f } Ĝ2 = {a, b, c, d} Dissimilarity factor (and co) Algorithm factor Number of clusters factor User factor Experiments To go further Need to compare partitions: rand index Two partitions z and ẑ A measure on basis of agreement vs. disagreement between object pairs Not limited to the same number of clusters between partitions Rand index [Rand 1971] A: #pairs of elements in x that are in the same subset in z and in the same subset in ẑ B: #pairs of elements in x that are in di↵erent subsets in z and in di↵erent subsets in ẑ C : #pairs of elements in x that are in the same subset in z and in di↵erent subsets in ẑ D: #pairs of elements in x that are in di↵erent subsets in z and in the same subset in ẑ rand(z, ẑ) = A + B A + B + C + D = nb. agree nb. agree + nb. disagree 2 {0, 1} Partitions are closer when rand is high Example

  Use bootstrap samples x (b) (b = 1, . . . , B) Obtain bootstrap partitions z (b) Deduce for instance confidence regions on centers µ through related centers µ (b) Be careful to the permutation of labelling! See the next lesson for more on the statistical properties (need formalizing). . . 12/68 Data factor Dissimilarity factor (and co) Algorithm factor Number of clusters factor User factor Experiments To go further Outline

  of the metric M (3/5) Animals represented by 13 Boolean features related to appearance and activity Large weight on the appearance features compared to the activity features: the animals were clustered into mammals vs. birds Large weight on the activity features: partitioning predators vs. non-predators Both partitions are equally valid, and uncover meaningful structures in the data The user has to carefully choose his representation to obtain a desired clustering E↵ect of the metric 4 M (5/5) Alternative: estimate M (k) by minimizing W M (k) (z) over (z,M (k) ) E↵ect of the linkage criterion (1/3) [A. Jain at al.. Data Clustering: A Review.] E↵ect of the linkage criterion (2/3) [P.-N. Tan at al. (2005). Introduction to data mining, second edition, Addison-Wesley, Chap.8] E↵ect of the linkage criterion (3/3) [P.-N. Tan at al. (2005). Introduction to data mining, second edition, Addison-Wesley, Chap.8]

K

  -means with M = I Not well-separated clusters Well-separated clusters

  Local maxima with K -means: explanation Run K -means from several random centers and keep the best W value Empty clusters with K -means Restart K -means when the empty cluster case occurs What about hierarchical clustering? No problem of starting point No problem of local maxima But the price is strong constraints on nested partitions (see previous lesson) Di↵erent values of K are valid! [P.-N. Tan at al. (2005). Introduction to data mining, second edition, Addison-Wesley, Chap.8] Dissimilarity factor (and co) Algorithm factor Number of clusters factor User factor Experiments To go further

  [A. Jain at al.. Data Clustering: A Review.]

  not well separated Prostate cancer data: chart of individuals sorted by distance to centers visualization In fact it is well separated Outline

  Use the reduced data set without individuals having missing data (n = 475) Use the completed data set where missing data are imputed 3 (n = 506)In both cases, use all mixed variables (not all details at this step, see next lesson)

	Dissimilarity factor (and co) Algorithm factor Number of clusters factor User factor			Experiments To go further
	Prostate cancer data: description 2 Prostate cancer data: partition according to retained variables Prostate cancer data: partition according to recoded variables
	475 patients from 506 (missing values have been discarded) categorical (raw) categorical (MCA)
	err=47.16% 8 quantitative variables, 4 categorical (some are ordinal) variables err=38.95% quantitative categorical (raw) mixing quali/quanti 1 2 1 2
	Two "evident" clusters for medical users: Stage 3 and Stage 4 of cancer 40 Continuous data 5 Categorical data err=9.46% err=47.16% err=8.63% 1 2 1 2 1 Stage 3 142 131 175 98 2 Stage 3 247 26 142 131 252 Stage 4 120 82 87 115 Data set completed data reduced data 21 Stage 4 19 183 120 82 20 182 err 12.8 8.1
	30
									4
	20
									3
	10
	-20 -10 0 Partition varies with retained variables as expected 2nd axis PCA 0 1 2 2nd axis MCA A general principle: categorical variables less informative than quantitative ones
	-30 However, categorical variables here improve quantitative ones -1
	-40
	-80 -50	-60	-40	-20	0	20	40	60	-2.5 -2	-2	-1.5	-1	-0.5	0	0.5	1
				1st axis PCA							1st axis MCA
	2 Byar and Green (1980)

MCA is equivalent to recoding categorical variables Raw data and MCA data are in a one-to-one mapping (no info. loss) It can however drastically impact clustering result It open the question of data units/coding to use Currently: let the user to choose the unit (prior or posterior choice) Next lesson: need formalizing to go further Prostate cancer data: partition according to missing data It is current to have a data "pretreatment" like missing data imputation Be careful: it can impact the clustering Imputation gives only an estimate data set x which is a "deteriorated" data set As a consequence it can lead to a "deteriorated" clustering result See next lesson to formalize this problem

Prostate cancer data: variable detail

We use the mice package:http://cran.r-project.org/web/packages/mice/mice.pdf

E↵ect of the metric M (4/5)

Figures from A.K. Jain (2008). Data Clustering:

Years Beyond K-Means.

Mirkin, Boris. (2011). Choosing the number of clusters. Wiley Interdisc. Rew.: Data Mining and Knowledge Discovery. 1. 252-260.

There are 3 other continuous descriptors we do not use

https://archive.ics.uci.edu/ml/machine-learning-databases/spambase/

Byar DP, Green SB (1980): Bulletin Cancer, Paris 67:477-488

Prostate cancer data: cluster weight

Some criteria to estimate K (1/3)

The first "elbow" on the W K curve [Hartigan, 1975] K = min

The Gap statistics measures the gap with uniformity [Tibshirani et al., 2001] Gap

K the within cluster sum of squares from a bth uniform data set with same range as the original data

Some criteria to estimate K (2/3)

Form of an ANOVA F-statistic [Calinski and Harabasz, 1974] K = arg max

Measure of how well the all x i are clustered [Kaufman and Rousseeuv, 1990] :

See experiments on the US census in Part 1

In particular, we can observe the e↵ect on the partition from:

the linkage criterion for HCA the algorithm (HCA vs. K -means)

Next lesson

To go further towards clustering evaluation, there is a need to further formalize. . .

Clustering : une vision unifiée pour une utilisation éclairée

Formalisation par modèles de mélange


Formalizing estimation Formalizing selection More advanced formalizing Experiments To go further The model-based clustering paradigm a cluster () a distribution It recasts all previous/next questions into model design/estimation/selection It takes benefits from all theoretical statistics environment How to choose the best metric M (k) ?

How to choose the number K of clusters?

Clusters of di↵erent sizes are they well estimated?

How to choose the data unit?

How to select features?

How to deal with mixed data?

How to deal with missing data?

How to deal with outliers?

. . .

Formalizing estimation

Formalizing selection More advanced formalizing Experiments To go further

What about empirical clustering?

Somewhere it works pretty well even if it has the previous mentioned limits Reformulate K -means: the hidden Gaussian assumption

W I (z) = n X i=1 K X k=1 z ik kx i µ k k 2 I = 2 n X i=1 K X k=1 z ik ln 2 6 6 6 4 1 K |{z} ! 1 (2⇡) d/2 |I| 1/2 exp ✓ 1 2 (x i µ k ) 0 I(x i µ k ) ◆ | {z } N d (µ k ,I)

Estimation of ✓ by complete-likelihood

Maximize the complete-likelihood over (✓, z)

`c (✓; x, z) = n X i=1 K X k=1 z ik ln {⇡ k p(x i ; ↵ k )} Equivalent to traditional methods Metric M = I M free M k free Gaussian model [⇡ I ] [⇡ C ] [⇡ k C k ]

Principle of EM and CEM

Initialization: ✓ 0 Iteration n o q:

Step E: estimate probabilities t q = {t ik (✓ q )}

Step C: classify by setting t q = MAP({t ik (✓ q )})

Step M: maximize

✓ q+1 = arg max ✓ `c (✓; x, t q )
Stopping rule: iteration number or criterion stability Gaussian M-step Note : low at the beginning but increase of the log-likelihood -80

n (q) k = n X i=n l +1 t ik (✓ (q) ) ⇡ (q+1) k = n (q) k n µ (q+1) k = 1 n (q) k n X i=1 t ik (✓ (q) )x i ! ⌃ (q+1) k = 1 n (q) k n X i=1 t ik (✓ (q) )(x i µ (q+1) k )(x i µ (
µ 1 µ 2 Likelihood -160 -140 -120 -100 -2 0 2 -2 -1 0 1 2 1EM essai 1 -2 0 2 -2 -1 0 1 2 essai 2 -2 0 2 -2 -1 0 1 2 essai 3 -2 0 2 -2 -1 0 1 2 essai 4 -2 0 2 -2 -1 0 1 2 essai 5 -2 0 2 -2 -1 0 1 2 essai

Categorical variables: latent class model

Categorical variables: d variables with m j modalities each, x j i 2 {0, 1} m j and Integer: Poisson mixture model

x jh i = 1 , variable j of x i takes level h Conditional independence: p(x i ; ↵ k ) = d Y j=1 m j Y h=1 (↵ jh k ) x jh i and ↵ jh k = p(x jh i = 1|z ik = 1) with ↵ k = (↵ jh k ; j = 1, . . . , d; h = 1, . . . ,
integer variables: d variables x j i 2 N Intra conditional independence: p(x int i ; ↵ int k ) = d Y j=1 (↵ j k ) x j i ↵ j k ! e ↵ j k
SPAM E-mail Database 3 n = 4601 e-mails composed by 1813 "spams" and 2788 "good e-mails"

d = 48 + 6 = 54 continuous descriptors 2
48 percentages that a given word appears in an e-mail ("make", "you'. . . ) 6 percentages that a given char appears in an e-mail (";", "$". . . )

Transformation of continuous descriptors into binary descriptors

x j i = ⇢ 1 if

Mixed data: conditional independence everywhere

The aim is to combine continuous, categorical and integer data

x 1 = (x cont 1 , x cat 1 , x int 1 )
The proposed solution is to mixed all types by inter-type conditional independence

p(x 1 ; ↵ k ) = p(x cont 1 ; ↵ cont k ) ⇥ p(x cat 1 ; ↵ cat k ) ⇥ p(x int 1 ; ↵ int k )
In addition, for symmetry between types, intra-type conditional independence

Discarded solutions

Suppress units and/or variables with missing data ) loss of information Imputation of the missing data by the mean or more evolved methods ) uncertainty of the prediction not taken into account

Retained solution

Use an integrated approach which allows to take into account all the available information to perform clustering

Assumption on the missingness mecanism

Missing At Randon (MAR): the probability that a variable is missing does not depend on its own value given the observed variables.

Missing data: maximum likelihood estimator

Observed log-likelihood. . .

`(✓; x O ) = n X i=1 log K X k=1 ⇡ k p(x O i ; ↵ k ) ! = ln 2 6 6 6 4 K X k=1 ⇡ k Z x M i p(x O i , x M i ; ↵ k )dx M i | {z
z + ik = P(Z ik = 1|x O i ; ✓) = ⇡ k (x O i ; ⌃ k ) P K `=1 ⇡ ` (x O i ; ⌃ `) x M + ik = E h X M i x O i , Z ik = 1; ✓ i = µ M ik + ⌃ MO ik ⇣ ⌃ OO ik ⌘ 1 (x O i µ O ik ).
⇡ + k = 1 n + k n X i=1 z + ik , µ + k = 1 n + k n X i=1 z + ik x + ik ⌃ + k = 1 n + k n X i=1 z + ik h (x + ik µ + k )(x + ik µ + k ) 0 + ⌃ + ik i where n + k = P n i=1 z + ik , x + ik = ✓ x O i x M + ik ◆ , ⌃ + ik = ✓ 0 O i 0 OM i 0 MO i ⌃ M + ik ◆ with 0 the d ⇥ d
null matrix, and

⌃ M + ik = ⌃ MO ik ⌃ O ik 1 ⌃ OM ik .
Interpretation of ⌃ M + ik Variance correction due to the under-estimation of variability caused by the imputation of missing data.

Missing data: SE algorithm

A SE algorithm estimates then (x M , z M ) Iteration nb q:

E-step: compute conditional probabilities p(x M , z M |D; ✓) S-step: draw (x M(q) , z M(q) ) from p(x M , z M |D; ✓)
Stopping rule: iteration number Properties simplicity because of conditional independence the mean/mode of the sequence (x M(q) , z M(q) ) estimates (x M , z M ) confidence intervals are also derived

⇡ 1 = ⇡ 2 = 0.5, µ 1 = (0, 0) 0 , µ 2 = (2, 2) 0 , ⌃ 1 = ⌃ 2 = I M = {spherical,

Some heuristics entropy-based criteria: examples

A fundamental decomposition of `(✓; x): for any "fuzzy partition" c = {c ik }

`(✓; x) = n X i=1 K X k=1 c ik ln{⇡ k p(x i ; ↵ k )} n X i=1 K X k=1 c ik ln t ik (✓) = `c (✓; x, c) + ⇠(✓; c) = complete-data log-likelihood+entropy NEC criterion (Normalized Entropy Criterion): retain m minimizing NEC K = 8 < : ⇠( ✓K ; t( ✓K )) `( ✓K ; x) `( ✓1 ; x) if K > 1 1 i f K = 1 CL criterion (Completed Likelihood): retain m maximizing CL = `c ( ✓; x, ẑ) = `( ✓; x) | {z } model adequacy ⇠( ✓; ẑ) | {z } partition evidence
Behaviour: not completely satisfactory but something happens. . .

Theoretical model selection criteria

The most widespread principle

BIC criterion: consistency

Consistency: BIC asymptotically selects the best

m ⇤ = arg inf m2M KL(p, p ✓ ⇤ m ) Theoretical illustration of consistency: m 1 ✓ m 2 , m 1 being the true model, ⌫ = ⌫ 2 ⌫ 1 , `= `( ✓2 ; D) `( ✓1 ; D), we have 2(BIC 2 BIC 1 ) + ⌫ ln(n) = 2 `d ! 2 ⌫ With µ = ⌫ and 2 = 2 ⌫ the mean and the variance of 2 ⌫ p( 2 ⌫ > ⌫ ln(n))  p(| 2 ⌫ µ| > ⌫ ln(n) µ)  2 ( ⌫ ln(n) µ) 2 n!1
! 0 by using the Chebyschev inequality. Thus, asymptotically, BIC will select m 1 Special case of K : be careful on the 2 approximation validity. . . Theoretical result:

Preliminaries: n = n(✓ ⇤ 2 ✓ ⇤p 2 ) 0 J(✓ ⇤ 2 )(✓ ⇤ 2 ✓ ⇤p 2 ), J(✓ ⇤
2 ) the Fisher matrix for a data unit calculated with the true parameter ✓2 and ✓ ⇤p 2 its projected value on the parameter subspace associated to the one component case,

µn = E[ 2 ⌫ ( n )] = ⌫ + n , 2 n = Var[ 2 ⌫ ( n )] = 2( ⌫ + n ) Asymptotically: by Chebishev inequality, with µn ⌫ ln n 2n ln 2 > 0 p(choose wrong model) = p(ICLbic2 < ICLbic1)  2 n (µn ⌫ ln n 2n ln 2) 2
Thus it goes towards 0 for well-separated groups

Experimental result: 100 samples from a univariate Gaussian mixture Large n: BIC behaviour (1/2)

⇡ 1 = ⇡ 2 , µ 1 = 0, µ 2 = µ, 2 1 = 2 2 =
The mixture density is wrong (as all models)

Mixtures allow to estimate any distribution by increasing the number of components (high flexibility) Units: changing the data units

Principle of data units transformation u:

u : X = X id ! X u x = x id = id(x) 7 ! x u = u(x) u is a bijective mapping to preserve the whole data set information quantity We denote by u 1 the reciprocal of u, so u 1 u = id Thus, id is only a particular unit u Often a meaningful restriction7 on u: it proceeds lines by lines and rows by rows u(x) = (u(x 1 ), . . . , u(xn)) with u(x i ) = (u 1 (x i1 ), . . . , u d (x id ))

Advantage to respect the variable definition, transforming only its unit u(x i ) means that u applied to the data set x i , restricted to the single individual i u j corresponds to the specific (bijective) transformation unit associated to variable j In addition, for symmetry between types, intra conditional independence for each Results:

New units u SZ and u AP are selected by ICL New units allow to select two groups and provides a lower error rate

1

  Interesting to understand why 2 Interesting to overcome their limits then In fact many empirical methods are hidden model-based clustering ones!

⇡

  Gaussian with variance matrix I and same cluster sample size (see later)Reformulate K -means: the hidden estimate choiceW I (z) = 2`c (µ; x, z) + cst + ẑW = arg minz W I (z) (ẑ `c , μ`c ) = arg max (z,µ) `c (µ; x, z) ) ẑW ⌘ ẑ`cEstimateMaximum of the so-called complete-likelihood (see later for its statistical properties) k p(x 1 ; ↵ k ) Mixture parameter:✓ = (⇡, ↵) with ↵ = (↵ 1 , . . . , ↵ K )Model: it includes both the family p(•; ↵ k ) and the number of groups Km = {p(x 1 ; ✓) : ✓ 2 ⇥}The number of free continuous parameters is given byThe clustering process in mixturesEstimation of ✓ by ✓ Estimation of the conditional probability thatx i 2 G k t ik ( ✓) = p(Z ik = 1|X i = x i ; ✓) = ⇡k p(x i ; ↵k ) p(x i ; ✓)Estimation of z i by maximum a posteriori (MAP) ẑik = I {k=arg max h=1,...,K t ih ( ✓)}

  borderline with [⇡ k C k ] Importance of model selection: bias/variance trade-o↵ Partition error rate: err(z 1 , z 2 ) 0 a distance-like between two partitions z 1 , z 2 Gap between true and model partition: Importance of model selection: illustration of the variance e↵ect 30 samples from a bivariate mixture with two components

  Some heuristics entropy-based criteria: CL illustrationInterpretation as a penalized within-cluster inertia criterion:CL ([p I ],K ) = nd 2 ln(W K ) n ln(K )

  Thus models leading to overlapping groups are more penalized (low z evidence) ICL criterion (Integrated Classification Likelihood): replace z by ẑ ICL = ln p(x, ẑ|m) BIC-like approximation of ICL: ln p(x, z|m) = ln p(x, z; ✓x,z|m) ⌫ 2 ln n + Op(1) In case of the right model m: ✓x,z a.s. ! ✓ ⇤ and ✓x a.s. ! ✓ ⇤ . Thus, for n large enough, ✓x,z ⇡ ✓x. Then, we take ẑ = MAP( ✓x) (or also ẑ = t( ✓x)). It gives ICLbic = lnp(x, ẑ; ✓x) ICL criterion: robustness to model misspecification A bivariate mixture of a uniform and a Gaussian cluster: non-Gaussian component: ⇡1 = 0.5, p 1 (x1) = 0.25 I [ 1,1] (x 1 ) I [ 1,1] (x 2 ) Gaussian component: ⇡2 = 0.5, µ2 = (3.3, 0) 0 , ⌃2 = I 50 simulated data sets of size n = 200 model with two groups and parameter ✓ ⇤ 2

  Model m: full mixed data x = (x cont , x cat , x ordi , x int , x rank ) (missing data are allowed also) are simply modeled by inter conditional independence p(x;↵ k ) = p(x cont ; ↵ cont k ) ⇥ p(x cat ; ↵ cat k ) ⇥ p(x ordi ; ↵ ordi k ) ⇥ . . .

  m j )

	Need to formalize	Formalizing estimation	Formalizing selection	More advanced formalizing	Experiments	To go further
				20/102

  i ✓ {1, . . . , d} the set of the observed variables from sample i x O i the observed data from sample i M

	Formalizing estimation	Formalizing selection	More advanced formalizing	Experiments	To go further
	Missing data: notations and MNAR assumption

O i the set of the missing variables for sample i µ O ik the sub-vector of µ k associated to index O i (the same for M i ) ⌃ OM ik the sub-matrix of ⌃ k associated to row O i and columns M i (the same for any other combination)

  }

	Need to formalize	Formalizing estimation	Formalizing selection	More advanced formalizing	Experiments	To go further
		EM and Gaussian case: E step
	✓ and ✓ + the parameters for two successive steps (idem for missing data)
				MAR assumption
					40/102

Table :

 : Raw units: 11% misclassified

	Need to formalize	Formalizing estimation	Formalizing selection	More advanced formalizing	Experiments	To go further
		Outline Outline Some remaining questions
	Need to formalize Need to formalize
	Formalizing estimation Formalizing estimation
	More on dependent data (like times series)
	Formalizing selection Formalizing selection High-dimensional data
	Missing not at random data (MNAR)
	More advanced formalizing More advanced formalizing . . .
	Experiments Experiments
				clusters
	To go further To go further		1 270	2 22
				23	191
				Table: New units: 9% misclassified
		81/102 100/102			79/102

101/102

Impossible to use BIC/ICL for kernel/spectral clustering (data set has changed. . . )

http://cran.r-project.org/web/packages/capushe/ 67/102

Reformulate K -means: elbow as a slope heuristics (2/3)

Possibility to relax this restriction, including for instance linear transformations involved in PCA (principal component analysis). But the variable definition is no longer respected.

This dataset is extracted from the fda package(Ramsay, Wickham, Graves, and Hooker 2018). It contains daily temperature and precipitation at 35 di↵erent locations in Canada averaged over 1960 to 1994.Data are contained in a list of 3 elements: tempav, a 365×35 matrix where one column corresponds to the daily temperature, precav a 365 × 35 matrix where one column corresponds to the daily precipitation and time, a vector containing the index of the day.13 Ramsay JO, Silverman BW(2005). Functional Data Analysis. Springer Series in Statis-tics, 2nd edition. Springer. ISBN 038740080X.

A SEM algorithm to estimate ✓ by maximizing the observed-data log-likelihood Initialisation: ✓ (0) Iteration nb q: E-step: compute conditional probabilities p(x M , z|D; ✓ (q) ) S-step: draw (x M(q) , z (q) ) from p(x M , z|x 0 ; ✓ (q) ) M-step: maximize ✓ (q+1) = arg max ✓ ln p(x O , x M(q) , z (q) ; ✓)

Stopping rule: iteration number Properties: simpler than EM and interesting properties! Avoid possibly di cult E-step in an EM Classical M steps

Avoids local maxima

The mean of the sequence (✓ (q) ) approximates ✓

The variance of the sequence (✓ (q) ) gives confidence intervals

Keep in mind

George E.P. Box (1987) "Essentially, all models are wrong, but some are useful" M = {m} will denote the set of competing models

The true distribution p is not necessarily in M Density estimation: AIC, BIC Clustering: ICL, CL, NEC Units: revisiting units as a modelling component Explicitly exhibiting the "canonical" unit id in the model

Thus the variable space and the probability measure are embedded

As the standard probability theory: a couple (variable space,probability measure)!

Changing id into u, while preserving m, is expected to produce a new modelling

A model should be systematically defined by a couple (u,m), denoted by p u m Units: interpretation and identifiability of p u m Standard probability theory (again): there exists a measure u 1 (m) s.t. 8

There exists two alternative interpretations of strictly the same model: Two points of view:

Statistician

The model p u m is not identifiable over the couple (m, u)

Practitioner

Freedom to choose the interpretation which is the most meaningful for him 8 This set is usually restricted to a single element Example for abs. continuous x and di↵erentiable u, the density transform in id is:


Modeling

Estimating Selecting Experiments To go further HD data: definition (2/2)

An attempt in the parametric case Dataset x = (x 1 , . . . , xn), x j described by d variables and a model m with ⌫ parameters, where n = o(g (⌫)), with g a given function Justification:

We consider the heteroscedastic Gaussian mixture with of true parameter ✓ ⇤ with K ⇤ components. We note ✓ the Gaussian MLE with K ⇤ components. We have g linear from the following result [Michel,08]: it exists constants , A and C such that

Ex[Hellinger 2 (p ✓ ⇤ , p ✓ K )]  C   ⌫ n n 2A ln d + 1 ln ⇣ 1 ^h ⌫ n A ln d i⌘o + 1 n .
But ⌫ can be high since ⌫ ⇠ d 2 /2, combined with potentially large constants.

HD density estimation: curse

A two-component d-variate Gaussian mixture:

⇡ 1 = ⇡ 2 = 1 2 , X 1 |z 11 = 1 ⇠ N d (0, I), X 1 |z 12 = 1 ⇠ N d (1, I)
Components are more and more separated when d grows: kµ 2 µ 1 k I = p d. . . A two-component d-variate Gaussian mixture:

⇡ 1 = ⇡ 2 = 1 2 , X 1 |z 11 = 1 ⇠ N d (0, I), X 1 |z 12 = 1 ⇠ N d (1, I)
Each variable provides equal and own separation information Theoretical error decreases when d grows: err theo = ( p d/2). . . Many variables provide no separation information Same parameter setting except:

FDA -4 -3 -2 -1 0 1 2 3 4 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 1st axis FDA 2nd axis FDA d=2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -3 -2 -1 0 1 2 3 4 1st axis FDA 2nd axis FDA d=20 -1.5 -1 -0.5 0 0.5 1 1.5 -4 -3 -2 -1 0 1 2 3 4 5 1st axis FDA 2nd axis FDA d=200 -1.5 -1 -0.5 0 0.5 1 1.5 -3 -2 -1 0 1 2 3 1st axis FDA
X 1 |z 12 = 1 ⇠ N d ((1 0 . . . 0) 0 , I)
Groups are not separated more when d grows: Many variables provide redundant separation information Same parameter setting except:

kµ 2 µ 1 k I = 1. . . -4 -3 -2 -1 0 1 2 3 4 -5 -4 -3 -2 -1 0 1 2 3 4 x1 x2
X j 1 = X 1 1 + N 1 (0, 1) (j = 2, . . . , d)
Groups are not separated more when d grows:

kµ 2 µ 1 k ⌃ = 1. . . -3 -2 -1 0 1 2 3 4 -6 -4 -2 0 2 4 6 x1
⇡ 1 = ⇡ 2 = 1 2 , X 1 |z 11 = 1 ⇠ N d (0, ⌃), X 1 |z 12 = 1 ⇠ N d (1, ⌃)
Each variable provides equal and own separation information But which kind of parsimony?

Remember that clustering is a way for dealing with large n

Why not reusing this idea for large d? Gaussian mixture of factor analysers (2/2) Use alternatively a lasso-like procedure for ranking quickly di↵erent sets of clustering related and clustering independent variables [Sedki et al.,14] crit

Co-clustering

⌃ k = B k B 0 k + ! k ⇤ k where B k is a loadings d ⇥ q non-square real matrix (1  q  qmax, qmax < d) ! k is a positive real number ⇤ k is a d ⇥ d diagonal positive definite matrix such that |⇤ k | = 1 Interpretation X 1 2 R d is generated by a latent variable Y 1 2 R q X 1|Y 1 ,Z 1k =1 = B k Y 1 |{z} factor +µ k + " k where Y 1 ? " k and Y 1 ⇠ Nq(0, I) and " k ⇠ N d (0, ! k ⇤ k ) 15/
Complexity ⌫ = (K 1) + Kd + Kq(d (q 1)/2) + Kd Model selection Models in competition: m = (q, K ) +
⌃ k = D k k D 0 k where D k is the orthogonal matrix of the eigenvectors of ⌃ k k is a diagonal matrix containing the eigenvalues of ⌃ k k = 0 B B B B B B B B B @ a k1 0 . . . 0 a k k 0 0 b k 0 . . . 0 b k 1 C C C C C C C C C A 9 = ; k 9 = ; (d k ) with a kj b k , for j = 1, ..., k and k < d Complexity ⌫ = (K 1) + Kd + K X k=1 k [d ( k + 1)/2] + K X k=1 k + 2K
; ✓) = ( K X k=1 ⇡ k p(x S 1 ; µ k , ⌃ k ) ) | {z } clustering variables ⇥ n p(x U 1 ; a + x R 1 b, C) o | {z } redundant variables ⇥ n p(x W 1 ; u, V) o | {z } independent
,⇢ = `(✓; x) K X k=1 d X j=1 |µ kj | ⇢ K X k=1 d X (j,j 0 ),j6 =j 0 |(⌃ 1 k ) jj 0 |
where ✓ full Gaussian parameters, x is x centered and Example in the continuous case:

x 2 R n⇥d
It could be possible to define a within-block inertia criterion

W (z, w ) = K X k=1 L X l=1 n X i=1 d X j=1 | {z } P i,j,k,l z ik w jl kx ij µ kl k 2
with µ kl the center of the block (k, l)

µ kl = 1 n kl X i,j z ik w jl x ij
where n kl = P ij z ik w jl is the sample size of the block (k, l)

But we know now that it hides some model-based assumptions. . . Distribution for di↵erent kinds of data [Govaert and Nadif, 2014] The pdf p(•; ↵z i w j ) depends on the kind of data x ij : Binary data:

x ij 2 {0, 1}, p(•; ↵ kl ) = B(↵ kl )
Categorical data with m levels:

x ij = {x ijh } 2 {0, 1} m with P m h=1 x ijh = 1 and p(•; ↵ kl ) = M(↵ kl ) with ↵ kl = {↵ kjh } Count data: x j i 2 N, p(•; ↵ kl ) = P(µ k ⌫ l kl ) 3 Continuous data: x j i 2 R, p(•; ↵ kl ) = N (µ kl , 2 kl )
3 The Poisson parameter is here split into µ k and ⌫ l the e↵ects of the row k and the column l respectively and kl the e↵ect of the block kl. Unfortunately, this parameterization is not identifiable. It is therefore not possible to estimate simultaneously µ k , ⌫ l and kl without imposing further constraints. Constraints P k ⇡ k kl = P l ⇢ l kl = 1 and

P k µ k = 1, P l ⌫ l =
MLE: ✓ = arg max ✓ `(✓; x)
Complete log-likelihood:

`c (✓; x, z, w) = lnp(x, z, w ; ✓) = X i,k z ik log ⇡ k + X k,l w jl log ⇢ l + X i,j,k,l
z ik w jl log p(x j i ; ↵ kl )

Be careful with asymptotics. . .

Q(✓, ✓ (q) ) = E [`c (✓; x, z, w)|x; ✓ (q) ] = X i,k p(z i = k|x; ✓ (q) ) | {z } t (q) ik ln ⇡ k + X j,l p(w i = l|x; ✓ (q) ) | {z } s (q) jl ln ⇢ l + X i,j,k,l p(z i = k, w j = l|x; ✓ (q) ) | {z } e (q) ijkl ln p(x ij ; ↵ kl )
M-step of EM (iteration q): classical. For instance, for the Bernoulli case, it gives

⇡ (q+1) k = P i t (q) ik n , ⇢ (q+1) l
= P j s (q) jl d , ↵ (q+1) kl
= P i,j e (q) ijkl x ij P i,j e (

MLE: variational EM (1/2)

Use a general variational result from [Hathaway, 1985] Maximizing `(✓; x) on ✓ is equivalent to maximize ˜c (✓; x, e) on (✓, e)

˜c (✓; x, e) = X i,k t ik ln ⇡ k + X j,l s jl ln ⇢ l + X i,j,
Q(✓, ✓ (q) ) ⇡ X i,k t (q) ik ln ⇡ k + X j,l s (q) jl ln ⇢ l + X i,j,k,l t (q) ik s (
(x ij = 1|z ik = 1) = ⌧ k = L X l=1 ↵ kl ⇢ l Thus marginal distribution of x ij is a mixture (indep. of x ij cond. z ik = 1) 0 @ X j x ij 1 A |z ik = 1 ⇠ B(d, ⌧ k )
Control of error on this partition mixture estimate ẑmix of binomial distributions p(ẑ mix 6 = z ⇤ )  2n exp n 1

8 d h min k6 =k 0 |⌧ k ⌧ k 0 | | {z } overlap io + K (1 min k ⇡ k ) n
We retrieve also consistency for very high dimension with constraint ln(n) = o(d)

  . . . but density estimation quality decreases with d

  variables where all parts are Gaussians S: set of variables useful for clustering U: set of redondant clustering variables, expressed with R ✓ S W : set of variables independent of clustering Trick Variable selection is recasted as a particular variable role 19in competition: m = (S, R, U, W , K ) ! combinatorics Use a backward stepwise algorithm guided by a model selection criterion: d ⇡ 30

  Notationsz i : the cluster of the row i w j : the cluster of the column j (z i , w j ): the block of the element x ij (row i, column j)z = (z 1 , . . . ,zn): partition of individuals in K custers of rows w = (w 1 , . . . , w d ): partition of variables in L clusters of columns (z, w): bi-partition of the whole data set x Both space partitions are respectively denoted by Z and W Restriction All variables are of the same kind (see discussion at the end) 25/69

  Theoretical error decreases when d grows: err theo = ( kµ 2 µ 1 k ⌃ 1 /2)

		Modeling						Estimating							Selecting			Experiments	To go further
					Some alternatives for reducing variance
	Dimension reduction in non-canonical space (PCA-like typically)
	Dimension reduction in the canonical space (variable selection)
	Model parsimony in the initial HD space (constraints on model parameters)
	Empirical error rate with the (true) intra-correlated model worse with d
	Empirical error rate with the (false) intra-independent model better with d!
		4											0.38
		3											0.36
		2											0.34
																					Empirical corr.
	x2	1										err	0.32								Empirical indep. Theoretical
		0											0.3
		-1											0.28
		-2											0.26
		-4 -3	-3	-2	-1	0	x1	1	2	3	4	5	0.24	1	2	3	4	5	d	6	7	8	9	10
																					13/69

  Of course maximizing `(✓; x) or ˜c (✓; x, e) are both intractable Idea: restriction on e to obtain tractability e ijkl = t ik s jl New variables are thus now t = (t ik ) and s = (s jl )As a consequence, it is a maximization of a lower bound of the max. likelihood

	Modeling	Estimating		Selecting	Experiments	To go further
	MLE: variational EM (2/2)
	Approximated E-step
		max ✓	`(✓; x) max ✓,t,s	˜c (✓; x, e)
				40/69

k,l e ijkl ln p(x ij ; ↵ kl ) where e = (e ijkl ), e ijkl 2 {0, 1}, P k,l e ijkl = 1, t ik = P j,l e ijkl , s jl = P i,k e ijkl

  Binary case: marginals seems so simple mixtures![Brault, 14]

	HD clustering	Modeling Modeling Modeling Modeling Modeling	Estimating Estimating Estimating Estimating Estimating	Selecting Selecting Selecting Selecting Selecting	Experiments Experiments Experiments Experiments Experiments	To go further To go further To go further To go further To go further
		Bayesian estimation: pitch Bayesian estimation: limiting degeneracy Block estimation: consistency Block estimation: non asymptotic properties (1/2) Block estimation: non asymptotic properties (2/2)
	[Brault, 14]
		Probability of x ij with no regard to the column membership is Bernoulli
			p
	Interest for avoiding degeneracy is the prior: it acts as a penalization term [Mariadassou and Matias, 12]
		Everything passes by the posterior distribution of ✓ p(✓|x) / p(x|✓) | {z } log-likelihood p(✓) | {z } D K (a,...,a) D L (a,...,a) model dependent prior Typical choices are Dirichlet for ⇡ and ⇢ (with independence between ⇡, ⇢, ↵) p(✓) = p(⇡) | {z } ⇥ p(⇢) | {z } ⇥ p(↵) | {z } ✓ n,d!1 ! ✓ ⇤ | {z } ) p(ẑ = z n,d!1 ! 1 | {z } we have seen that. . . exact bi-partition retrieval!
		Then, take (for instance) the MAP as a ✓ estimate (use a VEM like algo. . . ) The Dirichlet distribution is conjugate, thus easy calculus
	✓ = arg max ✓ Control degeneracy frequency with the a value: p(✓|x) a = 1: uniform prior, so ✓ is strictly the MLE (no regularisation) Thus we retrieve the HD clustering blessing. . .
		a = 1/2: Je↵reys prior, classical (no informative prior) but may favor degeneracy
		a = 4: a rule of thumb working well for limiting degeneracy frequency
			46/69 51/69	52/69

⇤ , ŵ = w ⇤ |x; ✓)

  Mixture of 1033 medical summaries and 1398 aeronautics summaries Lines: 2431 documents Columns: present words (except stop), thus 9275 unique words Data matrix: cross counting document⇥words ICL criterion: overcome both di culties ICL uses complete likelihood thus no intractability ICL = ln p(x, ẑ, ŵ) = ln p(x|ẑ, ŵ) + ln p(ẑ) + ln p(ŵ)Multinomial case (r levels):[Keribin et al., 2014]

	HD clustering	Modeling Modeling Modeling Modeling	Estimating Estimating Estimating Estimating					Selecting Selecting Selecting Selecting	Experiments Experiments Experiments Experiments	To go further To go further To go further To go further
		Illustration: document clustering (2/2) Outline Outline
		HD clustering HD clustering
		Modeling Modeling
		Poisson model
		Estimating Derive an exact (non-asymptotic) ICL version Deduce an asymptotic approximation of ICL Estimating
	Results with 2⇥2 blocks Selecting ICLbic = `c ( ✓; x, ẑ, ŵ) Selecting	K	2	1	ln(n)	L 1 2	ln(d)	KL(r 2	1)	ln(nd)
		Medline We can make a conjecture for the general case Cranfield
		Medline Cranfield ICLbic = `c ( ✓; x, ẑ, ŵ) Experiments Experiments	1033 0 K 1 2 ln(n)	0 1398 L 1 2 ln(d)	KL⌫↵ kl 2	ln(nd)
		To go further To go further
	Experiment illustrates previous theory: HD clustering is blessing
		54/69 56/69 68/69					53/69		55/69		59/69

Illustration: document clustering (1/2)

C. Biernacki, J. Jacques, C. Keribin (2023). A Survey on Model-Based Co-Clustering: High Dimension and

m = (K , L) typically, but not restricted to 57/69

[Robert, 2017] Algorithm Bi-KM1

We have already seen the SEM algorithm in Lesson 3 (thus we do not detail more) It limits dependency to starting point, so it limits local maxima The S-step: a draw (z (q) , w (q) ) ⇠ p(z, w |x; ✓ (q) ) instead an expectation But it is still intractable, thus use a Gibbs algorithm to approx. this draw

Approximated S-step

Two easy draws z (q) ⇠ p(z|w (q 1) , x; ✓ (q) ) and w (q) ⇠ p(w |z (q) , x; ✓ (q) ) Rigorously speaking, many draws within the S-step should be performed Indeed, Gibbs has to reach a stochastic convergence

In practice it works well while saving computation time

Block estimation: evaluation

Empirical error rate between blocks:

Rand index between blocks: it exists also a recent definition. . .

50/69

BIC criterion: two di culties

Di cult 1: which BIC definition because of the double asymptotic on n and d?

Di cult 2: the observed log-likelihood value is intractable

Could be estimated by harmonic mean but time consuming and high variance 58/69

We can obtain a BIC expression from ICLbic BIC = ICLbic ln p(ẑ, ŵ|x; ✓) et al., 2017] establish that asymptotically on n and d "`( ✓; x) = `c ( ✓; x, ẑ, ŵ ) 00

Thus, since BIC is consistent, ICL is also consistent 48 percentages that a given word appears in an e-mail ("make", "you'. . . ) 6 percentages that a given char appears in an e-mail (";", "$". . . )

Transformation of continuous descriptors into binary descriptors


Lopez Cataixa
email: clopezba@ull.edu.es

Freitas Rui
email: ruifreitas@docente.unicv.edu.cv

Sara S Ratao
email: ratao.sara@gmail.com

Brehmer Patrice
email: patrice.brehmer@ird.fr

Report of a Zoanthus Zone from the Cabo Verde Islands (Central Eastern Atlantic)

Keywords: Marine biodiversity, Maio Island, Molecular analyses, West Africa, Zoanthus sociatus, Zoanthus solanderi

Zoanthus is a genus of colonial cnidarians very common in shallow subtropical and tropical rocky and coral reef ecosystems around the world. Some species, such as Zoanthus solanderi and Z. sociatus, can coexist and cover extensive areas in the subtidal zone of some areas in the Caribbean. We report on a Zoanthus zone' in Western Africa in a bay in the north of Maio Island in the Cabo Verde Archipelago, discovered in October 2014. Molecular analyses showed that the Maio Island carpet is formed by Z. sociatus and Z. solanderi. Similar aggregations have been reported from the other side of the Atlantic Ocean in the Caribbean. This work demonstrates the presence of the Zoanthus zone' in the Central East Atlantic and adds two new records of Zoanthus species in the Cabo Verde Islands.

INTRODUCTION

Zoantharians (Anthozoa: Hexacorallia) are an order of colonial cnidarians with species found in most marine environments from temperate to tropical areas and from the intertidal zone to the deep sea [START_REF] Ryland | Revision of methods for separating species of Protopalythoa (Hexacorallia: Zoanthidea) in the tropical West Pacific[END_REF]. In particular, the zooxanthellate genera Zoanthus and Palythoa are common in subtropical and tropical shallow rocky and coral reef ecosystems [START_REF] Karlson | Reproductive patterns in Zoanthus sp. from Discovery Bay, Jamaica[END_REF][START_REF] Irei | Distribution patterns of five zoanthid species in Okinawa Island. Japan[END_REF]) usually as discrete colonies formed by hundreds of polyps. However, some zoanthids species, such as P. caribaeorum, are frequently found covering large areas in the subtidal zones [START_REF] Sebens | Intertidal distribution of zoanthids on the Caribbean coast of Panama: effects of predation and desiccation[END_REF][START_REF] Acosta | Disease in Zoanthids: dynamics in space and time[END_REF]. In 1981 Karlson and collaborators discovered, on the north coast of Jamaica, the unique place where it is known that local populations of Zoanthus solanderi (LeSueur 1818 sensu [START_REF] Karlson | Alternative competitive strategies in a periodically disturbed habitat[END_REF]) and Z. sociatus (Ellis 1768) can cover a huge portion of the intertidal zone. The area occupied by these two species is so huge that this portion of the intertidal has been called the 'Zoanthus zone' [START_REF] Karlson | Reproductive patterns in Zoanthus sp. from Discovery Bay, Jamaica[END_REF]. In October 2014 a huge Zoanthus aggregation was accidentally discovered in Porto Cais Bay, one of the northernmost and sheltered areas of Maio Island, Cabo Verde. This bay is characterized by being a flat area of warm shallow waters affected by semi-diurnal microtidal low current regimes with relatively low transparency due to terrigenous material accumulation in the water column from leaching. Porto Cais Bay is included in the protected Nature Park of the North of Maio Island (PNNM). This Nature Park is part of a network of protected areas implemented in 2003 with a significant expansion of marine area and final boundaries including a marine area were set in 2014 (Anonymous, 2014). . It is known that several coral species, including zoanthids [START_REF] Pax | Zoantharien aus Französisch Westafrika[END_REF][START_REF] Reimer | Zoanthids of the Cape Verde Islands and their symbionts: previously unexamined diversity in the Northeastern Atlantic[END_REF], have an amphi-atlantic distribution [START_REF] Nunes | Long distance dispersal and connectivity in amphi-Atlantic corals at regional and basin scales[END_REF], however, Z. solanderi and Z. sociatus have never been found in Cabo Verde. In fact, only one species of Zoanthus has been reported from the Archipelago until now, Z. aff. pulchellus (Duchassaing and Michelotti 1864) [START_REF] Reimer | Zoanthids of the Cape Verde Islands and their symbionts: previously unexamined diversity in the Northeastern Atlantic[END_REF], and there is no record that this species has the ability to grow disproportionately. On the other hand, Z. aff. pulchellus is characterized by having polyps partially or completely embedded in a well-developed coenenchyme [START_REF] Reimer | Zoanthids of the Cape Verde Islands and their symbionts: previously unexamined diversity in the Northeastern Atlantic[END_REF]) and the species located in Maio had liberate polyps.

In this study we report on this large aggregation showing that the Zoanthus zone in Jamaica is not an isolated event and some zoanthids species can overspread large areas. Moreover, we utilize both morphological and molecular analyses in order to identify the species of Zoanthus. Our aim is to demonstrate if the ability of cover huge areas is characteristic of Z. solanderi and Z. sociatus and therefore, we registered two new records of Zoanthus species in Cabo Verde, or if on the contrary, this massive grow is common in other Zoantharian species.

MATERIAL AND METHODS

Study site, specimen collection and initial in situ identification

This study was carried out around Maio Island, situated on the largest insular shelf in Cabo Verde Archipelago. The Zoanthus aggregation was discovered in Porto Cais Bay (15°18.95'N, 23°11.06'W) (Figure 1) at around 1.5 m depth in October 2014. For an accurate outer boundary of the Zoanthus aggregation a thorough survey was carried out in October 2016 by snorkeling around the outer edge while obtaining a total of 31 geo-referenced points (in UTM WGS 84) taken with a waterproof GPS receiver. The QGIS (2016) open source software was used to estimate the area as a polygon of the 'Zoanthus zone' in Porto Cais Bay, Maio Island based on geo-referenced points.

Species were initially identified as Zoanthus spp. due to the lack of sandy or muddy encrustations in the body wall. However, the specimens had polyps well clear and free of the coenenchyme, indicating that the component species of the aggregation were not Z. aff. pulchellus (Figure 2). Underwater digital images taken in situ were used for subsequent gross morphological analyses including polyp shape, numbers of tentacles, coenenchyme development, and colour (n = 30 polyps analyzed). While there are historical records of zoantharians in the East Atlantic (e.g., [START_REF] Van Der Land | Marine biota of the Cape Verde Islands[END_REF][START_REF] Morri | Shallow epibenthic communities of Ilha do Sal (Cape Verde Archipelago, eastern Atlantic)[END_REF], many species of Zoanthus remain hard to identify in the field due to high levels of intraspecific variation resulting from phenotypic plasticity of polyp shapes, colours, sizes, and oral disk colours [START_REF] Reimer | Reconsidering Zoanthus spp. diversity: molecular evidence of conspecifity within four previously presumed species[END_REF][START_REF] Ong | Morphologically plastic responses to shading in the zoanthids Zoanthus sansibaricus and Palythoa tuberculosa[END_REF]. Therefore, we identified specimens from Maio Island by utilizing a DNA barcoding approach.

DNA extraction and PCR amplification

For molecular analyses, polyps from both Zoanthus morphotypes (n = four colonies) were collected (Figure 2). DNA was extracted from 50mg of body tissue following the manufacturer's protocol of a DNeasy Blood & Tissue Kit (Qiagen. Tokyo). The nuclear internal transcribed spacer region of ribosomal DNA (ITS-rDNA) was amplified by polymerase chain reaction (PCR) using primers designed by [START_REF] Swain | Phylogeny-based species delimitations and the evolution of host associations in symbiotic zoanthids (Anthozoa, Zoanthidea) of the wider Caribbean region[END_REF]. Mitochondrial cytochrome oxidase subunit 1 (COI) was amplified using newly designed zoantharian-specific primers (COIZPF 5´-GAACTATGTTGGGGGACGAC 3', COIZPR 5`-GGCGGGGTCAAAGAAAGTAG -3') based on a sequence obtained with the general COI primers of [START_REF] Folmer | DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates[END_REF]. Products obtained by PCR were checked by 1% agarose gel electrophoresis and positive reactions sequenced in both directions at Fasmac (Kanagawa, Japan).

Phylogenetic analyses

DNA sequences generated in this study were aligned with previously reported sequences in GenBank [START_REF] Altschul | Basic local alignment search tool[END_REF]) using ClustalW (Thompson et al. 1994)

RESULTS

The zoantharian area was approximately estimated at 4 hectares with a perimeter of approximately 800 meters (Figure 1). The underwater video transect of the carpet (video available at https://youtu.be/CduwshatD60) revealed that below the colonies, the benthos was dominated by typical local stony coral species (i.e., Siderastrea radians, Porites spp., Millepora alcicornis)

attached to isolated pavements [START_REF] Moses | Pavements of Siderastrea radians on Cape Verde reefs[END_REF].

We observed two morphotypes within the 'Zoanthus zone'. All Zoanthus polyps were 'liberae', i.e., well free and clear of coenenchyme [START_REF] Pax | Studien an westindischen Actinien[END_REF], connected by a poorly developed coenenchyme. This grow form differ to the Zoanthus species already registered in Cabo Verde, Z.

aff. pulchellus, whose polyps are embedded in a well-developed coenenchyme [START_REF] Reimer | Zoanthids of the Cape Verde Islands and their symbionts: previously unexamined diversity in the Northeastern Atlantic[END_REF].

The number of tentacles varied from 36 to 58 and colour patterns of an oral disk varied from green and fluorescent green to yellow. According to previous studies these characteristics are similar to both Z. sociatus and Z. solanderi found in the Caribbean Sea [START_REF] Sebens | Intertidal distribution of zoanthids on the Caribbean coast of Panama: effects of predation and desiccation[END_REF][START_REF] Reimer | Species diversity of shallow water zoanthids (Cnidaria: Anthozoa: Hexacorallia) in Florida[END_REF].

The two morphotypes could be clearly distinguished by the sizes of the polyps. One morphotype had polyps with an average oral disc diameter of 3.2 ± 0.6 mm (n = 6 polyps examined), while the other morphotype was smaller 2.64 ± 0.53 mm (n = 11 polyps). The larger morphotype had white striped vertical markings on the body walls of polyps, similar to as previously reported for Z. solanderi [START_REF] Karlson | Alternative competitive strategies in a periodically disturbed habitat[END_REF], and the oral disk varied in colour (green, brown, grey) with white dots present. The smaller morphotype matched well with Z. sociatus, as this species has many colour morphotypes. We observed a grey oral disk with yellow or green fluorescent in the middle (Figure 2). The sizes of each morphotype also fit well within values previously reported for each species.

Both ITS-rDNA and COI sequences from the specimens each showed two different genotypes, one for each morphotype. ITS-rDNA sequences from specimens morphologically identified as Z. solanderi (ID 91, 92 and 93 on the phylogenetic tree) were identical to previously reported sequences of Z. gigantus from Japan (DQ442412, DQ442465), the Pacific sister species of Atlantic Z. solanderi [START_REF] Reimer | Species diversity of shallow water zoanthids (Cnidaria: Anthozoa: Hexacorallia) in Florida[END_REF]. COI sequences from these specimens were also similarly identical to Z. gigantus from Japan (AB252675, AB252673). Based on combined results, these three specimens were identified as Z. solanderi. Both ITS-rDNA and COI sequences from the other morphotype (specimen ID94 on the phylogenetic tree) grouped with previously reported Z. sociatus sequences (JX119131, JX119133, KT454373, KT454372) from Florida and Brazil, and thus this specimen was identified as Z. sociatus.

Similar Zoanthus aggregations were also subsequently observed in two more locations on more northerly Boavista Island in Cabo Verde (i.e., Derrubado and Gatas Bays; Carlos ML Afonso Pers. Comm.). However, field surveys of other sheltered bays off the western coast of Maio Island (Ponta Rica, Calhetinha, and two bays in Pau Seco) in October 2016 revealed no Zoanthus spp.

(triangular arrows in Figure 1).

DISCUSSION

Our results showed that the Zoanthus aggregation in Porto Cais Bay, Maio Island was composed of two species: Z. sociatus and Z. solanderi. It has been previously observed that these two species can reach high densities [START_REF] Karlson | Reproductive patterns in Zoanthus sp. from Discovery Bay, Jamaica[END_REF] while covering extensive areas of reef flats in the Caribbean Sea, and such areas have been called the 'Zoanthus zone' [START_REF] Sebens | Autotrophic and heterotrophic nutrition of coral reefs zoanthids[END_REF]. This work reports for the first time the presence of such a 'Zoanthus zone' in the central Eastern Atlantic.

While it is not known in this 'Zoanthus zone' is a recent occurrence or not, it is known that Zoanthus spp. can asexually reproduce [START_REF] Burnett | Patterns of genetic subdivision in populations of a clonal cnidarian, Zoanthus coppingeri, from the Great Barrier Reef[END_REF], and that they can expand their benthic coverage rapidly after disturbances [START_REF] Wee | Zoantharian abundance in coral reef benthic communities at Terengganu Islands, Malaysia[END_REF]. Polyclonal Zoanthus aggregations are theorized to be beneficial to success of sexual reproduction [START_REF] Karlson | Growth and survivorship of clonal fragments in Zoanthus solanderi Lesueur[END_REF] Future investigations should try to elucidate if such 'Zoanthus zones' are native to Cabo Verde or rather recently arrived due to an anthropogenic effect (e.g., invasive species, phase shift [START_REF] Cruz | Evidence of a phase shift to Epizoanthus gabrieli Carlgreen, 1951 (Order Zoanthidea) and loss of coral cover on reefs in the Southwest Atlantic[END_REF]).

  implemented in MEGA v.7 (Kumar, Stecher and Tamura 2015). The alignments were cleaned by eye and ambiguous sites were removed prior to analyses. A suitable outgroup to each alignment was added to the phylogenetic analyses. The best-fitting DNA model of nucleotide substitution was determined for each data set according to the Bayesian Information Criterion (Schwarz 1978) in Jmodeltest (Darriba et al. 2012). Phylogenetic trees were inferred using MrBayes (Huelsenbeck and Ronquist 2001) and RAxML (Ott et al. 2007) software through the Mobyle SNAP Workbench (Monacell and Carbone 2014). Posterior distributions of trees and parameters were estimated using Markov Chain Monte Carlo (MCMCs) for 5000000 generations sampled every 100 generations. The consensus trees were visualized with Fig.Tree v1.4.2 software (Rambaut 2009).

Fig. 1 :

 1 Fig. 1: Location of the Cabo Verde Archipelago in the eastern Atlantic Ocean and the large

Fig. 2 :

 2 Fig. 2: Underwater picture of Zoanthus sociatus and Z. solanderi observed inside the 'Zoanthus

  . Although the Pau Seco area is similar to Porto Cais Bay, both Pau Seco bays are more hydrodynamic (current and wind intensities) with steeper slopes and a more complex topography of the substrata (sand and rubble) in comparison to the 'Zoanthus zone' location found at Porto Cais Bay; these facts may explain the absence of Zoanthus spp. in Pau Seco. Detailed exploration around the Cabo Verde Archipelago and neighbouring islands is needed to confirm if the Zoanthus zone is a common occurrence in this region or a relatively unique ecosystem. However, until now, no such zone has been reported in any of the Cabo Verde Islands or along the West African coasts, and thus it is likely that this is either a relatively rarely occurring ecosystem and/or has appeared relatively recently in this part of the Atlantic Ocean. Long-term monitoring of the site at Porto Cais Bay is needed to understand if the 'Zoanthus zone' is expanding or not, and to elucidate the long-term dynamics of this phenomenon.

ACKNOWLEDGEMENTS

We thank the University of Cabo Verde for support in the prospecting survey and Nivaldo Ramos (FMB) for local assistance and field equipment. We are thankful the FMB local NGO for providing the opportunity to undertake this study. We are indebted to Mariano Hernández for his kindly help in primers designs and Carlos M. L. Afonso (UAlg, CCMar) for shared information on Zoanthus colonies. The Agencia Canaria de Investigación (ACIISI) and Fondo Social Europeo supported CL's internship in Okinawa, Japan where analyses were carried out. The senior author was

supported by JSPS Kiban B grant entitled "Global evolution of Brachycnemina and their Symbiodinium".

FINANCIAL SUPPORT

This research received no specific grant from any funding agency, commercial or not-for-profit sectors.

CONFLICT OF INTEREST

On behalf of all authors, the corresponding author states that there is no conflict of interest.