Keywords: COPD assessment test AI, artificial intelligence hounsfield unit SMOTE, synthetic minority oversampling technique DLCO, diffusion capacity for carbon monoxide KCO, carbon monoxide transfer coefficient

RESUME

Introduction

Le phénotype emphysème peut bénéficier d'un traitement spécifique, éventuellement préventif de sa progression. Sa reconnaissance précoce est donc cruciale surtout lorsque les accès à l'imagerie quantitative et/ou à la consultation pneumologique sont limités. Au cours de l'étude AIR mené chez 614 patients, 71 % des sujets étaient classés GOLD 1 et 2 dont les 2/3 présentaient un emphysème majoritairement méconnu. Nous avons souhaité vérifier s'il existait des clusters individualisables sur une partie de cette cohorte, en incluant des mesures d'imagerie quantitatives de l'emphysème, des données cliniques et spirométriques simples, sans mesure de la distension thoracique ni de la DLCO. Nous avons ensuite testé divers modèles de prédiction du phénotype emphysémateux à partir de toutes les variables, en excluant les mesures d'imagerie, comme cela a été précédemment démontré par une étude post-hoc sur une cohorte de BPCO plus sévères, dans COSYCONET.

Méthodes

Notre étude a concerné 143 des 614 participants du centre Niçois de l'essai AIR. Les données cliniques, fonctionnelles et TDM ont été colligées sur un clone de la consultation digitale Colibri-BPCO. La méthode de complétion basée sur les forêts aléatoires missForest a permis de compléter 9,4% de données manquantes. Pour le clustering nous avons utilisé une technique d'apprentissage non supervisée, l'Analyse Factorielle de Données Mixtes (AFDM), puis choisi le nombre idéal de clusters basé sur un algorithme de matrices de dissimilarité COSA (Clustering Objects on Subsets of Attributes). Enfin, pour l'analyse de prédiction, nous avons employé diverses techniques de forêts aléatoires dont HyperSMURF.

Résultats

Parmi nos 143 patients, il s'agissait majoritairement d'hommes (66,4%), âgés de 62,2 ± 6,9 ans, dont la plupart (57,3%) étaient fumeurs actifs à 57,2 ± 22,7 Paquets-Années. Des images L'algorithme prédictif de l'emphysème retrouvait une Exactitude (Accuracy) de 0,9643 (95% CI 0,88) à l'aide des variables âge, IMC, VEMS/CVF, grade et groupe GOLD (VEMS, Exacerbations dans l'année, CAT score) et DIRECT.

Conclusion

Le phénotype emphysème pourrait être prédit avec une quasi-certitude chez des sujets fumeurs avec une BPCO légère sans scanner thoracique, ni mesures de la distension thoracique et de la DLCO. Ce modèle reste à valider sur d'autres cohortes de sujets fumeurs avec divers niveaux de sévérité de la fonction respiratoire mesurée par spirométrie lente et forcée.

Introduction

Chronic obstructive pulmonary disease (COPD) is the third cause of death worldwide, accounting for 6% of total deaths in 2019 (WHO Global Health Estimates (1)). It is characterized by persistent respiratory symptoms and airflow limitation due to airway and/or alveolar abnormalities usually caused by significant exposure to noxious particles or gases and influenced by host factors [START_REF] Rabe | Global Strategy for the Diagnosis, Management, and Prevention of Chronic Obstructive Pulmonary Disease: GOLD Executive Summary[END_REF].

Since 2017, the GOLD committee released a classification wherein FEV1 is separated from clinical parameters (exacerbation history, modified Medical Research Council dyspnea mMRC and COPD Assessment Test CAT scores)(3). This classification often guides clinical practice, but lacks statistical support with respect to its predictive ability regarding clinically meaningful outcomes, such as mortality and response to treatment (4).

Therefore, in parallel with the GOLD Guidelines, some authors explored the heterogeneity of the COPD disease. Nonetheless, there is a great disparity in the literature as to the best way to split patients into different phenotypes. The number of COPD phenotypes identified varies in the literature, with the most frequently reported being: Alpha1-Antitrypsin Deficiency COPD (5), upper-lobe emphysema and low base-line exercise capacity [START_REF]A Randomized Trial Comparing Lung-Volume-Reduction Surgery with Medical Therapy for Severe Emphysema[END_REF] and frequent-exacerbator [START_REF] Hurst | Susceptibility to Exacerbation in Chronic Obstructive Pulmonary Disease[END_REF]. Some authors define up to seven phenotypes [START_REF] Mirza | Chronic Obstructive Pulmonary Disease Phenotypes: Implications for Care[END_REF] with specific aspects of the management for each.

Phenotyping enables a precision medicine in order to personalize COPD patients care, with tailored treatments. Recently, research has increasingly called for the integration of other determinants into phenotyping, such as physiological characteristics, comorbidities, pulmonary function tests, tomographic images and biomarkers (4). However, a universal and properly defined classification of phenotypes lacks due to major methodological challenges that remain in studies on COPD phenotyping at various stages : selection of patients, completeness of dataset, statistical analysis and follow-up (9).

To this complex task involving large amounts of data, where the patterns can be somehow cryptic to the human-based reasoning, Artificial Intelligence (AI), through its branch Machine Learning (ML) or Data Mining (DM), came to the rescue [START_REF] Exarchos | Review of Artificial Intelligence Techniques in Chronic Obstructive Lung Disease[END_REF]. To achieve these aims, ML models undergo multiple phases that can be summarized into : 1/ Acquisition of a worthy dataset, 2/ data preparation (including treatment of missing data) [START_REF] Stekhoven | MissForest--non-parametric missing value imputation for mixedtype data[END_REF], 3/ selection of the appropriate ML algorithms, 4/ the training phase where the model is trained to learn the dependencies among the input variables : learning can be either supervised wherein the algorithm is trained using labeled datasets to classify data and predict outcomes accurately, or unsupervised wherein algorithms discover hidden patterns in data without the need for human intervention, 5/ the testing phase where new instances are presented to the aforementioned model in order to assess its generalizability [START_REF] Arbet | Lessons and tips for designing a machine learning study using EHR data[END_REF].

In COPD, clustering studies in diagnosis, patients classification, management and prognosis [START_REF] Exarchos | Review of Artificial Intelligence Techniques in Chronic Obstructive Lung Disease[END_REF] have started emerging recently [START_REF] Weatherall | Use of cluster analysis to define COPD phenotypes[END_REF]. In the COPD phenotyping where AI is commonly used, many authors applied ML to large sets of COPD patients in order to unfold COPD clusters with different characteristics and therefore eventually different personalized treatment. For instance, Pikoula and al. ( 14) applied an unsupervised classification scheme cluster analysis method on more than 30,000 COPD patient in order to identify clinically meaningful patient subgroups using as input 16 clinical features (e.g. risk factors, comorbidities) with no chest imaging. They identity five COPD patient clusters with the following leading characteristics: 1/ anxiety/depression, 2/ severe airflow obstruction and frailty, 3/ cardiovascular disease and diabetes, 4/ obesity/atopy and 5/ low prevalence of comorbid conditions. Other authors analyzed chest computed tomography images of COPD patients with correlating them to the pulmonary function tests (15)( 16) [START_REF] Karabulut | Emphysema discrimination from raw HRCT images by convolutional neural networks[END_REF] with the aim of generating different emphysema phenotypes.

The emphysema phenotype has a specific management encompassing lung volume reduction (18) [START_REF] Criner | A Multicenter Randomized Controlled Trial of Zephyr Endobronchial Valve Treatment in Heterogeneous Emphysema (LIBERATE)[END_REF]. Therefore its detection is very important especially in smokers followed in primary care when a chest CT scan or specialist lung consultation are difficult to access. We therefore retrospectively analyzed a cohort of COPD patients with chest CT, with the aim of potentially identifying multidimensional clinical clusters including categorical CT variables (emphysematous vs. nonemphysematous) and then creating a predictive model of the emphysema phenotype using only standard clinical parameters, questionnaire scores and basic forced spirometry measures.

Methods

Participants

The AIR trial [START_REF] Marquette | Circulating tumour cells as a potential biomarker for lung cancer screening: a prospective cohort study[END_REF] was a multicenter study done in 21 university centers in France, testing whether circulating tumor cells could be used as a biomarker for lung cancer screening in 614 patients with COPD.

For our proof-of-concept study, we included all patients in the AIR trial who were recruited in Nice, France. Those without COPD were excluded. All patients therefore had persistent respiratory symptoms and fixed airflow limitation with a post-bronchodilator FEV1/Forced Vital Capacity ratio of less than 0,7, and the criteria of the National Lung Screening trial for lung cancer screening [START_REF] National | Reduced Lung-Cancer Mortality with Low-Dose Computed Tomographic Screening[END_REF] : aged 55-74 years, with at least 30 pack-year smoking history.

Clinical evaluation

Baseline data included anthropomorphic, clinical and anamnestic data, spirometric data as well as a chest CT scan.

