Frida Ben

Rais Lasram
email: frida.lasram@univ-littoral.fr

Tarek Hattab

Quentin Nogues

Grégory Beaugrand

Jean Claude Dauvin

Ghassen Halouani

François Le Loc'h

Nathalie Niquil

Boris Leroy

An open-source framework to model present and future marine species distributions at local scale

Highlights

-A new habitat modelling open-source framework for non-modeller ecologists.

-An explicit consideration of the third dimension and habitat filter in the modeling procedure.

-An automated procedure for data processing.

-The proposed framework avoids higher projections of species loss.

We refer here to all models that correlate species occurrences or abundance data with environmental spatial data layers to predict the suitability of any site for a given species. Shifts in species distributions are subsequently projected under different climate change scenarios. These SDMs are correlative and do not incorporate mechanistic information, such as species' physiological response to abiotic conditions. In the marine realm, there is a flourishing literature carried at global (e.g. [START_REF] Tittensor | Seamounts as refugia from ocean acidification for cold-water stony corals[END_REF], regional (e.g. Albouy et al. 2013) and local (e.g. [START_REF] Hattab | Towards a better understanding of potential impacts of climate change on marine species distribution: a multiscale modelling approach[END_REF]) scales that have highlighted future species range shift towards the poles (e.g. [START_REF] Morley | Projecting shifts in thermal habitat for 686 species on the North American continental shelf[END_REF], deepening (e.g. [START_REF] Pinsky | Marine Taxa Track Local Climate Velocities[END_REF], local extinctions (e.g. [START_REF] Jones | Multi-model ensemble projections of climate change effects on global marine biodiversity[END_REF], habitat fragmentation (e.g. Ben Rais [START_REF] Rais Lasram | The Mediterranean Sea as a'cul-de-sac'for endemic fishes facing climate change[END_REF], species invasions (e.g. [START_REF] Byrne | From pole to pole: the potential for the Arctic seastar Asterias amurensis to invade a warming Southern Ocean[END_REF]) and impacts on life-history traits (e.g. [START_REF] Genner | Body size-dependent responses of a marine fish assemblage to climate change and fishing over a century-long scale[END_REF], on abundances (e.g. [START_REF] Hermant | Impact of warming on abundance and occurrence of flatfish populations in the Bay of Biscay (France)[END_REF] and on trophic networks (e.g. [START_REF] Doney | Climate Change Impacts on Marine Ecosystems[END_REF]) due to changing climate in the next few decades. SDMs have been applied at either the species (e.g. Alabia et al. 2015) or community level [START_REF] Morley | Projecting shifts in thermal habitat for 686 species on the North American continental shelf[END_REF]) to assess species vulnerabilities for conservation purposes [START_REF] Jones | Predicting the impact of climate change on threatened species in UK waters[END_REF]) and fisheries management [START_REF] Asch | Future marine ecosystem drivers, biodiversity, and fisheries maximum catch potential in Pacific Island countries and territories under climate change[END_REF], and to couple with end-to-end models to simulate spatial patterns of predator-prey interactions [START_REF] Grüss | Producing distribution maps for a spatially-explicit ecosystem model using large monitoring and environmental databases and a combination of interpolation and extrapolation[END_REF]. Moreover, trophic ecologists are increasingly interested in using SDMs predictions as input for trophic models to investigate how food webs may reorganize due to individual species' responses to climate change [START_REF] Woodward | Climate change and freshwater ecosystems: impacts across multiple levels of organization[END_REF][START_REF] Chaalali | From species distributions to ecosystem structure and function: a methodological perspective[END_REF]. Although SDMs are useful and reliable tools to predict past, present and future species ranges, they rely on strong assumptions and are subject to many uncertainties [START_REF] Goberville | Uncertainties in the projection of species distributions related to general circulation models[END_REF]. These uncertainties are well known and occur at every stage of the modelling process (e.g. data collection, choice of algorithms and parameters, model evaluation), but few studies consider them explicitly [START_REF] Planque | Projecting the future state of marine ecosystems,"la grande illusion[END_REF][START_REF] Brun | Model complexity affects species distribution projections under climate change[END_REF]. Moreover, despite the ability to run SDMs in recent, reliable and user-friendly packages and the availability of data in huge global databases (e.g. GBIF, OBIS and FishBase for species occurrences; Bio-ORACLE and MarSpec for environmental variables), there are some misuses and proxies (e.g. data quality, selection of environmental drivers, parametrization, model selection, spatial and temporal scales) that are widely practiced and rarely addressed in the literature (see [START_REF] Jarnevich | Caveats for correlative species distribution modeling[END_REF][START_REF] Yalcin | Diversity and suitability of existing methods and metrics for quantifying species range shifts[END_REF][START_REF] Duffy | Explicitly integrating a third dimension in marine species distribution modelling[END_REF] for examples of reviews and recommendations). Although global georeferenced data have become more available for both species occurrence/abundance data and environmental variables, less concern has been given to assessing data quality in the marine realm. For example, sampling bias due to diverse and/or non-standardized monitoring surveys may increase the risk of undersampling (e.g. due to selectivity of a sampling method) or oversampling (e.g. in neritic regions or in areas with a long history of monitoring). Consequently, occurrence data often represent a biased sample of species populations, which may significantly alter models calibration [START_REF] Guillera-Arroita | Is my species distribution model fit for purpose? Matching data and models to applications[END_REF]. Careful assessment and preprocessing of data are thus recommended to mitigate effects of sampling bias associated with largescale occurrence data. Most studies in the marine realm have been based on two-dimensional SDMs even though species are fundamentally distributed in three dimensions. [START_REF] Duffy | Explicitly integrating a third dimension in marine species distribution modelling[END_REF] emphasized the importance of incorporating the third dimension in marine ecological models and showed that ignoring it may influence predictions greatly. Some studies (e.g. [START_REF] Weinmann | Heading for new shores: projecting marine distribution ranges of selected larger foraminifera[END_REF][START_REF] Gallardo | The importance of the human footprint in shaping the global distribution of terrestrial, freshwater and marine invaders[END_REF][START_REF] Asch | Future marine ecosystem drivers, biodiversity, and fisheries maximum catch potential in Pacific Island countries and territories under climate change[END_REF]) continue to use sea-surface temperature for both pelagic and benthic species despite the temperature gradient between the surface and the bottom. Moreover, species-environment relationships depend strongly on the scale at which the dependent (i.e. species occurrences) and independent (i.e. environmental data) variables are considered [START_REF] Cushman | Patterns in the species-environment relationship depend on both scale and choice of response variables[END_REF]. Scale issues are thus important to consider in SDMs especially for the downscaling process [START_REF] Da Re | Downscaling livestock census data using multivariate predictive models: Sensitivity to modifiable areal unit problem[END_REF]). It is for example important to choose the appropriate grid size and the relevant variables consistent with the scale at which the ecophysiological processes show greatest variation [START_REF] Austin | Improving species distribution models for climate change studies: variable selection and scale[END_REF]. . The common view is that large-scale processes determine local species diversity, whereas small-scale processes tend to limit the number of species that currently occupy a given site [START_REF] Pont | Modelling habitat requirement of European fishes: do species have similar responses to local and regional environmental constraints?[END_REF]. To address this issue, a two-stage hierarchical approach based on the concept of "hierarchical filters" is commonly (but not always) used in the terrestrial realm but more rarely in the marine realm (see [START_REF] Heino | Climate change and freshwater biodiversity: detected patterns, future trends and adaptations in northern regions[END_REF][START_REF] Hattab | Towards a better understanding of potential impacts of climate change on marine species distribution: a multiscale modelling approach[END_REF] for a review of underlying assumptions). In this filtering approach, SDMs are first built at the global scale to consider the full range of climatic conditions (temperature and salinity) species encounters (SDMs can then be termed as Bioclimatic Envelope Models (BEMs)). Then, a habitat filter considering factors at the local scale (such as seafloor type) is applied (the filter can then be termed as Habitat Models (HMs)) . This twostage approach improves predictions of species distributions at the local scale. Finally, as policy and decision making require spatially explicit information about the state of and future trends in species ranges, it is useful to generate current and future maps, along with response curves and model performances, for each species. A broad set of tools exists to model species distribution (e.g. BIOMOD [START_REF] Thuiller | BIOMOD a platform for ensemble forecasting of species distributions[END_REF], MAXENT [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF]), openModeller (de Souza Muñoz et al. 2011), ModEco [START_REF] Guo | ModEco: an integrated software package for ecological niche modeling[END_REF], ENMTool [START_REF] Warren | ENMTools: a toolbox for comparative studies of environmental niche models[END_REF], sdm [START_REF] Naimi | Sdm: a reproducible and extensible R platform for species distribution modelling[END_REF], NPPEN [START_REF] Beaugrand | A new model to assess the probability of occurrence of a species based on presence-only data[END_REF]) and dismo package [START_REF] Hijmans | Package 'dismo[END_REF]). These platforms include a comprehensive set of advanced algorithms that generally need to be supplemented by R scripts to pre-process data or to be coupled with GIS softwares, making them somewhat difficult for non-modeller ecologists to use. In this paper, we propose a framework that includes (i) a procedure for homogenizing occurrence data to eliminate the influence of sampling bias, (ii) a procedure for generating pseudo-absences, (iii) a hierarchical-filter approach (i.e. global Bioclimatic Envelope Models (BEM) combined with local Habitat Models HM), (iv) full incorporation of the third dimension by considering climatic variables at multiple depths and (v) building of maps that predict current and future ranges of marine species. This framework is available as a ready-to-use R script hosted on GitHub for non-modeller ecologists interested in investigating future species ranges with a user-friendly script. We investigated the robustness of the framework by applying it to 46 species representing a large part of the biomass in the Eastern English Channel (EEC). Projections were built for the middle and the end of this century under Representative Concentration Pathway (RCP) 2.6 (i.e. rapid mitigation of anthropogenic climate change, optimistic) and RCP8.5 (i.e. business-as-usual, pessimistic) scenarios of the IPCC (Intergovernmental Panel on Climate Change).