COPD symptoms, quality of life and disability were evaluated using respectively the modified Medical Research Council (mMRC) [START_REF] Bestall | Usefulness of the Medical Research Council (MRC) dyspnoea scale as a measure of disability in patients with chronic obstructive pulmonary disease[END_REF], CAT score and the Disability Related to COPD Tool (DIRECT) [START_REF] Regnault | Disability related to COPD tool (DIRECT): towards an assessment of COPD-related disability in routine practice[END_REF]. COPD's impact on mental health was assessed using Hospital Anxiety and Depression Scale. It is a combined questionnaire (depression and anxiety) with 14 questions.

Spirometry

The lung function parameters used in our study were: Forced Expiratory Volume in one second (FEV1), Forced Vital Capacity (FVC) and their ratio FEV1/FVC, as well as maximal Forced Inspiratory Flow (FIFmax) and inspiratory capacity. Plethysmography and carbon monoxide diffusion capacity were not performed as per the AIR protocol.

CT scan

At baseline, patients underwent Non-Contrast-Enhanced inspiratory CT scan.

Lung volumes were estimated and emphysema was quantified per lobe using a quantitative software, StratX (Thirona's lung quantification software).

Chest CT Phenotypes

We classified patients into four phenotypes, mainly based on the severity of emphysema on chest CT and also taking into account the annual exacerbation rate in the previous year before inclusion.

Patients with more than 15% (24)(25) [START_REF] Ding | Role of the Emphysema Index Combined with the Chronic Obstructive Pulmonary Disease Assessment Test Score in the Evaluation of Chronic Obstructive Pulmonary Disease. Sontake V, éditeur[END_REF] of voxels less than -950 HU at quantitative CT (27)(28) [START_REF] Gevenois | Comparison of computed density and macroscopic morphometry in pulmonary emphysema[END_REF], not having alternative abnormalities such as lung cysts, were defined as having Emphysema. A heterogeneity index was established to differentiate between low and high levels of emphysema heterogeneity. A difference of 15% in tissue destruction was used as a threshold to differentiate between homogenous and heterogenous emphysema [START_REF] Valipour | Patterns of Emphysema Heterogeneity[END_REF].

The exacerbators were defined as patients having at least 2 exacerbations or an hospitalization in the previous year [START_REF] Hurst | Susceptibility to Exacerbation in Chronic Obstructive Pulmonary Disease[END_REF]. Based on these two parameters, four phenotypes were defined:

-Emphysema phenotype : emphysema without exacerbation -Exacerbator phenotype : Exacerbator without emphysema -Mixed phenotype : Exacerbators with emphysema -Mild phenotype : none exacerbator non emphysematous patients.

Patients with Asthma-COPD overlap were not included since reversibility was an exclusion criterion in the AIR trial.

Statistical analysis Managing missing data

Our data set was composed for each patient by the observed values of two types of variables, continuous and categorical and presents several missing values (9,4%). Assuming that these missing values over the variables in the data set occur completely at random (the process of missing entries is not influenced by the data or the data generating process), we impute the missing data by simply handing it over to missForest [START_REF] Stekhoven | MissForest--non-parametric missing value imputation for mixedtype data[END_REF], which is a nonparametric imputation method for basically any kind of data. The missForest algorithm can cope with mixed-type of variables, nonlinear relations, complex interactions and high dimensionality. It only requires the observation (i. e. the rows of the data set supplied to the algorithm) to be pairwise independent. The algorithm is based on random forest [START_REF] Breiman | Using Iterated Bagging to Debias Regressions[END_REF] and is dependent on its R implementation random Forest by Andy Liaw and Matthew Wiener [START_REF] Liaw | Classification and Regression by randomForest[END_REF]. Put simple: for each variable missForest fits a random forest on the observed part and then predicts the missing part. The algorithm continues to repeat these two steps until a stopping criterion is met or the user specified maximum of iterations is reached [START_REF] Stekhoven | MissForest--non-parametric missing value imputation for mixedtype data[END_REF].

Cluster analysis

To enrich the description of the data, we then performed on the imputed data set an unsupervised cluster analysis according to the recorded characteristics of the study participants patients. To achieve this task, we combined three exploratory data analysis methods, a factor analysis of mixed data (FAMD) (a kind of principal component method dedicated to analyze a data set containing both continuous and categorical variables (see Lê, S, Josse, J. and Husson, F. ( 2008)) [START_REF] Lê | FactoMineR: An R Package for Multivariate Analysis[END_REF], hierarchical clustering (see e.g. Murtagh F and Contreras P (2012)) [START_REF] Murtagh | Algorithms for hierarchical clustering: an overview[END_REF] and partitioning (35).

Principal component methods are used as preprocessing step, to clustering. From principal components, which are the linear amalgamation of the original variables, distances between individuals are derived. Using then these distances (compatible with mixed data) we performed a hierarchical clustering method, with an optimal choice of the possible number of clusters, based on a dissimilarity matrix derived by the COSA (Clustering Objects on Subsets of Attributes) algorithm proposed by [START_REF] Friedman | Clustering objects on subsets of attributes (with discussion)[END_REF] [START_REF] Friedman | Clustering objects on subsets of attributes (with discussion)[END_REF]. The identified clusters are characterized by considering between-cluster differences of patients' recordings using appropriate statistical tests adapted to the types of recordings. All analyses were performed using R software. Evaluation of emphysema on CT scan showed that 107 patients of 143 (74,8%) had emphysema.

Cluster generating features

On average, 27,82% (SD 13,4%) of the lung was affected by emphysema. For 47,7% of emphysema patients (51 of 107), emphysema was heterogenous and for 52,9% of them (27 of 51) it is predominantly in upper lobe. Seven patients (4.9%) had chest CT scans that were not analyzable for emphysema quantification using StratX.

When we focus on phenotypes, 98 of 143 patients (68,5%) had emphysema phenotype. The second was mild phenotype (25 of 143 (17,5%)). Only 4 patients had the Exacerbator phenotype (2,8%) and 9 (6,3%) the mixed phenotype.

The Hopkins' statistic value was 0,68 suggesting clustered data (Supplementary appendix 1). The optimal cluster analysis separated our population in two clusters (Figure 1 and2). Since our population was split in two clusters mainly depending on the presence of emphysema, we tested multiple non-tomographic variables in order to try to predict the Emphysema Phenotype in patients without the need for a chest CT. The statistical analysis highlights a very robust model with nine variables (Table 2) which can predict with a high power the Emphysema phenotype (Sensitivity: 100%, Specificity: 100%, Positive predicted value: 100%, Negative predicted value: 100%, Accuracy: 0,9643 with (95% CI 0,8766; 1). In this model, we count two common Pulmonary Function Test (FEV1/FVC, FEV1), two common pulmonary questionnaires (CAT and DIRECT), the GOLD Classification, age and Body Mass Index (BMI).

Discussion

We have demonstrated that using an unsupervised ML-based clustering, we were able to recognize a cluster of emphysematous patients which is concordant with a clinical-based phenotyping. To our surprise, when we eliminated tomographic data from the dataset and we analyzed a predictive model of emphysema phenotype versus non emphysema phenotype, we were still able to predict with an impressive 96 % accuracy the presence of emphysema phenotype using 9 readily available clinical data that we can rearrange into only four:

-Age -BMI: weight / height 2 -GOLD score: FEV1/FVC, FEV1, Exacerbations in the previous year, CAT score -Direct questionnaire.

COPD is now regarded as an umbrella of different sub-types of diseases presenting with respiratory symptoms, fixed obstruction and tomographic changes [START_REF] Stolz | Towards the elimination of chronic obstructive pulmonary disease: a Lancet Commission[END_REF]. Therefore, endophenotyping is becoming the cornerstone of understanding and management of COPD using tailored treatments. Nonetheless, COPD phenotyping is very complex and that's related to the fact that phenotypes are rarely mutually exclusive and often overlap [START_REF] Lopez-Campos | Moving Towards Patient-Centered Medicine for COPD Management: Multidimensional Approaches versus Phenotype-Based Medicine-A Critical View[END_REF]. In their recent commission towards the elimination of COPD, authors introduce a distinction between COPD types (genetic, connected to early life events, post-infectious, smoking and environmental exposure) those articulate with 10 various clinical phenotypes [START_REF] Stolz | Towards the elimination of chronic obstructive pulmonary disease: a Lancet Commission[END_REF]. The current tendency is therefore clearly towards sub-typing and clustering this heterogeneous disease in order to better tailor treatment. It is now clear that AI-based clustering could serve to find homogenous groups of patients deserving specific treatments [START_REF] Pikoula | Identifying clinically important COPD sub-types using data-driven approaches in primary care population based electronic health records[END_REF][START_REF] Peng | Classification and Quantification of Emphysema Using a Multi-Scale Residual Network[END_REF][START_REF] Sorensen | Quantitative Analysis of Pulmonary Emphysema Using Local Binary Patterns[END_REF][START_REF] Karabulut | Emphysema discrimination from raw HRCT images by convolutional neural networks[END_REF]. Nonetheless, until now, clustering COPD patients has not allowed for a significant modification in our understanding of phenotyping or treatment. As long as we are unable to universally phenotype patients, and to understand the various COPD endo-phenotypes pathophysiologies, it will continue to be hard to develop phenotype-specific treatments. Emphysema phenotype is one that is clearly commonly accepted [START_REF] Burgel | Clinical COPD phenotypes: a novel approach using principal component and cluster analyses[END_REF] and specific care is already available (lung volume reduction [START_REF] Herth | Endoscopic Lung Volume Reduction: An Expert Panel Recommendation -Update 2019[END_REF]) or being tested (potential role of Metformin in Emphysema treatment ( 43)). It is therefore important to be able to identify this phenotype even in the absence of chest CT scan.

Contrary to most clustering studies of COPD patients which usually report at least 3 clusters (4), we chose to see, based on what seemed to be the best number of clusters, that is two clusters, if we could isolate the emphysema phenotype. This was the case. Clustering distinguished emphysematous from non-emphysematous patients when tomographic data were included in the dataset. These two clusters are clinically distinguishable, which is not always obvious in clusters found in the literature due to the mathematical construction of clustering method (4). What was surprising, is that our model was still able to distinguish these classes without any tomographic data. Besides emphysema, we found other significant differences between our two clusters such as cardiovascular comorbidities (46,2% in cluster 1 versus 23,1% in cluster 2; p-value = 0,031)

and the HAD depression score (5,2 versus 3,7; p-value = 0,046). It should be noted that cardiovascular comorbidities seem usually to be more associated with chronic bronchitis phenotype [START_REF] Wedzicha | Mechanisms and impact of the frequent exacerbator phenotype in chronic obstructive pulmonary disease[END_REF][START_REF] Domenech | High risk of subclinical atherosclerosis in COPD exacerbator phenotype[END_REF]. In our clusters, it is more frequent in the "emphysema cluster" but that is very probably because the other cluster is mainly comprised of patients with a mild phenotype and not exacerbators. It is also remarkable that the two clusters do not seem to be distinguishable to the basic reasoning since there is no statistical difference between the two groups as per the prediction variables used in the model.

Moreover, contrary to some clustering studies, we decide to not select a priori variables to clusterise. For instance, Chubachi and al [START_REF] Chubachi | Identification of five clusters of comorbidities in a longitudinal Japanese chronic obstructive pulmonary disease cohort[END_REF] One could wonder if adding plethysmographic and diffusion data to this prediction model could have provided more precise prediction model with information on homogeneity and heterogeneity of emphysema or its severity. This needs to be tested in a population with these spirometric data, which was not the case of the AIR trial cohort.

Our study has some limitations. First, it is monocentric, and retrospective and therefore suffered from a certain amount of missing data (9,4%), managed by a non-parametric imputation method (recommended in clustering analysis ( 4)). Our definition of tomographic emphysema is not consensual. We chose a threshold of fifteen percent which is the most common in literature [START_REF] Lim | Mixed Phenotype of Emphysema and Airway Wall Thickening Is Associated with Frequent Exacerbation in Chronic Obstructive Pulmonary Disease Patients[END_REF][START_REF] Loh | A novel CT-emphysema index/FEV1 approach of phenotyping COPD to predict mortality[END_REF].

For instance, some authors have used the cutoff of 6 % of voxels of the parenchyma with < -910

Hounsfield Units [START_REF] Lynch | CT-Definable Subtypes of Chronic Obstructive Pulmonary Disease: A Statement of the Fleischner Society[END_REF].

At this stage, an external validation study would allow to test the generalizability of our conclusions especially to more severe COPD patients and to exacerbators which were a minority in our population.

Conclusion

To conclude, the present study demonstrates that emphysema as a COPD phenotype could be predicted without chest CT. This is of significant importance to physicians with limited technical resources regarding respiratory diseases [START_REF] Stolz | Towards the elimination of chronic obstructive pulmonary disease: a Lancet Commission[END_REF]. Even when a chest CT is performed, emphysema quantification is not readily available on reports, and to an unexperienced eye, it can be overlooked.

Our study brings therefore new data in phenotyping COPD patient with clustering method and offers a new tool in order to predict with high accuracy emphysema phenotype without chest CT.

The D index is a graphical method of determining the number of clusters.

In the plot of D index, we seek a significant knee (the significant peak in D index second differences plot) that corresponds to a significant increase of the value of the measure. Among all indices: 9 proposed 2 as the best number of clusters 7 proposed 3 as the best number of clusters 2 proposed 4 as the best number of clusters 2 proposed 13 as the best number of clusters 3 proposed 14 as the best number of clusters 1 proposed 15 as the best number of clusters Conclusion: According to the majority rule, the best number of clusters is 2

7 ( 3 ,

 73 d'emphysème étaient retrouvées à des degrés variables chez 74,8 % des sujets. Deux clusters ont été isolés à partir des variables anamnestiques, cliniques, spirométriques, scanographiques et de questionnaires HAD, CAT et DIRECT. Dans le cluster 1, 97,4% des individus présentaient des lésions d'emphysème tandis que le cluster 2 ne renfermait pas d'emphysémateux. Seuls les comorbidités cardiovasculaires (46,2% vs 23,1% ; p= 0.031) et le score HAD-D (5,2 (3,7%) vs 3,3%) ; p=0.045) différenciaient significativement les clusters 1 et 2.

  Categorical variables: Sex, cardiovascular comorbidity, diabetes comorbidity, psychiatric comorbidity, cancer comorbidity, obstructive sleep apnea syndrome comorbidity, chronic respiratory failure comorbidity, work status, smoking status, pneumococcal vaccine status, flu vaccine status, classification of daily activities, treatment with inhaled corticosteroids, long acting bronchodilator inhalers, long acting muscarinic antagonist, frequent cough, frequent expectoration, , exacerbation severity, Gold grade (1,2,3,4), Gold group (A,B,C,D), right, left and total lung emphysema heterogeneity. Continuous variables: Age, smoking pack-years, duration of smoking in years, weight, height, BMI, leukocyte, neutrophil, eosinophil counts, FEV1, FVC, FEV1/FVC, FIFmax, slow vital capacity, inspiratory capacity, inspiratory reserve volume, peak expiratory flow, right upper lobe emphysema -910 UH, middle lobe emphysema -910 UH, right lower lobe emphysema -910 UH, left upper lobe emphysema -910 UH, left lower lobe emphysema -910 UH, right upper lobe emphysema -950 UH, middle lobe emphysema -950 UH, right lower lobe emphysema -950 UH, left upper lobe emphysema -950 UH, left lower lobe emphysema -950 UH, inspiratory volume in right upper lobe, inspiratory volume in middle lobe, inspiratory volume in right lower lobe, inspiratory volume in left upper lobe, inspiratory volume in left lower lobe, inspiratory volume in right lung, inspiratory volume in left lung, inspiratory volume in total lung, mMRC score, CAT score, DIRECT score, Epworth score, HAD-A, HAD-D and HAD-A+D score.

Figure 1 :

 1 Figure 1: Cluster with COSA method (most robust and significative method for our database)

Figure 2 :

 2 Figure 2: Cluster in the first two components

Table 2 : Emphysema Phenotype Prediction model

 2

  decided to use only comorbidities data while Yoon and al[START_REF] Yoon | Prediction of first acute exacerbation using COPD subtypes identified by cluster analysis[END_REF] selected a combination of lung functions, comorbidities and demographic data. Our clustering method is robust due to the use of three exploratory data analysis methods (FAMD, Hierarchical and partionning clustering).Concerning our predictive Emphysema Phenotype model, nine non-tomographic parameters (two common Pulmonary Function Tests, two common pulmonary questionnaires, Gold classification, exacerbation evaluation, and some basic information such as age and BMI) enable to predict if a patient belonged to the emphysema phenotype or not with a high degree of certainty. This finding could have a significant impact on COPD management and can help phenotyping patients based on clinical and basic spiromectric data only. Indeed, none of the pulmonary function tests variables reliably evaluating distension or emphysema (residual volume, total lung capacity, Diffusing Capacity for Carbon monoxide (DLCO)) were part of our prediction model. Therefore, the prediction variables we used do not inform about the presence of Emphysema phenotype or not, except for the BMI which is more often low in emphysematous patients[START_REF] Shin | Association of body mass index and COPD exacerbation among patients with chronic bronchitis[END_REF][START_REF] Ji | Differences in the Outcome of Patients with COPD according to Body Mass Index[END_REF].

	Kellerer and al (50) did a similar trial and found, in a COPD cohort of 436 patients
	(COSYCONET), variables which could predict emphysema in COPD patients (specific questions
	of the CAT score and St George's respiratory Questionnaire, mMRC score and lung function
	respiratory variables including diffusion of carbon monoxide). Emphysema phenotype was only
	qualitatively assessed and was not measured using percentage of attenuation. Their sensitivity for
	predicting emphysema was 61 % when combining questionnaires, spirometry and KCO. Our
	model seems to better predict emphysema without the need for diffusion evaluation and by
	including DIRECT score.

The population of our study are COPD patients recruited in an outpatient setting to perform lung cancer screening. Some of them were not known to have COPD before the AIR trial. Therefore, on average they are slightly symptomatic and have overall mild disease. Our results are applicable to the patients recruited in only one center, Nice, France. It would be interesting as a start, to try to externally validate the model on the rest of the cohort (471 patients). In case the model keeps a good accuracy of prediction, a validation on a more severe COPD patient population could be attempted.

SERMENT D'HIPPOCRATE

Au moment d'être admis à exercer la médecine, je promets et je jure d'être fidèle aux lois de l'honneur et de la probité.

Mon premier souci sera de rétablir, de préserver ou de promouvoir la santé dans tous ses éléments, physiques et mentaux, individuels et sociaux.

Je respecterai toutes les personnes, leur autonomie et leur volonté, sans aucune discrimination selon leur état ou leurs convictions. J'interviendrai pour les protéger si elles sont affaiblies, vulnérables ou menacées dans leur intégrité ou leur dignité.

Même sous la contrainte, je ne ferai pas usage de mes connaissances contre les lois de l'humanité. J'informerai les patients des décisions envisagées, de leurs raisons et de leurs conséquences.

Je ne tromperai jamais leur confiance et n'exploiterai pas le pouvoir hérité des circonstances pour forcer les consciences.

Je donnerai mes soins à l'indigent et à quiconque me les demandera.

Je ne me laisserai pas influencer par la soif du gain ou la recherche de la gloire.

Admis dans l'intimité des personnes, je tairai les secrets qui me seront confiés.

Reçu à l'intérieur des maisons, je respecterai les secrets des foyers et ma conduite ne servira pas à corrompre les moeurs.

Je ferai tout pour soulager les souffrances. Je ne prolongerai pas abusivement les agonies. Je ne provoquerai jamais la mort délibérément.

Je préserverai l'indépendance nécessaire à l'accomplissement de ma mission.

Je n'entreprendrai rien qui dépasse mes compétences. Je les entretiendrai et les perfectionnerai pour assurer au mieux les services qui me seront demandés. J'apporterai mon aide à mes confrères ainsi qu'à leurs familles dans l'adversité.

Que les hommes et mes confrères m'accordent leur estime si je suis fidèle à mes promesses ; que je sois déshonoré et méprisé si j'y manque.

Prediction model

Finally, to address the prediction the Emphysema Phenotype of a patient participating in the study, we first performed a variable selection by feature importance using CATBOOST (a recent gradient boosting algorithm adapted to mixed-type data, see Prokhorenkova and al. ( 2017)) [START_REF] Prokhorenkova | CatBoost: unbiased boosting with categorical features[END_REF] and developed a predictive model using an ensemble regression and classification algorithm which takes into account the unbalance of the Emphysema Phenotype in the data. The final prediction is derived using a recent hyper-ensemble of SMOTE (Synthetic Minority Oversampling Technique) under sampled random forests (HyperSMURF) method, which is based on resampling techniques and a hyper-ensemble approach (see Schubach and al.) [START_REF] Schubach | Imbalance-Aware Machine Learning for Predicting Rare and Common Disease-Associated Non-Coding Variants[END_REF]. All analyses were performed using R software.

Data presentation

Continuous data are presented as mean (± SD) and categorical data as number (percentage).

Results

The health records of 143 COPD patients were retrospectively analyzed (Table 1). Patients were predominantly male (95 of 143 66,4%), aged 62,2 (SD 6,9) years, most of whom (57,3%) were active smokers with an average of 57,2 (SD 22,7) pack-years. The main comorbidities were cardiovascular (60 of 143, 42%) and overweight (BMI 25,8 Kg/m2 SD 5,1, 25 of 143 (35%)).

Our population had a relatively conserved lung function (FEV1pred average 72,6% SD 21,8%) mainly GOLD 1 and 2 (41,3% and 38,5% respectively). Most patients were not frequent exacerbators (21,7% GOLD A and 67,8% GOLD B). Their dyspnea was evaluated as mild to moderate (mMRC scores of 0, 1, 2 accounted for almost 90% of our population) and had a mild disability (DIRECT 8,3 on average SD 6,4). HAD scores were below eleven, the value above which mental status disturbance is highly likely (Anxiety: 5,9 SD 3,7 Depression 4,9 SD 3,7).

APPENDIX Appendix 1: Estimation of the optimal number of clusters according some different technics

The Hubert index is a graphical method of deeming the number of clusters. In the plot of Hubert index, we seek a significant knee that corresponds to a significant increase of the value of the measure i.e the significant peak Hubert index second differences plot


Alana Correia-Martins

Réjean Tremblay

Béatrice Bec

Cécile Roques

Ariane Atteia

Angélique Gobet

Marion Richard

Masami Hamaguchi

Toshihiro Miyajima

Masakazu Hori

Gilles Miron

Stéphane Pouvreau

Franck Lagarde

Failure of bivalve foundation species recruitment related to trophic changes during an extreme heatwave event

Keywords: Climate change, Phenology, Extreme Heat Wave, Bivalves, Pacific Oyster, Crassostrea gigas, Reproduction, Larval Ecology, Cascade of Environmental Effects, Trophic Changes

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

INTRODUCTION

Coastal lagoons provide a wide range of ecosystem services [START_REF] Chapman | Management of coastal lagoons under climate change[END_REF][START_REF] Villamagna | Capacity, pressure, demand, and flow: A conceptual framework for analyzing ecosystem service provision and delivery[END_REF][START_REF] Kermagoret | How does eutrophication impact bundles of ecosystem services in multiple coastal habitats using state-and-transition models[END_REF], associated with biodiversity, including bivalves which are of great ecological interest and high commercial value for some of them. Bivalves also have an important regulatory functions in the ecosystem thanks to their capacity to extract particles, to regenerate and store nutrients and to form hard biogenic structures [START_REF] Smaal | Goods and Services of Marine Bivalves[END_REF]). However, because coastal lagoons are shallow and have limited exchange with the ocean, they are highly sensitive to eutrophication, heat waves, hypoxia and acidification, as well as to the effects of global climate change [START_REF] Lloret | Is coastal lagoon eutrophication likely to be aggravated by global climate change ? 78[END_REF][START_REF] Lu | Major threats of pollution and climate change to global coastal ecosystems and enhanced management for sustainability[END_REF], Thomas et al. 2018). An atmospheric heat wave is defined as five consecutive days with a maximum temperature 5°C above the 1976-2005 normal [START_REF] Jouzel | Scénarios régionalisés: édition 2014 pour la métropole et les régions d'outre-mer[END_REF]. Summer 2019 was characterized by two heat waves of exceptional intensity over France, including the Thau Basin, one from June 24 to July 7, the other from July 21 to 27. The absolute heat record for France (46 °C) was measured in Vérargues in the Hérault administrative department (Météo-France 2019), which includes the Thau basin.

Marine heatwaves (MHW) are extreme events defined as abrupt but prolonged periods of high sea surface temperatures that can occur anywhere, at any time [START_REF] Scannell | Frequency of marine heatwaves in the North Atlantic and North Pacific since 1950[END_REF][START_REF] Schlegel | Nearshore and offshore co-occurrence of marine heatwaves and cold-spells[END_REF][START_REF] Hobday | Categorizing and Naming Marine Heatwaves[END_REF]. High water temperatures increase the metabolic requirements of bivalves [START_REF] Filgueira | Bivalve aquaculture-environment interactions in the context of climate change[END_REF][START_REF] Thomas | Assessing the sensitivity of bivalve populations to global warming using an individual-based modelling approach[END_REF]. Even if temperatures remain within the species' thermal range, high temperatures combined with salinity and/or food variations, can negatively impact the life cycle of bivalves [START_REF] Filgueira | Bivalve aquaculture-environment interactions in the context of climate change[END_REF].

Several studies suggest that global changes are disrupting plankton communities and their nutritional values by affecting the abundance, size and diversity of primary producers [START_REF] Klauschies | Climate change effects on phytoplankton depend on cell size and food web structure[END_REF], Sommer et al. 2012[START_REF] Trombetta | Water temperature drives phytoplankton blooms in coastal waters[END_REF]. Generally, elevated temperatures affect phytoplankton cell size with a shift from larger to smaller species [START_REF] Bec | Phytoplankton seasonal dynamics in a Mediterranean coastal lagoon: Emphasis on the picoeukaryote community[END_REF][START_REF] Trombetta | Water temperature drives phytoplankton blooms in coastal waters[END_REF]. Adult bivalves can assimilate small phytoplanktonic particles [START_REF] Sonier | Picophytoplankton contribution to Mytilus edulis growth in an intensive culture environment[END_REF]. However, the efficiency of the capture is regulated by the morphology of their gills, and is generally low when small particles as picoplankton are present [START_REF] Rosa | Selective Capture and Ingestion of Particles by Suspension-Feeding Bivalve Molluscs: A Review[END_REF]. Larvae feed through a less selective velum [START_REF] Bower | Atlas of anatomy and histology of larvae and early juvenile stages of Japanese scallop Patinopecten yessoensis[END_REF]. Marine phytoplankton species are major producers of long-chain polyunsaturated essential fatty acids (EFA) but are predicted to decrease due to ocean warming [START_REF] Hixson | Climate warming is predicted to reduce omega-3, long-chain, polyunsaturated fatty acid production in phytoplankton[END_REF][START_REF] Colombo | A fundamental dichotomy in long-chain polyunsaturated fatty acid abundance between and within marine and terrestrial ecosystems[END_REF]. The fatty acids docosahexaenoic acid (22:6ω3; DHA), eicosapentaenoic acid (20:5ω3; EPA) and arachidonic acid (AA) are essential for the growth and survival of marine invertebrates, particularly during their metamorphosis from pelagic larvae to benthic juveniles and ultimately, their recruitment success [START_REF] Gagné | Lipid requirements of the scallop Pecten maximus (L.) during larval and post-larval development in relation to addition of Rhodomonas salina in diet[END_REF][START_REF] Bassim | Predicting growth and mortality of bivalve larvae using gene expression and supervised machine learning[END_REF]. Since EFAs are poorly biosynthesized by marine animals, their intake depends on their food [START_REF] Glencross | Exploring the nutritional demand for essential fatty acids by aquaculture species[END_REF][START_REF] Costa | Essential Fatty Acid Assimilation and Synthesis in Larvae of the Bivalve Crassostrea gigas[END_REF]. Thus, both the right size and the right fatty acid composition of larval food are essential for the recruitment success of bivalves.

The aim of this study was to identify the environmental factors and trophic conditions (Supplementary material table1 & table 2) associated with the recruitment failure of the Pacific oyster, Crassostrea gigas, during a heat wave. We compared two contrasted years (2017 no heat wave and 2019 heat wave) in four sites in the Thau lagoon, France (Fig. 1). We hypothesize that the heat waves, characterized by high temperatures and high salinity, have a negative impact on oyster recruitment due to poor larval feeding conditions caused by changes in plankton diversity.

MATERIALS AND METHODS

Experimental design

Oyster recruitment was monitored from July 24 to August 21, 2017, and from July 2 to 29, 2019, at four experimental sites in the Thau lagoon (southern France; Fig 1 .). The average depth of the lagoon is 4 m. The lagoon covers an area of 7 500 ha (19 km x 4.5 km) of which 20% is used for shellfish culture (oysters and mussels). The lagoon is connected to the Mediterranean Sea via a network of channels through Sète Harbor [START_REF] Fiandrino | Spatial patterns in coastal lagoons related to the hydrodynamics of seawater intrusion[END_REF]. Two experimental sites were located inside the shellfish farming areas (Marseillan and Bouzigues) while two others were located outside the shellfish farming areas (Meze and Listel) (Fig 1 .).

Oyster analyses

Three sets of oyster collectors were submerged vertically 2 m below the surface at each of the four study sites in the Thau lagoon to collect young settlers (pediveligers settled on collectors, metamorphosed juveniles, and newly metamorphosed juveniles). The collectors were installed once the oyster's larval supply reached a density of 10 000 larvae/m 3 (VELYGER network 44 ).

The collectors located inside the shellfish culture areas were suspended from existing farming structures. Those outside the area were suspended using a mooring system [START_REF] Lagarde | Recruitment of the Pacific oyster Crassostrea gigas in a shellfish-exploited Mediterranean lagoon : discovery, driving factors and a favorable environmental window[END_REF][START_REF] Lagarde | Duality of trophic supply and hydrodynamic connectivity drives spatial patterns of Pacific oyster recruitment[END_REF]. Each collector was made of 44 white PVC plastic plates (15 cm diameter; surface area of 250 cm 2 ) stacked on a 110 cm tube. Two weeks after their immersion, three plates per collector were harvested [at the top (i.e., the 5 th from the surface), in the middle (the 22 nd ) and at the bottom (the 39 th )] and data were pooled to assess the abundance of young settlers and fatty acid (FA) content (μg larva -1 ). A similar sampling procedure was used four weeks after the collectors were immersed to assess the abundance of juveniles.

The abundance of young settlers and juveniles was assessed on the upper surface of each plate using standard 15 cm 2 sub-units. Depending on the abundance, 3 to 12 sub-units were randomly selected for counting and the resulting replicates were averaged to obtain the total number of individuals per plate. Recruitment was evaluated from the abundance of juveniles and metamorphosis from the ratio of juvenile to young-settler abundances. Size at metamorphosis was estimated by measuring the prodissoconch II (PII) [START_REF] Martel | Prodissoconch morphology, planktonic shell growth, and site at metamorphosis in Dreissena polymorpha[END_REF]. A maximum of 60 spats were removed from each plate sampled after the fourth week after immersion, and placed on a plasticine flange fixed on a microscope blade. Observations were made under the widerange zoom lens of a high-resolution digital microscope Keyence (VHX 2000E, 1 μm resolution, HDR images), and the maximum dorsoventral axis was measured. This measurement corresponds to the distance between the umbo and the most distant part of the clear demarcation formed by a growth line delimiting the PII from the dissoconch shell.

The fatty acid (FA) composition of young settlers was determined using a pool of 77 to 212 individuals per replicate (2-3 replicates per site depending on pediveliger abundances). Samples were preserved in vials with 3 mL of dichloromethane methanol (CH2Cl2:MeOH, 2:1 v:v), closed with a Teflon-lined cap under nitrogen atmosphere and stored at -80 °C until analysis.

Lipids were extracted by grinding in dichloromethane methanol using a modified Folch procedure [START_REF] Parrish | Determination of Total Lipid, Lipid Classes, and Fatty Acids in Aquatic Samples[END_REF]. Fatty acid methyl esters (FAME) were prepared using sulfuric acid and methanol (2:98 v:v) at 100 °C for 10 min and using 19:0 as internal standard [START_REF] Lepage | Improved recovery of fatty acid through direct transesterification without prior extraction or purification[END_REF]. Samples were purified on an activated silica gel with 1 mL of hexane ethyl acetate (v:v) to eliminate free sterols. FAME were analyzed in the full scan mode (ionic range: 50-650 m:z) on a Polaris Q ion trap coupled to a Trace GC Ultra gas chromatograph (Thermo Scientific) equipped with a TriPlus autosampler, a PTV injector and an ITQ900 mass detector (Thermo Scientific). An Omegawax 250 (Supelco) capillary column was used for separation using high purity helium. Xcalibur v.2.1 software (Thermo Scientific) was used for FAME identification and quantification with the standards (Supelco 37 Component FAME Mix and Supelco menhaden oil). Unknown peaks were identified according to their mass spectra with emphasis on FA trophic makers.

Environmental measurements

Environmental factors were measured once a week (supplementary files table 1 and table 2) just after immersion of the collectors until the plates were harvested, i.e., a total of five weeks.

Temperature (°C), salinity (PSU) and dissolved oxygen concentrations (mg L -1 ) were measured at a depth of 1 m and at the bottom of the water column with an Oxi1970i WTW oximeter and an LF 197-S WTW conductivity meter.

Potential food for oysters is expressed as the concentration of total suspended particulate matter varying in size from 0.7 and 20 μm (TPM0.7-20µm, mg L -1 ). It consisted of inorganic (PIM0.7-20µm, mg L -1 ) and organic particulate matter (POM0.7-20µm, mg L -1 ). Once a week, three replicate water samples were collected at a depth of 1 m using a Ruttner Standard Water Sampler (Hydro-Bios Apparatebau) and stored at 4 °C for less than 2 hours before filtration for the measurement of the concentrations (mg mL -1 ) of pico and nano-seston. In 2017, 500-mL subsamples of 1-L samples were used for filtration, while 1-L subsamples of 2-L samples were used in 2019. Water samples were first filtered by gravity through a Nuclepore membrane (20 μm pore size).

Fractioned water samples were then filtered using a vacuum pressure pump (0.3 bar) on pre-weighed (Mettler Toledo XP6 microbalance) pre-combusted (at 500 °C) Whatman 25 mm GF/F filters (0.7 μm pore size). The GF/F filters were rinsed with an isotonic seawater solution of ammonium formate (38 g L -1 distilled water) to eliminate salt deposits and stored in Millipore™ PetriSlide™ containers at -25°C. The filters were dried at 70 °C for 24 h, weighed and the concentration of total particulate matter TPM0.7-20µm was determined. The filters were then combusted at 500 °C for 5 h and reweighed to determine the concentration of particulate inorganic matter (PIM0.7-20µm, mg L -1 ). The concentration of particulate organic matter (POM0.7-20µm, mg L -1 ) was the difference in weight between the dried and the combusted filter. To determine the FA content of the pico-and nano-seston (μg.mg TPM0.7-20µm -1 ), 1-L water samples collected in 2017 and 2-L water samples collected in 2019 were filtered as described above without addition of ammonium formate solution. GF/F filters were stored in 3 ml of CH2Cl2:MeOH (2:1 v:v) under a nitrogen atmosphere in vials with a Teflon-lined cap and stored at -80 °C. The mass of total fatty acids in the seston (MTFA; µg mg -1 POM) and its composition (% fatty acids) were obtained as already described for oysters, with lipid extraction carried out by sonification rather than grinding.

Plankton diversity was collected in 1-L samples in 2017 and in 2-L samples in 2019 collected weekly with a Ruttner Standard Water Sampler (Hydro-Bios Apparatebau) at each sampling site. This sampling strategy provided 40 observations (4 sites x 5 weeks x 2 years).

Phytoplankton was characterized using the standard Utermöhl method NF- EN-152014, 2006 were evaluated in 200 ml seawater samples filtered [START_REF] Bec | Phytoplankton seasonal dynamics in a Mediterranean coastal lagoon: Emphasis on the picoeukaryote community[END_REF][START_REF] Bec | Distribution of picophytoplankton and nanophytoplankton along an anthropogenic eutrophication gradient in French Mediterranean coastal lagoons[END_REF] on Whatman GF/F membranes (0.7 μm pore size) with a vacuum pressure pump (<10 cm Hg). Filters were stored in glass tubes at -20 °C until analysis. To determine the contribution of picophytoplankton (<3 μm), nanophytoplankton (3 to 20 μm) and microphytoplankton (>20 μm), two out of three samples were size-fractioned beforehand by gravity through Nuclepore membranes (3 and 20 μm pore size). Filters were ground in acetone (90%) and extracted at 4 °C for 24 h in the dark.

Pigment contents were measured with a spectrofluorometer (Perkin-Elmer LS50b) [START_REF] Neveux | Spectrofluorometric assay of chlorophylls and phaeopigments using the least squares approximation technique[END_REF] and are expressed in μg chl a L -1 . Concentrations of picocyanobacteria (<1 μm), autotrophic picoeukaryotes (<3 μm), nanophytoplankton (3-20 μm) and bacteria were estimated using a FACSCalibur flow cytometer Becton Dickinson methods [START_REF] Marie | Enumeration and cell cycle analysis of natural populations of marine picoplankton by flow cytometry using the nucleic acid stain SYBR Green I[END_REF][START_REF] Bec | Distribution of picophytoplankton and nanophytoplankton along an anthropogenic eutrophication gradient in French Mediterranean coastal lagoons[END_REF]. Seawater samples (1-ml) were analyzed; abundances are expressed in cells per liter. Total picophytoplankton abundances were assessed by summing picocyanobacteria and photosynthetic picoeukaryote abundances. Fluorescent beads (0.94 μm; 2 and 3 µm, Polysciences) were added to each sample. To measure bacterial abundances, seawater samples were fixed with prefiltered (0.2 μm) buffered formaldehyde (2% final concentration) and stored in liquid nitrogen. The procedure was slightly modified as higher concentrations of fluorochromes (SYBR Green I) were used [START_REF] Bouvy | Plankton communities in the five Iles Eparses (Western Indian Ocean) considered to be pristine ecosystems[END_REF]. The fixed samples were incubated with SYBR Green I (Molecular Probes) at a final concentration of 1/375 at 4 °C for 15 min in the dark. Stained bacterial cells excited at 488 nm were determined according to their side-scattered light and green fluorescence collected using a 530/ 30 nm filter. Fluorescent beads (0.94 μm; Polysciences) were added to each sample.

Protozooplankton (heterotrophic flagellates) abundances were determined using the standard 2006 Utermöhl method NF-EN-152014, and are expressed in cells per liter. Until used for heterotrophic flagellate analysis, 30-ml seawater samples were preserved with 2.5-ml of prefiltered (0.2 μm) formaldehyde and kept at 4 °C in the dark. Before counting, 10 ml subsamples were stained with 4',6-diamidino-2-phenylindole (DAPI) for a final concentration of 2.5 μg ml -1 . Heterotrophic flagellates were counted by size class (2-5 μm, 5-10 μm and >10 μm) under an epifluorescence microscope (Olympus AX70) with UV illumination [START_REF] Sherr | Staining of heterotrophic protists for visualization via epifluorescence microscopy[END_REF].

Territorial competition

Percent cover of tubeworm (Ficopomatus enigmaticus) on plates sampled in the fourth week after immersion (6 plates per site) was estimated to assess territorial competition with oyster juveniles, but only during the 2019 sampling season, as no tubeworms were observed in 2017.

Photographs of each plate were taken with a GoPro HERO4 Silver camera equipped with a macro pro filter (San Mateo, CA, USA) and the % of tubeworms recovered on the plate was estimated using Image-Pro Insight 9.1 software (MediaCybernetics, Rockville, MD, USA).

Statistical analyses

All PERMANOVA analyses were performed with Primer 7 and Permanova+1 (version 7.0.13) software. A two-way PERMANOVA (n perm.: 9999) was conducted using a Euclidian distance matrix to test the effect of year (2 fixed levels) and sampling site (4 fixed levels) on size at metamorphosis, total and essential fatty acid contents in young settlers and on all the environmental variables measured, except the oxygen level, which was added as a third factor (depth) in the analysis. Homogeneity was evaluated using the permutation analysis of multivariate dispersion (PERMDISP) routine. When significant PERMANOVAs were observed, post hoc multiple comparison tests were carried out. Multivariate analyses of total FA composition in young settlers and in seston, including a posteriori pairwise comparison, were done using distance-based permutational multivariate analysis of variance (PERMANOVA, 9999 permutations) based on Euclidian dissimilarities with year (2 fixed levels) and sampling site (4 fixed levels) as sources of variation. Variations in FA composition, expressed in percentages, were visualized using non-metric multidimensional scaling (n-MDS).

The similarity percentages (SIMPER) procedure was performed on untransformed data to identify the FAs that explained the most dissimilarity between significant different levels.

RESULTS

Oyster recruitment

Recruitment numbers showed dramatic annual variability with great success at some sites in 2017 but an overall near-zero recruitment level at all sites in 2019 (Fig. 2a,b). In 2017, the metamorphosis survival rate, expressed as the ratio of juvenile to young settler abundances per plate, also showed marked spatial variability. The ratio of juvenile (123 ± 9 ind. plate -1 ) to young-settler abundances per plate (49 ± 6 ind. plate -1 ) was 2.5 in Bouzigues, suggesting up to 100% successful metamorphosis by competent larvae and the arrival of competent larvae from elsewhere. However, in the other sites, recruitment level decreased by 24% (94 ± 16 juveniles plate -1 ) in Meze, 90% (13 ± 2 juveniles plate -1 ) in Listel, and 97% (4 ± 2 juveniles plate -1 ) in Marseillan. A poorer supply of larvae (6 ± 2 young-settlers plate -1 ) was observed in Marseillan, but the metamorphosis survival rate was 0.6. However, in Meze and Listel, the low recruitment rates were not linked to the supply of larvae, as young settler abundances were higher in Meze (328 ± 71 young settler plate -1 , with a metamorphosis survival rate of 0.3) and in Listel (670 ± 65 young settler plate -1 , with a metamorphosis survival rate of 0.02) than in Bouzigues. Failure characterized the 2019 oyster recruitment season: low abundances of young settlers were observed in Meze (116 ± 5 ind. plate -1 ) and in Listel (31 ± 2 ind. plate -1 ), with almost 3 and 22 times fewer individuals than in 2017, respectively. This trend was not observed in Bouzigues (84 ± 9 ind. plate -1 ) or in Marseillan (45 ± 3 ind. plate -1 ) in 2019. Instead, young settlers were respectively 2 and 7 times higher in 2019 than in 2017. However, two weeks later, almost no juveniles were observed on the plates (average 0.14 ± 0.06), regardless of the sites, pointing to a general oyster recruitment failure in 2019.

The size of juveniles at metamorphosis (PII length) was established in all samples, except samples from Bouzigues in 2019 (Fig 2c,d), in which no metamorphosis of young settlers to juveniles was observed. PII individuals sampled in 2019 were 5.1% smaller (mean 262 ± 1 µm) than those sampled in 2017 (mean 276 ± 1 µm). Differences among sites were only observed in 2017, when PII sizes in Bouzigues were 2.7% smaller than those in Meze (p = 0.02), Listel (p = 0.01) and Marseillan (p = 0.03).

No differences in total fatty acid (TFA) contents were observed in young settlers in the four sites and the two years. The overall TFA average was 51 ± 19 ng lavae -1 (p > 0.06). The sum of essential fatty acids (EFA) corresponded to about 10% of TFA with an effect of year × site (df=3 and 19, pseudo-F=6.47, p=0.007), as individuals in Listel (p=0.02) and Marseillan (p=0.006) had 5 times lower TFA contents in 2017 than in 2019. The fatty acid composition of young settlers varied with the year × site interaction (df=3 and 19, pseudo-F=2.34, p=0.017), as individuals sampled in Listel (p=0.047) and Marseillan (p=0.044) had different profiles in the two years (Supplementary material Fig. 1). According to SIMPER analysis, the interannual differences observed at these two sites were linked to DHA (22:6n3), EPA (20:5n3), AA (20:4n6), 18:2n6, 18:0 and 16:0 explained more than 83% of the average dissimilarity in the fatty acid profiles. DHA, EPA and AA levels in young settlers sampled in 2019 were twice higher than in 2017, while the levels of 18:2n6 were five times lower in 2019 than in 2017, except for the Meze and Bouzigues sites (p > 0.09).

Physico-chemical parameters

Average water temperatures were 2.6 o C higher and salinity was 0.34 S higher in 2019 than in 2017 (Fig 3a,b, Supplementary Table 3 andSupplementary Table 4). A site effect was also observed for salinity in the Thau lagoon. Salinity increased from east to west: the mean value at Marseillan was 0.68 S higher than at Bouzigues. Conversely, no effect of site on temperature was observed. There was a site × year effect on oxygen concentration (Supplementary Table 5).

No difference was observed among sites in 2017 (c). The lowest values were observed in Bouzigues in 2019 (p = 0.001) near the bottom of the lagoon (21.8% lower than in 2017).

Oxygen concentrations varied with water depth, lower values generally being observed near the bottom (Fig. 3c).

Potential food for oyster larvae

Concentrations of TPM0.7-20 (), PIM0.7-20 () and POM0.7-20 () were more than twice higher in 2019 than in 2017 (Fig. a,b,c, Supplementary Table 6, 7 and8). Significant differences among the four sites were only observed in POM0.7-20 concentrations. In both years, POM0.7-20 concentrations in Marseillan were 0.7 and 0.8 times lower than in Listel and Meze (p = 0.01 and 0.03 respectively). An effect of year × chl-a biomasses fraction was observed (Supplementary Table 9). Mean nanophytoplankton and picophytoplankton biomasses (p = 0.0001 and p = 0.0004 respectively) were 3 times higher in 2019 (Fig. 4d,e) than in 2017. A site × year effect was also observed, chl-a biomass values were 45% lower in Bouzigues than in Listel (p=0.01) and Meze (p = 0.004) in 2017. In 2019, biomasses in Marseillan were 62% lower than at the other sites (p < 0.02). Interannual variability in chl-a was only found in Bouzigues with 3 times more biomass in 2019 (p = 0.0007) than in 2017. Similar patterns were observed for chlb and chl -c biomass, with twice as much chl-b in the samples collected 2019 samples than in the samples collected in 2017 (0.069 µg L -1 versus 0.026 µg L -1 ; p=0.0001), and a more than two-fold increase in chl-c (0.103 ug L -1 versus 0.046 ug L -1 ), particularly in Listel (p=0.039) and Bouzigues (p=0.0003).

Flow cytometry data showed an effect of the year factor on cells smaller than 3 µm (Fig. ).

Abundances of picoeukaryotes (<3 μm) (Supplementary Table 10), picocyanobacteria (<1 μm) (Supplementary Table 11 and12) and bacteria (Supplementary Table 14) were higher in 2019 than in 2017. However, nanophytoplankton (3-20 μm) abundances decreased by 39% in 2019 (Supplementary Table 13). The abundance of total heterotrophic flagellates did not vary significantly among sites or between years, mean value 2 866 ± 291 cell mL -1 . Dinoflagellate and diatom abundances were affected by the year factor (df=1 and 35, pseudo-F=5.64, p=0.023), total values decreased by 60% in 2019 compared to 2017. These variations were linked to a 93% decrease in Chaetoceros abundance from 184 715 ± 66 846 to 12 483 ± 3 540 cells L -1 (Simper contribution: 77%, df=1 and 35, pseudo-F=8.73, p=0.0001) and a decrease that led to the disappearance of Skeletonema in Listel and Meze between 2017 and 2019. Diatom taxa were fewer in number at all sites sampled in 2019 with a maximum of 13 compared to 21 taxa identified in 2017. A marked increase in Pseudo-nitzschia (19 920 ± 10 513 to 50 562 ± 13 652 cells L -1 ) with a Simper contribution of 8% and (df=1 and 35, pseudo-F=8.73, p=0.0001), Leptocylindrus (Simper contribution 7%), Thalassionema, and Cylindrotheca (1 837 ± 222 to 18 712 ± 12 010 cells L -1 ) was observed in 2019 compared to 2017. This trend is especially expressed in Bouzigues (Fig. 6). This result also reflects the higher diversity of dinoflagellate taxa observed in 2019 (16 taxa) than in 2017 (12 taxa).

TFA contents in the TPM0.7-20 samples were twice higher in 2019 (19.2 µg mg TPM0.7-20 -1 ) than in 2017 (9.9 µg mg TPM0.7-20 -1 ; df=1 and 61, pseudo-F=17.1, p=0.0002) with no differences among sites and year × site effects. The fatty acid composition of the TPM0.7-20 samples differed between years (df=3 and 76, pseudo-F=3.08, p=0.0001; Fig. S2) and, as determined by the SIMPER analysis, explained 97% of the differences in the levels of 18:1n9, 18:0, 16:1, 18:2n6, 16:0, 14:0, 20:5n3 and 22:6n3. Twenty-six percent of the difference observed between years was related to 18:1n9, a FA that was twice as abundant in 2017 (up to 24.1% of the TFA) compared to 2019. The dissimilarity in the FA profiles observed between years was also explained by higher values of 18:2n6 (representing up to 10.8% of TFA), and EPA (7%) in 2017. 18:2n6 and EPA were, respectively, 11.3% and 5% higher in 2017 than in 2019. The most abundant FAs in the TPM0.7-20 samples in 2019 were 16:1 and DHA which explained, respectively, 13% and 4.3% of the dissimilarity shown by the SIMPER analysis.

Territorial competition by worms

The percent cover of tubeworms (Ficopomatus enigmaticus) on plates in 2019 showed a marked increase in this species. Differences were observed among the sites (df=3 and 33, pseudo-F=157, p=0.0001). Results showed similar cover of tubeworms (93.6 ± 1.5%) in Listel and Bouzigues and a lower cover in Meze (83.2 ± 2.6%) (p < 0.032) and in Marseillan (23.6 ± 3.7%) (p < 0.0001).

DISCUSSION

The aim of this study was to identify the environmental and trophic drivers of the decline in the recruitment of the Pacific oyster, Crassostrea gigas, in association with a heat wave. Our hypothesis that a heat wave has a negative effect on oyster recruitment by altering plankton diversity was confirmed. While oyster recruitment was normal in 2017, an unprecedented failure was observed in summer 2019 in the Mediterranean Thau lagoon. The atmospheric conditions during a heat wave have a strong direct effects on marine and lagoon environments that supply a variety of ecosystem services and valuable host species [START_REF] Sarà | Multiple climatedriven cascading ecosystem effects after the loss of a foundation species[END_REF]. Temperature and salinity conditions are key ecological and physiological factors for Crassostrea larvae (His et al. 1989b, Baldwin & Newell 1995a[START_REF] Devakie | Salinity-temperature and nutritional effects on the setting rate of larvae of the tropical oyster, Crassostrea iredalei (Faustino)[END_REF][START_REF] Troost | Effects of an increasing filter feeder stock on larval abundance in the Oosterschelde estuary (SW Netherlands)[END_REF]. In controlled experimental settings, the entire larval life of C. gigas, including metamorphosis, showed a high tolerance to temperatures ranging from 17 °C to 32 °C at a salinity level of 34, with low mortality (≤10%) and the maximum growth rate at 32 °C [START_REF] Rico-Villa | Influence of food density and temperature on ingestion, growth and settlement of Pacific oyster larvae, Crassostrea gigas[END_REF]). The physiological limits of temperature tolerance were therefore not reached in our experimental conditions and temperature was not the origin of the failure in this case. Salinity did not drop below 38 in either the 2017 or 2019 recruitment season, and intermittently reached more than 40 in 2019. Crassostrea gigas is an estuarine organism that tolerates a wide range of salinity [START_REF] Nell | Effects of Salinity on the Growth and Survival of Sydney Rock Oyster (Saccostrea commercialis) and Pacific Oyster (Crassostrea gigas) Larvae and Spat[END_REF], but no information is available in the literature on the upper salinity tolerance of the larval stage in real conditions. The high salinity in 2019 could represent the physiological salinity threshold for oyster larvae. Our results showed that the larval shell (prodissoconch) at the time of metamorphosis (PII) was smaller in 2019, suggesting a reduction in larval growth or faster achievement of metamorphosis competence in high salinity years. In agreement with [START_REF] Nell | Effects of Salinity on the Growth and Survival of Sydney Rock Oyster (Saccostrea commercialis) and Pacific Oyster (Crassostrea gigas) Larvae and Spat[END_REF] who reported an optimal salinity range for larval growth up to 27 and a very marked growth reductions at 31-39 [START_REF] Nell | Effects of Salinity on the Growth and Survival of Sydney Rock Oyster (Saccostrea commercialis) and Pacific Oyster (Crassostrea gigas) Larvae and Spat[END_REF], the smaller observed PII size could be related to growth limitation under high salinity. Interestingly, these authors reported no significant effect of salinity on larval survival between 19 and 39 but the growth rate of larvae decreased markedly from 30 S [START_REF] Helm | Experiments in the hatchery rearing larvae[END_REF]. Upper tolerance limits of oysters to high salinity ranging from 35 S to 45 S should thus be further tested in the laboratory including interactions with high temperature and different nutritional inputs (His et al. 1989a).

The heat wave that occurred in 2019 resulted in large quantities of particulate matter and chlorophyll biomass, but their quality appeared to be unfavorable for oyster recruitment. The failure of oyster recruitment in 2019 could thus be linked to the change in the phytoplankton community with low abundance of forage diatoms and high abundance of picoplanktonic prokaryotes and eukaryotes, flagellates, and of the diatoms Pseudo-Nitzschia and Cylindrotheca. However, the trophic environment was not characterized by a planktonic community poor in fatty acids, and it was in fact richer than in 2017. Pediveliger larvae accumulated the same quantity of fatty acids in 2019 as in 2017, but metamorphosis failures were observed at all sites. We suggest that this failure may be linked to inappropriate trophic conditions, which in turn, are mainly linked to the size of picoplankton species. These species are poorly retained by the newly developed gills of young metamorphosed juveniles. Except for larvae, the retention efficiency of bivalves for particles < 3-4 µm is low (Baldwin & Newell 1995b[START_REF] Rosa | Selective Capture and Ingestion of Particles by Suspension-Feeding Bivalve Molluscs: A Review[END_REF]. Our results suggest that the overabundance of small particles (picoplanktonic prokaryotes and eukaryotes) could be critical for larval settlement and metamorphosis. Higher chlorophyll biomass was observed in the nanophytoplankton fraction during the heat wave of 2019 than in 2017, indicating changes in the phytoplankton community.

The heat wave was characterized by the increasing abundances of picocyanobacteria [START_REF] Bec | Phytoplankton seasonal dynamics in a Mediterranean coastal lagoon: Emphasis on the picoeukaryote community[END_REF][START_REF] Collos | Oligotrophication and emergence of picocyanobacteria and a toxic dino fl agellate in Thau lagoon , southern France[END_REF], Derolez et al. 2020b) and decreasing abundances of nanophytoplankton. The Thau lagoon began an oligotrophication trajectory in the early 2000s [START_REF] Collos | Oligotrophication and emergence of picocyanobacteria and a toxic dino fl agellate in Thau lagoon , southern France[END_REF], Derolez et al. 2020a). This process caused a community shift due to a reduction in nutrient loads since the 1970s thanks to improved wastewater treatment in the watershed aimed at halting eutrophication (EC 1991a(EC b, 2000)). The reduction in nutrient loads has been amplified by a decrease in total rainfall since the 2000s due to climate change (Derolez et al. 2020a). Our results corroborate evidence that the proportion of small taxa, like picoplankton, in the phytoplankton community, is increasing in coastal, marine and freshwater ecosystems in response to global warming [START_REF] Daufresne | Global warming benefits the small[END_REF][START_REF] Mousing | Global patterns in phytoplankton community size structureevidence for a direct temperature effect[END_REF][START_REF] Pinckney | Phytoplankton community structure and depth distribution changes in the Cariaco Basin between 1996 and 2010[END_REF]. Small phytoplankton cells have been observed to dominate in oligotrophic environments [START_REF] Irwin | Scaling-up from nutrient physiology to the size-structure of phytoplankton communities[END_REF]).

The high temperatures and high salinity in 2019 had a negative impact on trophic conditions for larvae. However, recruitment failure appeared to be more linked to the ecological limitations of the larvae at the time of their metamorphosis than to their physiological state. At the same time, high temperatures and high salinity stimulated the development of the annelid Ficopomatus enigmaticus, triggering a shift in community composition that is destructive for oyster recruitment. We consequently hypothesize that these annelids are important territorial competitors [START_REF] Heiman | Non-native Ecosystem Engineer Alters Estuarine Communities[END_REF][START_REF] Mcquaid | Alien reef-building polychaete drives long-term changes in invertebrate biomass and diversity in a small, urban estuary[END_REF][START_REF] Peria | Tolerance to salinity and thermal stress by larvae and adults of the serpulid annelid Ficopomatus enigmaticus[END_REF] and trophic competitors of oyster larvae [START_REF] Davies | The filtration activity of a serpulid polychaete population Ficopomatus enigmaticus (Fauvel) and its effects on water quality in a coastal marina[END_REF][START_REF] Bruschetti | Grazing effect of the invasive reef-forming polychaete Ficopomatus enigmaticus (Fauvel) on phytoplankton biomass in a SW Atlantic coastal lagoon[END_REF], 2018[START_REF] Pan | Top-Down Effects of an Exotic Serpulid Polychaete on Natural Plankton Assemblage of Estuarine and Brackish Systems in the SW Atlantic[END_REF] in shallow water and brackish habitats.