Model framework

Although designed for the ecosystem of the EEC, our script can be applied to other ecosystems by replacing the habitat data files provided with the script with similar data files for another region. To illustrate application of the script, we considered 46 species representing a large part of the biomass of the EEC [START_REF] Raoux | Benthic and fish aggregation inside an offshore wind farm: Which effects on the trophic web functioning?[END_REF]), with at least 300 occurrences at the global scale and 100 occurrences at a regional scale. We provide the script and environmental data online, at https://github.com/TarekHattab/SDM. For the sake of reproducibility, the script comes with a Docker container [START_REF] Merkel | Docker: lightweight linux containers for consistent development and deployment[END_REF] including everything needed to run the script. The model framework is described in Figure 2.

Step 1: Setting global parameters

The script is designed to manage automatically all of the framework's processing steps, described below. Users must set four parameters: (i) species' scientific names, (ii) species' vertical habitats, (iii) the algorithms to be used and (iv) the value of K in K-fold cross-validation. Our script cleans the global dataset by deleting erroneous occurrences (e.g. species occurring on land) and occurrences outside a user-defined period. In this study, we deleted records before 1955 and ultimately retained 1 944 154 occurrences from the five databases for the 46 species considered. Additional scientific survey data are available at the regional scale for the EEC and the Bay of Biscay. Thus, the global dataset was supplemented by scientific survey data from IBTS (1990-2016, http://dx.doi.org/10.18142/17), CGFS (1988-2016, http://dx.doi.org/10.18142/11) and EVHOE (1997-2006, http://dx.doi.org/10.18142/8) available from the DATRAS database (http://www.ices.dk/marinedata/data-portals/Pages/DATRAS.aspx) and also from other surveys of benthic and demersal taxa, such as LANCIE (1997)(1998)(1999)(2000)(2001)(2002), PECTOW(2001), BENTHOSEINE (1998BENTHOSEINE ( -1999) ) and CAMANOC (http://dx.doi.org/10.17600/14001900) (Appendix A). Ultimately, we retained 148 096 occurrences at the regional scale.

Step 2: Data acquisition 2.2.1 Species occurrence data

Current and projected climatic data

To calibrate BEMs, users can consider temperature and salinity climatologies from the global database WOD 2013 V2 (https://www.nodc.noaa.gov/OC5/woa13/), with a spatial resolution of 0.25°. These climatologies represent mean decadal temperatures and salinities for 1955-1964, 1965-1974, 1975-1984, 1985-1994, 1995-2004 and 2005-2012 for 40 depth layers. These variables were bilinearly interpolated at a 1/12 th -degree spatial resolution (5 arcmin) and aggregated vertically by calculating mean temperature and salinity for the first 50 m of depth to calibrate pelagic species models, for the first 200 m of depth for benthopelagic species models and for the last 50 m of depth for benthic and demersal species models. For projections, the script performs climate projections for 2041-2050 and 2091-2100 under two RCP2.6 (strong mitigation) and RCP8.5 (business-as-usual) scenarios. These projections are derived from three General Circulation Models (GCMs) of the Coupled Model Intercomparaison Project, Phase 5 (CMIP5) [START_REF] Taylor | An overview of CMIP5 and the experiment design[END_REF]): (i) ESM2G (Geophysical Fluid Dynamics Laboratory; NOAA), (ii) CM5A-MR (Institut Pierre-Simon-Laplace) and (iii) ESM-MR (Max-Planck-Institut für Meteorologie). To avoid potential bias of the GCMs compared to historical data, the script calculates anomalies between 2041-2050/2091-2100 and 1955-2012. These anomalies are then added to the observed mean temperatures and salinities for the same period in WOD 2013 data. Current and projected climatic data are provided with the script.

Habitat data

Seafloor type, bathymetry, slope and aspect of the slope (i.e. eastness and northness) are used to build HMs for benthic and demersal species. Data are available from EMODnet Bathymetry (http://www.emodnet-bathymetry.eu/) and EMODnet Seabed Habitats (http://www.emodnetseabedhabitats.eu/) at a spatial resolution of 250 m. In our study, we aggregated EUNIS seafloor substrate classes into six simplified categories (mud, sand, muddy sand, sandy mud, coarse sediment and gravel). Habitat data grids are provided with the script.

2.3

Step 3: Data pre-processing Once occurrence data have been downloaded, the script performs a spatiotemporal match between climatologies and occurrence of each species based on geographic coordinates of occurrences and their corresponding decade (to consider climate trends over the past 58 years) and the correct vertical layer, which corresponds to the vertical habitat of the species. To address the risk of oversampling, and the ensuing over-representation of environmental features [START_REF] Kramer-Schadt | The importance of correcting for sampling bias in MaxEnt species distribution models[END_REF], we first applied an environmental filtering procedure [START_REF] Beaugrand | A new model to assess the probability of occurrence of a species based on presence-only data[END_REF][START_REF] Varela | Environmental filters reduce the effects of sampling bias and improve predictions of ecological niche models[END_REF] to species occurrences to assign the same weight to over-and under-sampled regions.

To do this, the script creates a grid containing all combinations of temperature and salinity at the global scale for the period 1955-2012 with a resolution of 0.3°C × 0.3 salinity (which represents 100 × 100 bins). Species occurrences are then projected onto this grid, and only one occurrence is retained in each cell of the grid The datasets we used were obtained from online databases of occurrences including presence-only data, but the BEMs used in the script require both presence and absence data. The script therefore generates pseudo-absences to better characterize the environmental conditions experienced by a species within the environmental background [START_REF] Hattab | Towards a better understanding of potential impacts of climate change on marine species distribution: a multiscale modelling approach[END_REF]. Pseudo-absences were generated outside the convex hull of presences in the environmental space [START_REF] Cornwell | A trait-based test for habitat filtering: convex hull volume[END_REF], Getz & Wilmers 2004). The convex hull was defined as the smallest convex area in the environmental space containing all species records. A restricted convex hull was defined as a convex hull excluding occurrence points in the 1 st and 99 th percentiles for each environmental parameter in order to exclude observations in the most extreme environmental conditions. Thus, pseudo-absences are randomly selected by the script outside the environmental space described by this restricted convex hull, which is a proxy of the suitable environmental conditions (Figure 1). In HMs (SDMs used to characterize habitat), pseudo-absences were generated similarly but were based on the "Surface Range Envelope" [START_REF] Thuiller | BIOMOD a platform for ensemble forecasting of species distributions[END_REF]. In this approach, rectilinear hypervolumes describing the 1 st and 99 th quantiles of the cells containing presences are used instead of the convex hull. Based on the statistical theory of model-based designs (d-designs; [START_REF] Hengl | Spatial prediction of species' distributions from occurrence-only records: combining point pattern analysis, ENFA and regression-kriging[END_REF]), the script generates as many pseudo-absences as presences. Figure 1. Generation of pseudo-absences. Grey cells represent the environmental space as a combination of temperature and salinity at the global scale. Green cells represent the presence of a given species. The grid represents the restricted convex hull describing the 1 st and 99 th quantiles of the cells containing presences. Red cells represent pseudo-absences, which are generated randomly outside the convex hull.

Step 4: Bioclimatic envelope models BEMs

Filtered presence data and the generated pseudo-absences were used to build BEMs that consisted of SDMs applied at the global scale to characterize mainly the bioclimatic envelope of each species according to its temperature and salinity preference. For each species, the script applies eight modelling techniques belonging to four model categories using the BIOMOD multi-model platform [START_REF] Thuiller | BIOMOD a platform for ensemble forecasting of species distributions[END_REF]): multiple regressions (Generalized Linear Model, Generalized Additive Model, Multiple Adaptive Regression Splines), regression trees (Boosted Regression Tree, Random Forest, Classification Tree Analysis), discriminant analysis (Flexible Discriminant Analysis) and learning techniques (Artificial Neural Network). For each species, the eight models were calibrated using a random sample of the initial data (67%), while the remaining 33% were used to evaluate the model.

Users can perform a 3-fold cross validation procedure and assess model performance using both the True Skill Statistic (TSS; [START_REF] Allouche | Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS)[END_REF]) and Continuous Boyce Index (CBI; [START_REF] Hirzel | Evaluating the ability of habitat suitability models to predict species presences[END_REF]. Note that the script uses mean temperature and salinity from the first 50 m of depth to calibrate pelagic species models, from the first 200 m of depth for benthopelagic species models and from the last 50 m of depth for benthic and demersal species models. For each species, 24 predictions of current climate envelopes (8 algorithms × 3 permutations) were calculated. Only models with CBI > 0.5 for all three permutations were retained (i.e. models that are resilient to occurrence permutations) for projections of future potential bioclimatic envelopes for each species according to RCP2.6 and RCP8.5 scenarios. Finally, the script generates ensemble suitability maps by calculating the mean suitability of the 24 predictions weighted by the CBI to estimate modelbased uncertainty [START_REF] Thuiller | BIOMOD a platform for ensemble forecasting of species distributions[END_REF]). These maps could then be transformed into binary maps using a probability threshold that optimized the TSS. In addition, the script generates uncertainty maps by calculating the standard deviation of predictions.

Step 5: Habitat models HMs

HMs are SDMs used to characterize local habitat variability. Because habitat variables can be strongly intercorrelated, the script can remove multicollinearity by using the ordination method of [START_REF] Hill | Principal component analysis of taxonomic data with multi-state discrete characters[END_REF]. The resulting orthogonal and uncorrelated principal components are used as predictive variables in the HMs. The HMs, used only for benthic and demersal species, are built using the same method used to develop the BEMs.

Step 6: Combining BEMs and HMs

For benthic and demersal species, users can extract the predicted bioclimatic envelope for their study area from the projected BEM maps. This represents the first hierarchical filter. Next, the habitat map is used to filter out habitats deemed unsuitable by the HM. The extracted cells can be resampled on a grid (we used one with 250 m resolution) using bilinear interpolation. Thus, an area is considered suitable for a species if both the first filter (BEM) and second filter (HM) predict it as suitable.

Step 7: Visualizing predictions

The final stage of the framework automatically plots current and future species ranges for both periods (2041-2050 and 2091-2100) under both scenarios (RCP2.6 and RCP8.5) using only BEM for pelagic and benthopelagic species and combined BEM and HM for benthic and demersal species. Maps are plotted in both suitability index and binary format. Response curves, model performances and number of occurrences can also be provided for each species.

Focus on the Eastern English Channel

We applied this modelling framework to the EEC, an epicontinental sea with relatively shallow water (maximum 50 m), subject to a variety of hydrodynamic forces, characterized by a strong tidal range (up to 9 m) and under the influence of river inputs, mainly from the Seine River on the French side [START_REF] Carpentier | Channel Habitat Atlas for marine Ressource Management, final report/Atlas des Habitats des Ressources Marines de la Manche Orientale[END_REF]) (Figure 3). We considered 46 species that dominate total biomass in the EEC: 15 pelagic and benthopelagic (13 bony fishes and 2 cephalopods), to which we applied BEMs, and 31 benthic and demersal (12 bony fishes, 5 cartilaginous fishes, 6 molluscs, 4 echinoderms, 3 annelids and 1 cnidarian), to which we applied BEMs filtered with HMs (Appendix B). By the end of the century (2091-2100), the GCMs ESM2G, CM5A-MR and ESM-MR predicted increases in temperature in the EEC of 0.54, 0.79 and 0.73°C, respectively, under RCP2.6 and 1.5, 2.4 and 1.7°C, respectively, under RCP8.5 . For the same period, for ESM2G, CM5A-MR and ESM-MR, salinity anomalies were 0.12, -0.67 and -0.06, respectively, under RCP2.6 and -0.41, -1.45 and -1.11, respectively, under RCP8.5. Figure 3. The Eastern English Channel, study area for the application of the script 3.1 Potential changing patterns of species distributions according to the proposed approach By comparing predictive maps of current species richness to those of future projections in response to changing climate, potential impact of climate change on species' assemblages can be assessed. Species richness in the EEC currently ranges from 19-46 species. Under RCP2.6, by 2041-2050, species richness was projected to decrease by a mean of 7.8% (i.e. -3.5 species for the entire EEC) and by a maximum of 13.0% (i.e. -6.0 species). The Bay of Seine and the central EEC would lose the most species. Under RCP2.6, by 2091-2100, the decrease would be a mean of 11.5% (i.e. -5.3 species) and a maximum of 19.6% (i.e. -9.0 species). As expected, more species would be lost under RCP8.5, reaching 17.4% of the initial pool by 2041-2050 (-8.0 species) and 54.4% by 2091-2100 (-25.0 species) (Figure 4). Overall, for both scenarios and both periods, no increase in species range was projected, as none of the species considered was at its northern limit. Species range was projected to be maintained in the EEC for 42 species by 2041-2050 and 39 species by 2091-2100 under RCP2.6 and 41 species by 2041-2050 and 21 species by 2091-2100 under RCP8.5. Under RCP2.6, by 2041-2050, only three benthic species (the fishes Pleuronectes platessa and Limanda limanda and the gastropod Buccinum undatum) and one benthopelagic species (the fish Gadus morhua) were projected to lose 100% of their current range in the EEC (i.e. local extinction). By 2091-2100, two other benthic species (the fish Platichthys flesus and the sea urchin Psammechinus miliaris) and one benthopelagic species (the fish Merlangius merlangus) were projected to lose 50-75% of their ranges. Under RCP8.5, by 2041-2050, two benthic species (the fish P. flesus and the sea star Asterias rubens), in addition to those projected to go locally extinct under RCP2.6, were projected to lose 50% and 75% of their ranges, respectively. Under RCP8.5, by 2091-2100, 19 species, both benthic/demersal and pelagic, were projected to go locally extinct.

Classic approach vs. the proposed framework

According to the complete approach developed in our methodological framework (namely the proposed framework), current species richness was projected to decline by a mean of 3.6% by 2041-2050 and 5.3% by 2091-2100 under RCP2.6 and by 4.4% by 2041-2050 and 16.0% by 2091-2100 under RCP8.5. These rates could be mitigated by the arrival of species that are currently at their northern limit in the English Channel. According to the classic approach (i.e. considering only BEMs and sea surface temperature and salinity), species richness was projected to decline by a mean of 9.3% by 2041-2050 and 6.5% by 2091-2100 under RCP2.6 and by 11.6% by 2041-2050 and 40.8% by 2091-2100 under RCP8.5. Thus, not considering the vertical temperature or salinity gradient or the type of seafloor resulted in higher projections of species loss (Figure 5). Moreover, considering both scenarios and both periods, the classic approach projected lower range loss (by a mean of 32%) for eight species (mainly benthic or demersal) but projected higher range loss (by a mean of 21%) for 11 species (also mainly benthic or demersal) (Figure 6). The proposed framework and the classic approach yielded similar projections for eight species, mainly benthopelagic (e.g. Pollachius pollachius, Trisopterus minutus, Trisopterus luscus, Gadus morhua).