This study demonstrates, for the first time, an ecological process leading to the recruitment failure of the Pacific oysters due to an extreme heat wave. The oligotrophication trajectory of our study site combined with the effects of high water temperatures promoted variations in the phytoplankton communities that benefit picophytoplankton including cyanobacteria, that are likely unfavorable for the successful larval development of oysters until their juvenile metamorphosis [START_REF] Lagarde | Recruitment of the Pacific oyster Crassostrea gigas in a shellfish-exploited Mediterranean lagoon : discovery, driving factors and a favorable environmental window[END_REF]. The present study thus reveals the ecological limits of the reproductive process of the Pacific oyster in the context of a heat wave in a Mediterranean lagoon. The heat wave phenomenon observed in 2019 severely disrupted the reproductive cycle of oysters in the Thau lagoon. In this context, the oyster nursery function within an oyster farming ecosystem can only be achieved or maintained when pico-, nano-and microphytoplankton communities are present and abundant and oysters can find favorable areas for larval development and optimize their recruitment. This study provides evidence that, in the conditions caused by a heat wave, the ecological limits of Pacific oyster larvae are narrower than their physiological limits. The effects of climate change, particularly the warming of waters in semi-enclosed basins, will certainly lead to problems in larval harvesting in the near futures.

The information presented in this paper should help adapt oyster aquaculture, including husbandry practices, to a future marked by climate change.

Data and code availability

All data used in the current study and scripts used in our analysis are publicly available or were obtained by the corresponding author. This research benefited from the VELYGER Database:

The Oyster Larvae Monitoring French Project (http://doi.org/10.17882/41888) and REPHY and M.H. thank FREA, JSPS and Campus France/Ministry of Foreign Affairs for funding the scientific exchange needed for this study and the development of all the ideas presented by our group.

  in 10 mL seawater samples. Abundances are expressed as the number of individuals per liter in 52 diatom taxa and 38 dinoflagellate taxa. Chlorophyll a (Chl-a), b (Chl-b) and c (Chl-c) biomasses

Fig. 1 .

 1 Fig. 1. The four sampling sites in the Thau lagoon. Marseillan and Bouzigues are located in the

Fig. 2 .

 2 Fig. 2. Variability of oyster recruitment and prodissoconch II size according to the years 2017

Fig. 3 .Fig. 4 .Fig. 5 :

 345 Fig. 3. Physico-chemical monitoring in 2017 (no heat wave) and 2019 (heat wave). (a) Mean temperature (°C ± SE) per year (n = 40), (b) mean salinity (PSU ± SE) per year (n = 40) and per sampling site (n = 20) and (c) mean dissolved oxygen concentration (mg L -1 ± SE) according to the position of the sample in the water column (n = 40) and per year and sampling site (n = 10). Stars indicate significant differences in parameter average per year (* p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001). Different letters indicate significant differences between sites according to post hoc multiple comparison tests after PERMANOVA.

Fig. 6 :

 6 Fig. 6: Heatmap of microphytoplankton genera with changes in abundances in 2017 (no heat wave) and 2019 (heat wave). Average phytoplankton abundance (cells L -1 ) per taxon and sampling site in (a) 2017 (n = 5) and (b) in 2019 (n = 4).

We thank Hélène Cochet, Serge Mortreux, Anaïs Crottier, Clarisse Hubert, Nabila Guenineche, Gabriel Devique, Hervé Violette, Elise Hatey, Camille Gianaroli and Nathalie Gauthier for logistic support in the field and technical contributions in the laboratory.

ACKNOWLEDGEMENTS

This research was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC-Discovery Grant no.299100) to R.T and by the Ressources Aquatiques Québec Research Network (Fonds de Recherche du Québec-Nature et Technologies, #2014-RS-171172) and MITACS (FR37656) for the internships of ACM during her stay at IFREMER. This research was also fund by the VELYGER network. F.L. and R.T. thank the RECHAGLO international research group, co-funded by Ifremer and the Department of Fisheries and Oceans (DFO), and the Institut France-Québec pour la coopération scientifique en appui au secteur Maritime (IFQM) for encouragement, support, and exchanges with Canada. F.L., M.H. T.M.

AUTHOR CONTRIBUTIONS

A.C.M. was involved in investigation, methodology, writing, data curation, formal analysis, and visualization. R.T. and F.L. were involved in conceptualization, funding acquisition, investigation, methodology, writing, data curation, formal analysis, visualization, and project administration. S.P. was involved in conceptualization, funding acquisition, investigation, methodology, writing and project administration. B.B was involved in conceptualization, funding acquisition, investigation, methodology, writing, data curation, formal analysis, and visualization. C.R. contributed to funding acquisition, methodology, writing, data curation and formal analysis. A.A and A.G. contributed to writing and interpretation. G.M. contributed to funding acquisition, investigation, methodology, writing and formal analysis. M.R., M.Ho, M.Ha. and T.M. contributed to conceptualization, investigation, methodology and writing.

COMPETING INTERESTS

The authors declare no competing interests. 3 0,3865 0,12883 0,063154 0,9764 9951 0,9754 sitexposition 3 2,573 0,85767 0,42042 0,7357 9950 0,7371 yearxposition 1 1,1045 1,1045 0,54142 0,4681 9828 0,473 sitexyearxposition 3 0,0865 0,028833 0,014134 0,9977 9955 0