Figure 5. Projected species richness anomalies between current climate conditions and projections for both periods (2041-2050 and 2091-2100), both scenarios (RCP2.6 and RCP8.5) and both approaches (the proposed framework with both vertical climatic conditions and habitat and the classic approach with surface climatic conditions only) Figure 6. Projected range loss for both periods (2041-2050 and 2091-2100) and both scenarios (RCP2.6 and 8.5) according to the proposed framework (with both vertical climatic conditions and habitat) and the classic approach (with surface climatic conditions only). Plots show species for which a range loss was projected for at least one period, scenario or method. The framework's projections of the 46 species revealed that the classic approach (i.e. considering only BEMs without vertical climatic conditions or seafloor type) projected higher loss of species richness. In contrast, no general trend was observed in the projected range loss for individual species. Nevertheless, the classic approach may have projected higher range losses more often than it predicted lower range losses, and mainly for benthic or demersal species, which depend on a specific type of seafloor. The classic approach projected range loss of benthic and demersal species due to the loss of suitable temperature conditions at the surface. However, since temperature fluctuates more at the surface than at the bottom, temperature conditions could remain suitable at the depth range of the species modelled, thus explaining why the classic approach projected higher range losses. Our results highlight that ignoring the vertical gradient of climatic conditions, even in shallow seas, and the habitat filter leads to inconsistent results. This is even more critical in ecosystems like the EEC that are strongly structured by the seafloor substrate [START_REF] Dauvin | The food web in the lower part of the Seine estuary: a synthesis of existing knowledge[END_REF][START_REF] Garcia | Assessment of benthic ecosystem functioning through trophic web modelling: the example of the eastern basin of the English Channel and the Southern Bight of the North Sea[END_REF]. Similar studies exist for marine areas next to the EEC, but none used such a complete framework. Moreover, they are based on the outdated IPCC SRES scenarios and consider only a few species (e.g. only 4 benthic species in the western English Channel; [START_REF] Rombouts | Potential changes in benthic macrofaunal distributions from the English Channel simulated under climate change scenarios[END_REF])) or only benthic species (75 benthic species in the North Sea; [START_REF] Weinert | Modelling climate change effects on benthos: Distributional shifts in the North Sea from 2001 to 2099[END_REF])), making comparison of results difficult. There are several research perspectives to refine marine species distribution models. Among them, incorporating mechanistic information such as species' physiological response to abiotic conditions (e.g., [START_REF] Feng | Physiology in ecological niche modeling: using zebra mussel's upper thermal tolerance to refine model predictions through Bayesian analysis[END_REF]) is a promising avenue to increase the realism of predicted responses to environmental conditions. In addition, marine species can have complex life cycles with phases that may require distinct models, and thus may be good candidates for infra-specific models (Smith et al. 2018). Another important perspective to consider is the biotic interactions in the process of model building, for example by constructing joint species distribution models.

Conclusion

Our results reinforce the need to consider species' positions in the water column and their dependence on a given type of seafloor to avoid biased predictions of species distributions or even misleading results (e.g. [START_REF] Duffy | Explicitly integrating a third dimension in marine species distribution modelling[END_REF]. In this paper, using the EEC as a case study, we went further and demonstrated that the simple and classic approach that considers sea surface climatic conditions as a proxy of bottom conditions and ignores habitat variables predicts higher species loss. Species distribution maps are informative, and many atlases contain them. However, they remain descriptive and show only current ranges (e.g. CHARM Atlas of the English Channel; [START_REF] Martin | The Channel habitat atlas for marine resource management (CHARM): an aid for planning and decision-making in an area under strong anthropogenic pressure[END_REF], Carpentier et al. 2009). As assessing species' future distributions is increasingly required for conservation biology and fisheries management, routinely generating maps using a relevant methodological framework would be useful. We developed a complete and operational framework, from data acquisition and preparation to a twostage modelling procedure, that considers habitat type, temperature and salinity at species' depth ranges. It creates maps illustrating current and projected ranges as well as response curves for each species. Our script can be adapted for individual needs (e.g., integrating other maps of climatic variables in the same data format) and is suitable for non-modeller ecologists to easily and widely address issues related to changes in species ranges due to climate change.

Acknowledgment

  Species occurrences used in BEMs can be downloaded from five global biogeographic databases:  Ocean Biogeographic Information System (OBIS): http://www.iobis.org/  Global Biodiversity Information Facility (GBIF): http://www.gbif.org/  iNaturalist: http://www.inaturalist.org  VertNet (vertebrate biodiversity networks): http://vertnet.org/  Ecoengine (UC Berkeley's Natural History Data): https://ecoengine.berkeley.edu/

Figure 2 .

 2 Figure 2. Computational framework and data processing steps 313 314 315

Figure 4 .

 4 Figure 4. Current and projected species richness in the Eastern English Channel for a pool of the 46 most abundant species. (a) current, (b) projection under RCP2.6 by 2041-2050, (c) projection under RCP2.6 by 2091-2100, (d) projection under RCP8.5 by 2041-2050 and (e) projection under RCP8.5 by 2091-2100.3.2 Classic approach vs. the proposed frameworkAccording to the complete approach developed in our methodological framework (namely the proposed framework), current species richness was projected to decline by a mean of 3.6% by 2041-2050 and 5.3% by 2091-2100 under RCP2.6 and by 4.4% by 2041-2050 and 16.0% by 2091-2100 under RCP8.5. These rates could be mitigated by the arrival of species that are currently at their northern limit in the English Channel. According to the classic approach (i.e. considering only BEMs and sea surface temperature and salinity), species richness was projected to decline by a mean of 9.3% by 2041-2050 and 6.5% by 2091-2100 under RCP2.6 and by 11.6% by 2041-2050 and 40.8% by 2091-2100 under RCP8.5. Thus, not considering the vertical temperature or salinity gradient or the type of seafloor resulted in higher projections of species loss (Figure5). Moreover, considering both scenarios and both periods, the classic approach projected lower range loss (by a mean of 32%) for eight species (mainly benthic or demersal) but projected higher range loss (by a mean of 21%) for 11 species (also mainly benthic or demersal) (Figure6). The proposed framework and the classic approach yielded similar projections for eight species, mainly benthopelagic (e.g. Pollachius pollachius, Trisopterus minutus, Trisopterus luscus, Gadus morhua).

  This work is directed by the University of Caen and France Energies Marines (FEM) and benefits from FEM and State financing managed by the National Research Agency under the Investments for the Future program bearing the reference ANR-10-IED-0006-12 References Alabia ID, Saitoh S-I, Igarashi H, Ishikawa Y, Usui N, Kamachi M, Awaji T, Seito M (2015) Future projected impacts of ocean warming to potential squid habitat in western and central North Pacific. ICES J Mar Sci 73:1343-1356. Albouy C, Guilhaumon F, Leprieur F, Ben Rais Lasram F, Somot S, Aznar R, Velez L, Le Loc'h, F, Mouillot D (2013) Projected climate change and the changing biogeography of coastal Mediterranean fishes. J Biogeogr 40:534-547. Albouy C, Velez L, Coll M, Colloca F, Le Loc'h F, Mouillot D, Gravel D (2014) From projected species distribution to food-web structure under climate change. Glob Change Biol 20:730-741.


Sarah Dellière

Emmanuel Dudoignon

Sofiane Fodil

Sebastian Voicu

Magalie Collet

Pierre-Antoine Oillic

Maud Salmona

François Dépret

Théo Ghelfenstein-Ferreira

Benoit Plaud

Pierre- Antoine Oillic

Benjamin Chousterman

Stéphane Bretagne

Elie Azoulay

Alexandre Mebazaa

Bruno Megarbane

MD Alexandre Alanio
email: alexandre.alanio@pasteur.fr

Risk factors associated with Covid-19-associated pulmonary aspergillosis in ICU patients: a French multicentric retrospective cohort

Keywords: Aspergillus, COVID-19, azithromycin, SARS-CoV-2, Critical Care, Coronavirus, Corticosteroids

à 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

Although pulmonary invasive fungal disease is typically described in the immunocompromised host, invasive pulmonary aspergillosis (IPA) has been increasingly reported in critically ill patients including patients without classical risk factors of immunosuppression [START_REF] Tudesq | Invasive Pulmonary Aspergillosis in Nonimmunocompromised Hosts[END_REF]. In acute respiratory distress syndrome (ARDS) patients, ~12.5% of the patients had IPA as shown by random post-mortem histopathological examination of lung tissue [START_REF] De Hemptinne | ARDS: a clinicopathological confrontation[END_REF]. Coronavirus disease 2019 (COVID-19)-associated pulmonary aspergillosis (CAPA) has been increasingly reported [START_REF] Alanio | Prevalence of putative invasive pulmonary aspergillosis in critically ill patients with COVID-19[END_REF][START_REF] Bartoletti | Epidemiology of invasive pulmonary aspergillosis among COVID-19 intubated patients: a prospective study[END_REF][START_REF] White | A national strategy to diagnose COVID-19 associated invasive fungal disease in the ICU[END_REF]. Whether the use of immunomodulatory therapies such as corticosteroids prescribed to dampen detrimental inflammatory response and antibiotics to treat and/or prevent bacterial superinfections is responsible for increased susceptibility of COVID-19 patients to pulmonary invasive fungal disease remains incompletely studied [START_REF] Schauwvlieghe | Invasive aspergillosis in patients admitted to the intensive care unit with severe influenza: a retrospective cohort study[END_REF][START_REF] Erjavec | Trends in invasive fungal infections, with emphasis on invasive aspergillosis[END_REF]. The aim of our study was to evaluate the incidence of IPA and the risk factors associated with IPA in severe COVID-19 patients admitted to the intensive care unit (ICU), and evaluate the impact of IPA on patient's outcome.

J o u r n a l P r e -p r o o f

Method

We conducted a retrospective observational cohort study. All successive COVID-19 patients admitted to the four ICUs of our two university hospitals between March 15 th and May 1 st 2020 with a positive SARS-CoV-2 PCR (Cobas® SARS-CoV-2 Test, Roche) and ≥1 respiratory sample (bronchoalveolar lavage (BAL), tracheal aspirate, sputum) sent to the mycology department were included (Figure 1). Of note, the 27 first patients included were previously partially analyzed [START_REF] Alanio | Prevalence of putative invasive pulmonary aspergillosis in critically ill patients with COVID-19[END_REF]. On respiratory sample, culture, galactomannan (GM) (BAL only) and Aspergillus quantitative polymerase chain reaction (qPCR) were systematically performed. In concomitantly received blood sample, GM, ß-D-glucan (BDG) and Aspergillus qPCR were performed on serum/plasma. Patients were classified as probable IPA according to European Organization for Research and Treatment of Cancer (EORTC) and the Mycoses Study Group Education and Research Consortium (MSGERC) consensus criteria in immunocompromised patients [START_REF] Donnelly | Revision and Update of the Consensus Definitions of Invasive Fungal Disease From the European Organization for Research and Treatment of Cancer and the Mycoses Study Group Education and Research Consortium[END_REF] and according to the consensus case definition proposal for influenza/COVID-19-associated pulmonary aspergillosis (CAPA) in ICU patients otherwise (Table S1) [START_REF] Verweij | Review of influenza-associated pulmonary aspergillosis in ICU patients and proposal for a case definition: an expert opinion[END_REF]. An extensive list of clinical data was collected as part of the initial protocol (Table 1). The cumulative dose of corticosteroids, azithromycin and ß-lactams were determined as the total dose of drug received prior to the day of sampling. Prescription of azithromycin >1500mg and ß-lactams >3 days were predefined as exposure variables. Azithromycin was systematically prescribed before or on the day of admission to the ICU, except for 2 patients for whom it was introduced at day 4 and 5, respectively. If no fungal infection was diagnosed the latest sample was used. Corticosteroids doses were quantified as dexamethasone-equivalent [START_REF] Lionakis | Glucocorticoids and invasive fungal infections[END_REF].

Culture of respiratory specimens were performed, as previously described [START_REF] Alanio | New therapeutic strategies for invasive aspergillosis in the era of azole resistance: how should the prevalence of azole resistance be defined?: Table 1[END_REF]. For Aspergillus qPCR, DNA was extracted from 1mL of plasma or from bead-beaten pellet of the respiratory sample and resuspended in 1000µL of DNA-free water using the Qiasymphony DSP virus/Pathogen Mini kit (Qiagen,) and a QIAsymphony apparatus (Qiagen). PCR assay was previously reported [START_REF] Alanio | Circulating Aspergillus fumigatus DNA Is Quantitatively Correlated to Galactomannan in Serum[END_REF]. GM and BDG detection were performed using Platelia Bio-Rad kit (BioRad Laboratories) and Fungitell assay (Cape Cod Diagnostics) according to the manufacturer respectively.

Statistics

J o u r n a l P r e -p r o o f interquartiles [Q1-Q3] as appropriate. Univariate analyses were performed to assess an association between clinical factors and IPA using Fisher's exact, Chi-2 and Wilcoxon tests as appropriate. Odds ratios (OR) with 95%-confidence intervals (IC95%) were calculated for each significant variable based on univariate logistic regression. All analyses were performed using R software, version 3.5.3 (http://www.r-project.org).

Ethical statements

Our institutional ethics committee approved the study (IDRCB, 2020-A00256-33; CPP, 11-20-20.02.04.68737).

J o u r n a l P r e -p r o o f 6

Results

A total of 366 patients with positive SARS-CoV2 qPCR were admitted to the four intensive care units between March 15 th and May 1 st 2020 among which 246 were intubated and mechanically ventilated (Figure 1). The mycology department received 193 respiratory samples from 108 patients, whose conditions deteriorated despite appropriate initial care. Patient characteristics are described in When comparing patients who developed probable IPA (n=21) or not (n=87), no significant differences were observed regarding general population characteristics and severity upon admission (Table 1). Prescription of hydroxychloroquine (n=34) did not differ between both groups. Administration of azithromycin for more than 3 days (cumulative dose ≥1500mg) was associated with probable IPA (OR, 3.1; IC95%, 1.1-8.5; p=0.025) (Figure S1). Of note, 34 patients received azithromycin, which was prematurely discontinued on day 1 or day 2 in 8 patients because of QT interval prolongation. Administration of high-dose corticosteroids was not significantly associated with IPA (11.5 vs. 28.6%; p=0.08), although the cumulative dose ≥100 mg tended to be higher among IPA patients (OR, 3.7; IC95%, 1.0-9.7). Details for the incidence of IPA among patients who received azithromycin and/or corticosteroids is available in Table S3. Mortality was significantly higher in the probable IPA group (15/21, 71.4% vs. 32/87, 36.8%; p<0.01). Further details on each IPA patients are available in Table S4.

Discussion

J o u r n a l P r e -p r o o f

In our study, the incidence of IPA in COVID-19 patients was 5.7% in all ICU patients and 8.5% in those mechanically ventilated yet may be underestimated considering only patients with clinical worsening were tested.

A cumulative azithromycin dose ≥1500mg was associated with IPA. Azithromycin have in vitro antiviral effect and is a broad spectrum antibiotic with immunomodulatory properties which could have also both prevented bacterial superinfections and reduced inflammation [START_REF] Oldenburg | Azithromycin for severe COVID-19[END_REF]. A recent metaanalysis found an increased mortality when hydroxychloroquine is associated to azithromycin [START_REF] Fiolet | Effect of hydroxychloroquine with or without azithromycin on the mortality of COVID-19 patients: a systematic review and meta-analysis[END_REF].

Azithromycin-related impact on the risk of secondary infections have been only incompletely studied.

Azithromycin has been shown to decrease serum interleukin-6 and induce delayed down-regulation of neutrophil oxidative burst and increased apoptosis up to 28 days after 3 doses of azithromycin (i.e. 1500 mg) in humans [START_REF] Culić | Azithromycin modulates neutrophil function and circulating inflammatory mediators in healthy human subjects[END_REF]. Neutrophils and oxidative burst represent the first and most important immune system barrier against aspergillosis [START_REF] Latgé | Aspergillus fumigatus and Aspergillosis in 2019[END_REF]. Furthermore, azithromycin may promote Aspergillus colonization by altering the lung microbiome [START_REF] Dickson | Macrolides, inflammation and the lung microbiome: untangling the web of causality[END_REF].

Corticosteroids are known to increase susceptibility to invasive fungal disease due to complex quantitative and qualitative immune deregulation [START_REF] Lionakis | Glucocorticoids and invasive fungal infections[END_REF]. High dose corticosteroid, although not precisely defined, was previously found to be associated with CAPA [START_REF] Bartoletti | Epidemiology of invasive pulmonary aspergillosis among COVID-19 intubated patients: a prospective study[END_REF][START_REF] White | A national strategy to diagnose COVID-19 associated invasive fungal disease in the ICU[END_REF]. Although not significantly associated to IPA in our study, probably because of insufficient statistical power, a trend was observed after a cumulative dose of ≥100 mg dexamethasone-equivalent (OR, 3.7; IC95, 1.0-9.7), Recent studies such as the RECOVERY trial and various meta-analyses showed that corticosteroid administration is beneficial in COVID-19 patients requiring hospitalization [START_REF] Sterne | Association Between Administration of Systemic Corticosteroids and Mortality Among Critically Ill Patients With COVID-19[END_REF]. Interestingly, the cumulative dose of dexamethasone using the RECOVERY trial regimen does not exceed 60 mg.

Although no severity score or variables, including oxygenation parameters, were found to be associated to IPA occurrence, we cannot rule out the extension of the lesion to be an associated risk factor. Indeed, the extension of lung lesions quantified by computed tomography was a predictor of COVID-19 severity and early death [START_REF] Ruch | CT lung lesions as predictors of early death or ICU admission in COVID-19 patients[END_REF]. Susceptibility to IPA in previously immunocompetent critically ill patients is most likely multifactorial. In ARDS patients, epithelial damage, impaired mucocilliary clearance and temporary immune deregulation, starting with excess release of danger-associated molecular patterns (DAMPs) J o u r n a l P r e -p r o o f secondary to COVID-19 damages, may be initiating factors [START_REF] Arastehfar | COVID-19 Associated Pulmonary Aspergillosis (CAPA)-From Immunology to Treatment[END_REF]. The addition of known or suspected risk factors, such as corticosteroid or azithromycin, further inhibiting neutrophils and innate immune response may tilt the balance in favor of IPA development. The risk of IPA associated with corticosteroids compellingly depends on its cumulative dose although cut-offs are not clearly defined and depends on the underlying host factor [START_REF] Donnelly | Revision and Update of the Consensus Definitions of Invasive Fungal Disease From the European Organization for Research and Treatment of Cancer and the Mycoses Study Group Education and Research Consortium[END_REF]. Our findings raise questions regarding the possible connection between azithromycin use and the observed increased susceptibility to IPA, which needs to be further explored.

Funding

No external funding was received for the present study

Obesity n (%)

J o u r n a l P r e -p r o o f 14 Figure 1 .

 141 Figure 1. Study flowchart.

Table S2 .

 S2 Male/female sex ratio was 4.4 and median age was 61 years. Twenty-one patients developed probable IPA according to CAPA criteria stricto sensu (n=19) and EORTC/MSGERC definitions (n=2; one solid organ transplant recipient and one myeloma patient). Overall, incidence was 5.7% (21/366) in severe COVID-19 patients admitted to the ICU and 8.5% (21/246) in those mechanically ventilated. IPA incidence in patients whose conditions worsened despite appropriate care was 19.4% (21/108). The median times from symptom onset to IPA diagnosis and from ICU admission to IPA diagnosis were 16 days (10-23) and 6 days (1-15), respectively.

Table 1 .

 1 Comparison of severe COVID-19 patients with and without probable invasive pulmonary

	251
	252	aspergillosis
	253
			Total	Without IPA	With IPA
						OR	CI95%	p-value
			(n=108)	(n=87)	(n=21)
		Male n (%)	88 (81.5)	72 (82.8)	16 (76.2)	0.7	0.2-2.1	-
	Age median [Q1-Q3]	62 [56-68]	62 [56-68]	63 [56.75-68.25]	-	-	0.63∘
		Mechanical
			105 (97.2)	85 (97.7)	20 (95.2)	0.5	0.04-5.3
		ventilation n (%)
				COVID risk factors
		HTA n (%)	64 (59.3)	50 (57.5)	14 (66.7)	1.5	0.5-4.0	-
		Diabetes n (%)	40 (37.0)	31 (

  Bronchoalveolar lavage; BMI: body mass index; CI95%: 95% confidence interval; IPA: invasive pulmonary aspergillosis; LOS: length of stay; OR: odd ratio; PMN: polymorphonuclear; SAPSII, Simplified Acute Physiology Score II; SOFA, Sequential Organ Failure Assessment; ∘Wilcoxon test

	48H n (%) Eculizumab n(%)	6 (5.6)	4 (4.6)		2 (9.5)	2.2	0.4-12.8
	Creatininemia Tocilizumab n(%)	4 (3.7)	2 (2.3)		2 (9.5)	4.5	0.6-33.8
		103.34 (74.01) 92.64 ( 47.19) 149.85 (132.96)	-	-	0.08∘
	(mg/dL) mean (SD)	Therapy with cumulative dose before sampling
	D-dimer median [Q1-Azithromycin >1500	2395 [1193-	2325 [1163-	2515 [1610-
		26 (24.1)	17 (19.5)		9 (42.9)	-3.1	-1.1-8.5	0.63∘
	Q3] mg total dose n (%)	4635]	4563]		10917]
	Dexamethasone	755.11
	LDH mean (SD)	16 (14.8)	759.49 (303.61) 740.20 (350.52) 10 (11.5) 6 (28.6)	-3.1	-1.0-9.8	0.80∘
	>1000mg n (%)	(312.15)
	SAPS2 mean (SD) Any ß-lactam > 3	39.93 (14.40)	40.4 (14.6)		38.1 (13.8)	-	-	0.58∘
	) 28 [25-29] Other patient characteristics 28 [26-32] 3 (3.4) 2 (9.5) J o u r n a l P r e -p r o o f 0.6 BMI median [Q1-Q3] 28 [25-31] -Asthma n (%) 5 (4.6) 2.9 SOFA mean (SD) 6.02 (3.79) 5.8 (3.6) 7.1 (4.5) -Severity during hospitalization Nadir PaO2/FiO2 mean (SD) 79.75 (37.21) 81.54 (39.01) 72.50 (28.40) -ECMO n (%) 10 (9.3) 9 (10.3) 1 (4.8) 0.4 Renal replacement therapy n (%) 38 (35.2) 30 (34.5) 8 (38.1) 1.2 Vasopressors n (%) 89 (82.4) 70 (80.5) 19 (90.5) 2.3 Specific COVID therapy Lopinavir-ritonavir 16 (14.8) 10 (11.5) 6 (28.6) 3.1 days n (%) 90 (83.3) 74 (85.1) 16 (76.2) 0.6 Respiratory sample characteristics PaO2/FiO2 at sampling mean (SD) 173.69 (91.70) 173.18 (96.05) 175.66 (74.40) -% BAL macrophages mean (SD) 31.23 (21.94) 31.00 (22.76) 32.17 (20.21) -% BAL PMN mean (SD) 47.37 (30.92) 47.08 (31.99) 48.50 (28.89) -% BAL Lymphocytes mean (SD) J o u r n a l P r e -p r o o f 20.63 (18.94) 21.17 (19.10) 18.50 (19.92) -	0.1-2.9 -0.5-18.9 --0.1-3.6 0.4-3.1 0.5-10.9 0.9-9.8 0.2-1.8 ----	-0.70∘ -0.28∘ 0.50∘ ---0.612∘ 0.836∘ 0.959∘ 0.795∘
	COPD n (%) n(%)	2 (1.9)	2 (2.3) Outcome	0 (0.0)	0.8	0.04-17.2	-
	Immunocompromised Hydroxychloroquine Mortality n (%)	47 (43.5)	32 (36.8)		15 (71.4)	4.3	1.5-12.1	<0.01 ◊
		10 (9.3) 34 (31.5)	8 ( 9.2) 27 (31.0)		2 (9.5) 7 (33.3)	0.6 1.1	0.1-2.9 0.4-3.1	-
	patient n (%) n(%) LOS days mean (SD) 24.33 (18.88) 25.13 (19.18)	21.05 (17.60)	-	-	0.313∘
	Long term Azithromycin +
		11 (10.2)	8 ( 9.2)		3 (14.3)	1.6	0.4-6.8	-
	corticosteroids n (%) Hydroxychloroquine	29 (26.9)	22 (25.3)		7 (33.3)	1.4	0.5-4.1
	n(%)		Severity at admission
	Immunoglobulins	173.47
	PaO2/FiO2 mean (SD)	3 (2.8)	169.63 (125.85) 187.79 (114.74) 3 (3.4) 0 (0.0)	-0.6	-0.03-11.3	0.34∘
	n(%)	(123.19)
	Vasopressors in first Sarilumab n(%)	65 (60.2) 1 (0.9)	52 (59.8) 1 (1.1)		13 (61.9) 0 (0.0)	1.1 4.3	0.4-2.9 0.3-71.8

BAL:

Acknowledgement

We thank Dr Pierre Gazeau for providing follow-up information regarding patients transferred to Brest University Hospital. We thank all staff, nurses, lab technicians who were essential to patients care in intensive care and medical mycology from Saint Louis Hospital and Lariboisière Hospital.

Conflict of interest

The authors declare no conflict of interest related to the content of the present study

Authors contribution

Writing-Original draft: SD and AA. Writing-Review & Editing: All. Conceptualization: SD, AA and SB. Investigation: SD, ED, SF, SV, MC. Data curation: SD, ED, SF, SV, MC, TFG. Formal analysis: PAO, MS. Visualization: SD, AA, BM, EA, AM. Supervision: AA, BM, EA, AM.

J o u r n a l P r e -p r o o f


Alexandre Alanio
email: alexandre.alanio@pasteur.fr

Sarah Dellière

Sebastian Voicu

Stéphane Bretagne

Bruno Mégarbane

The presence of Pneumocystis jirovecii in critically ill patients with COVID-19

Keywords: COVID-19, pneumocystosis, Pneumocystis jirovecii

1

Letters to the Editor

Dear Editor,

We read with interest the recent review on co-infections in coronavirus disease-2019 (COVID-19) patients [START_REF] Lansbury | Co-infections in people with COVID-19: a systematic review and meta-analysis[END_REF] and believe that fungal co-infections as evaluated from selected studies are underestimated.

Severe acute respiratory syndrome coronavirus-2 (SARS-CoV2) is still spreading pandemically.

Approximately 5-10% of COVID-19 patients may require intensive care unit (ICU) management and 30% may develop secondary pneumonia without identified etiology. [START_REF] Huang | Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China[END_REF] Hospital-acquired bacterial or fungal superinfections, as described in critically ill patients with Influenza virus, can be suspected. [START_REF] Schauwvlieghe | Invasive aspergillosis in patients admitted to the intensive care unit with severe influenza: a retrospective cohort study[END_REF] Since pneumocystosis is usually reported in patients with T-cell immunodepression, [START_REF] Alanio | ECIL guidelines for the diagnosis of Pneumocystis jirovecii pneumonia in patients with haematological malignancies and stem cell transplant recipients[END_REF] less attention has been paid to Pneumocystis jirovecii in non-immunocompromised ICU patients although it accounts for 7% of the co-infections reported in those admitted with Influenza. [START_REF] Beumer | Influenza virus and factors that are associated with ICU admission, pulmonary co-infections and ICU mortality[END_REF] Interestingly, COVID-19 patients may develop lymphocytopenia and acute respiratory distress syndrome (ARDS) requiring adjunctive steroids and/or immunomodulatory therapies, well-known susceptibility factors for developing pneumocystosis. [START_REF] Beumer | Influenza virus and factors that are associated with ICU admission, pulmonary co-infections and ICU mortality[END_REF] We designed this observational cohort study to investigate the prevalence of P. jirovecii acid nucleic detection in respiratory specimens sampled to identify co-infections in COVID-19 patients in the ICU.

All consecutive patients admitted to the ICU between 2020/03/15 and 2020/05/01 with a positive SARS-CoV-2 PCR (Cobas ® SARS-CoV-2 Test, Roche, France) and ≥1 respiratory sample (bronchoalveolar lavage (BAL), tracheal aspirate, sputum) sent to the mycology department. This study was part of the COVID-ICU and French COVID-19 cohort registries. Our institutional ethics committee approved the study (IDRCB, 2020-A00256-33; CPP, 11-20-20.02.04.68737). When possible, signed informed consent was obtained from the patients or the next of kin.

Whenever possible, bronchoalveolar lavage (BAL) was performed in the middle lobe by a trained pneumologist in the ICU using at least 120 ml saline (yield, ~50%). Upon reception, specimens including BAL fluids and dTT-treated aspirations (dTT 1X at 37°C for 15min) were centrifuged, suspended in 200µL of water and submitted to extraction (whole nucleic acids extraction) using the GeneLead-VIII extractor-thermocycler™ (Precision System Science, Japan). P. jirovecii reverse transcriptase quantitative PCR (RTqPCR) was performed to amplify mtSSU and mtLSU RNA and DNA of P. jirovecii using the new R-DiaPnJ kit™ (Diagenode, Belgium). Serum -D-glucan was tested using the Fungitell kit™ (Cape Cod Inc, US) as recommended by the manufacturer. Data are presented as median [25 th -75 th percentiles] or percentages as appropriate. Comparisons were performed using Mann-Whitney or exact Fisher tests as required. P-values ≤0.05 were considered as significant.

One hundred-and-eight successive HIV-negative COVID-19 patients (Male/Female sex ratio, 4.4; age, 62 years [56-68]) with the usual risk factors for severe COVID-19 presentation were included (Table1).

All except three patients were intubated on admission. Thirty-four patients (31.4%) who developed ARDS received at least one day of corticosteroids before BAL sampling. Respiratory samples included 80 BALs (74.1%), 22 tracheal aspirates (20.4%), 4 sputa (3.7%) and two bronchial aspiration fluids (1.9%).

In 10/108 patients (9.3%), P. jirovecii RTqPCR was positive. Median delay between sampling and ICU admission was 2days [1-2]. The median quantitative cycle value was 32. 6 [30.8-34.7]. Serum -D-glucan were measured in nine patients and was negative (>80pg/mL) in seven patients. Clinical characteristics of the patients carrying P. jirovecii did not significantly differ from the other patients except for lower plasma D-dimer (1,270ng/mL [750-2,390] vs 2,610ng/mL [1,405-4,700], P=0.03) and more frequent lopinavir/ritonavir administration (40.0% vs 12.2%, P=0.04), while long-term corticosteroid prescription tended to be more frequent (30.0% vs 8.2%, P=0.06). Of note, among our 10 P. jirovecii carriers, five concomitantly met the criteria for COVID-19-associated pulmonary aspergillosis. [START_REF] Alanio | Prevalence of putative invasive pulmonary aspergillosis in critically ill patients with COVID-19[END_REF] Out of these patients, four (40%) received co-trimoxazole as prophylaxis (80/400mg once daily) whereas six including four who rapidly improved did not. One co-trimoxazole-treated and two nontreated patients died while the seven remaining patients were discharged. Mortality was similar in both groups.

We found an unexpectedly high proportion of critically ill COVID-19 patients detected with P. jirovecii (10/108 patients; 9.3%), similarly to previous findings in influenza patients (3/45; ~7%). [START_REF] Beumer | Influenza virus and factors that are associated with ICU admission, pulmonary co-infections and ICU mortality[END_REF] The presence of P. jirovecii in the healthy adult population has been measured using oropharyngeal wash samples obtained by gargling and examined by conventional or nested PCR methods. [START_REF] Medrano | Pneumocystis jirovecii in general population[END_REF] However, experts agree that the reported prevalence (~20%) has been overestimated due to technical issues such as contamination with amplicons responsible for false positives. [START_REF] Alanio | ECIL guidelines for the diagnosis of Pneumocystis jirovecii pneumonia in patients with haematological malignancies and stem cell transplant recipients[END_REF] In our center managing almost exclusively immunocompromised patients, prevalence of qPCR-positive respiratory specimens with fungal load as low as in our COVID-19 patients, is ~13% (unpublished data), as reported elsewhere. [START_REF] Schauwvlieghe | Invasive aspergillosis in patients admitted to the intensive care unit with severe influenza: a retrospective cohort study[END_REF] COVID-19 patients mostly exhibited marked lymphopenia and alterations in lymphocyte functions, [START_REF] Liu | Longitudinal characteristics of lymphocyte responses and cytokine profiles in the peripheral blood of SARS-CoV-2 infected patients[END_REF] likely explaining the high-rate of P. jirovecii detection.

Since serum -D-glucan is advocated in pneumocystosis diagnosis, [START_REF] Alanio | ECIL guidelines for the diagnosis of Pneumocystis jirovecii pneumonia in patients with haematological malignancies and stem cell transplant recipients[END_REF] we measured its concentrations in four of our five P. jirovecii RTqPCR-positive patients and obtained low values (<120pg/mL) in accordance with the low nucleic acids fungal loads in the lung alveoli. [START_REF] Karageorgopoulos | Accuracy of β-D-glucan for the diagnosis of Pneumocystis jirovecii pneumonia: a meta-analysis[END_REF] Of note, in two out of our nine tested P. jirovecii RTqPCR-negative patients, higher B-D-glucan concentrations (450 and 500pg/ml) lead to the diagnosis of pulmonary aspergillosis, another fungal infection of risk in COVID-19 patients. [START_REF] Alanio | Prevalence of putative invasive pulmonary aspergillosis in critically ill patients with COVID-19[END_REF] Although a recent meta-analysis questioned its sensitivity in non-HIV patients, 10 -D-glucan has been widely used to rule out pneumocystosis because of its high negative predictive value. This finding may support the hypothesis that our patients were carrying P. jirovecii, yet not being infected per se. Thus, although interesting in the context of invasive fungal infections diagnosis, serum -D-glucan should be interpreted with caution when excluding the diagnosis of pneumocystosis.

Here, four out of ten P. jirovecii RTqPCR-positive patients received co-trimoxazole as prophylactic regimen, based on the treating physician's decision. Whether a positive result should be an indication to consider administering co-trimoxazole, at least at prophylactic dosage in COVID-19 patients remains questionable.

Our study limitations include the relatively small number of patients, the bi-center setting, and the short study period. However, to the best of our knowledge, this is the first study evaluating the prevalence of P. jirovecii in COVID-19 patients. Because we focused on critically ill COVID-19 patients, P. jirovecii prevalence in less severe patients remains to be determined.

In conclusion, an unexpectedly high proportion of P. jirovecii-positive pulmonary samples is observed in critically ill COVID-19 patients. Based on our findings, we advocate systematically searching for P. jirovecii in deep respiratory specimens in these patients. We believe that this strategy may be useful in limiting enhanced inflammation due to the presence of P. jirovecii in the lung and avoiding inter-patient P. jirovecii transmission.

  10. Del Corpo O, Butler-Laporte G, Sheppard DC, Cheng MP, McDonald EG, Lee TC. Diagnostic accuracy of serum (1-3)-β-D-glucan for Pneumocystis jirovecii pneumonia: a systematic review and meta-analysis. Clin Microbiol Infect 2020;S1198-743X(20)30301-3. doi:10.1016/j.cmi.2020.05.024

Table 1

 1 Characteristics of hundred-and-eight critically ill COVID-19 patients according to Pneumocystis jirovecii detection in the respiratory samples. Data are presented as percentages or medians [25 th -75 th percentiles].

	Comparisons were performed using Mann-Whitney (∘) or Fisher exact tests (*), as appropriate.
		Total	No detection of P. jirovecii	Detection of P. jirovecii	P
		(N=108)	(N=98)	(N=10)
	Male gender, N (%)	88 (81.5%)	80 (81.6%)	8 (80.0%)	1*
	Age (years)	62 [56-68]	62 [56-68]	59 [46-68]	0.40∘
	COVID-19 risk factors
	Past hypertension, N (%)	64 (59.3%)	58 (59.2%)	6 (60.0%)	1*
	Diabetes, N (%)	40 (37.0%)	37 (37.8%)	3 (30.0%)	0.74*
	Obesity, N (%)	35 (32.4%)	32 (32.7%)	3 (30.0%)	1*
	Coronary disease, N (%)	15 (13.9%)	14 (14.3%)	1 (10.0%)	1*
	Body-mass index (kg/m 2 )	28 [25-31]	28 [25-31]	28 [27-32]	0.61∘
	Other remarkable comorbidities
	Asthma, N (%)	5 (4.6%)	4 (4.1%)	1 (10.0%)	0.39*
	Chronic obstructive pulmonary disease, N (%)	2 (1.9%)	2 (2.0%)	0 (0.0%)	1*
	Immunocompromised patient, N (%)	10 (9.3%)	10 (10.2%)	0 (0.0%)	0.59*
	Long-term corticosteroids, N (%)	11 (10.2%)	8 (8.2%)	3 (30.0%)	0.06*
	Biological data of interest on admission
	PaO2/FiO2 (mmHg)	137 [83-247]	134 [83-239] 177 [108-253] 0.60∘

Table 1

 1 Characteristics of hundred-and-eight critically ill COVID-19 patients according to Pneumocystis jirovecii detection in the respiratory samples. Data are presented as percentages or medians [25 th -75 th percentiles].

	Comparisons were performed using Mann-Whitney (∘) or Fisher exact tests (*), as appropriate.
		Total	No detection of P. jirovecii	Detection of P. jirovecii	P
		(N=108)	(N=98)	(N=10)
	Male gender, N (%)	88 (81.5%)	80 (81.6%)	8 (80.0%)	1*
	Age (years)	62 [56-68]	62 [56-68]	59 [46-68]	0.40∘
	COVID-19 risk factors
	Past hypertension, N (%)	64 (59.3%)	58 (59.2%)	6 (60.0%)	1*
	Diabetes, N (%)	40 (37.0%)	37 (37.8%)	3 (30.0%)	0.74*
	Obesity, N (%)	35 (32.4%)	32 (32.7%)	3 (30.0%)	1*
	Coronary disease, N (%)	15 (13.9%)	14 (14.3%)	1 (10.0%)	1*
	Body-mass index (kg/m 2 )	28 [25-31]	28 [25-31]	28 [27-32]	0.61∘
	Other remarkable comorbidities
	Asthma, N (%)	5 (4.6%)	4 (4.1%)	1 (10.0%)	0.39*
	Chronic obstructive pulmonary disease, N (%)	2 (1.9%)	2 (2.0%)	0 (0.0%)	1*
	Immunocompromised patient, N (%)	10 (9.3%)	10 (10.2%)	0 (0.0%)	0.59*
	Long-term corticosteroids, N (%)	11 (10.2%)	8 (8.2%)	3 (30.0%)	0.06*
	Biological data of interest on admission
	PaO2/FiO2 (mmHg)	137 [83-247]	134 [83-239] 177 [108-253] 0.60∘

Table 1 2

 1 Serum BAL, bronchoalveolar lavage; ECMO, extracorporeal membrane oxygenation; ICU, intensive care unit; SOFA, Sepsis-6 related Organ Failure Assessment; SAPS II, Simplified Acute Physiology Score II.

	creatinine (µmol/L)	80 [64-111]	80 [63-111]	80 [67-104] 0.99∘
	Plasma D-dimer (ng/mL)	2,2395 [1,193-4,635]	2,610 [1,405-4,700]	1,270 [750-2,390]	0.03∘
	Serum lactate dehydrogenase (IU/L)	687 [540-901] 687 [568-903] 708 [436-893] 0.65∘
	Bronchoalveolar lavage characteristics
	% BAL macrophages	28 [15-46]	27 [14-42]	51 [49-55]	0.13∘
	% BAL polymorphonuclear cells	37 [26-81]	46 [26-81]	29 [18-32]	0.54∘
	% BAL Lymphocytes	13 [6-32]	14 [5-34]	13 [10-23]	0.81∘
	Specific anti-COVID-19 therapy, N (%)
	Azithromycin, N (%)	34 (31.5%)	30 (30.6%)	(40.0%)	0.72*
	Hydroxychloroquine, N (%)	34 (31.5%)	30 (30.6%)	(40.0%)	0.72*
	Hydroxychloroquine + Azithromycin, N (%)	29 (26.9%)	26 (26.5%)	(30.0%)	1*
	Lopinavir-ritonavir, N (%)	16 (14.8%)	12 (12.2%)	(40.0%)	0.04*
	Polyvalent immunoglobulins, N (%)	3 (2.8%)	3 (3.1%)	0 (0.0%)	1*
	Sarilumab, N (%)	1 (0.9%)	1 (1.0%)	0 (0.0%)	1*
	Eculizumab, N (%)	6 (5.6%)	4 (4.1%)	(20.0%)	0.10*
	Tocilizumab, N (%)	4 (3.7%)	4 (4.1%)	0 (0.0%)	1*
	Dexamethasone, N (%)	53 (49.1%)	46 (46.9%)	(70.0%)	0.19*
	Dexamethasone cumulative dose >100mg, N (%)	16 (14.8%)	15 (15.3%)	(10.0%)	1*
	Severity during hospitalization and outcome
	SAPS II on admission	37 [31-49]	38 [31-51]	34 [28-37]	0.16∘
	SOFA on admission	6 [3-8]	6 [3-8]	5 [2-7]	0.43∘

[START_REF] Schauwvlieghe | Invasive aspergillosis in patients admitted to the intensive care unit with severe influenza: a retrospective cohort study[END_REF][START_REF] Alanio | ECIL guidelines for the diagnosis of Pneumocystis jirovecii pneumonia in patients with haematological malignancies and stem cell transplant recipients[END_REF][START_REF] Beumer | Influenza virus and factors that are associated with ICU admission, pulmonary co-infections and ICU mortality[END_REF][START_REF] Alanio | Prevalence of putative invasive pulmonary aspergillosis in critically ill patients with COVID-19[END_REF][START_REF] Medrano | Pneumocystis jirovecii in general population[END_REF][START_REF] Liu | Longitudinal characteristics of lymphocyte responses and cytokine profiles in the peripheral blood of SARS-CoV-2 infected patients[END_REF] 6[START_REF] Schauwvlieghe | Invasive aspergillosis in patients admitted to the intensive care unit with severe influenza: a retrospective cohort study[END_REF][START_REF] Alanio | ECIL guidelines for the diagnosis of Pneumocystis jirovecii pneumonia in patients with haematological malignancies and stem cell transplant recipients[END_REF][START_REF] Beumer | Influenza virus and factors that are associated with ICU admission, pulmonary co-infections and ICU mortality[END_REF][START_REF] Alanio | Prevalence of putative invasive pulmonary aspergillosis in critically ill patients with COVID-19[END_REF][START_REF] Medrano | Pneumocystis jirovecii in general population[END_REF][START_REF] Liu | Longitudinal characteristics of lymphocyte responses and cytokine profiles in the peripheral blood of SARS-CoV-2 infected patients[END_REF] 5[START_REF] Huang | Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China[END_REF][START_REF] Schauwvlieghe | Invasive aspergillosis in patients admitted to the intensive care unit with severe influenza: a retrospective cohort study[END_REF][START_REF] Alanio | ECIL guidelines for the diagnosis of Pneumocystis jirovecii pneumonia in patients with haematological malignancies and stem cell transplant recipients[END_REF][START_REF] Beumer | Influenza virus and factors that are associated with ICU admission, pulmonary co-infections and ICU mortality[END_REF][START_REF] Alanio | Prevalence of putative invasive pulmonary aspergillosis in critically ill patients with COVID-19[END_REF][START_REF] Medrano | Pneumocystis jirovecii in general population[END_REF] 0.43∘

Acknowledgments

The authors would like to thank Mrs. Alison Good (Scotland, UK) for her helpful review of the manuscript.

Funding

None. The case investigations, analysis, and manuscript preparation were completed as part of official duties at the university hospital.

Declaration of Competing Interest

The authors declare that they have no competing interests.

Ethics approval and consent to participate

This study was part of the French COVID-19 cohort registry conducted by the REACTing consortium (REsearch and ACTion targeting emerging infectious diseases) and directed by INSERM (Institut national de la santé et de la recherche médicale) and ISARIC (International Severe Acute Respiratory and Emerging Infection Consortium). Our institutional ethics committee approved the study (N°, IDRCB, 2020-A00256-33; CPP, 11-20-20.02.04.68737).

Availability of data and materials

Drs. Alanio, Mégarbane and Bretagne conceived of and designed the study. Drs. Voicu, Azoulay and Mégarbane managed the patients. Drs. Alanio, Dellière and Bretagne performed the microbiological analysis. All authors acquired, analyzed the data and interpreted the results. Drs. Alanio, Mégarbane and Bretagne drafted the manuscript. All authors participated to the critical revision of the manuscript for important intellectual content.

Dr. Alanio has full access to all data and takes responsibility for the data integrity and its analysis accuracy.

Consent for publication

All the authors agree to publish.