Marcin Pietras
email: mpietras@man.poznan.pl

Marta Kolanowska

Marc-André Selosse

Claus Baessler

Quo vadis? Historical distribution and impact of climate change on the worldwide distribution of the Australasian fungus Clathrus archeri (Phallales, Basidiomycota)

Keywords: Clathrus archeri, Devil's finger fungus, Biological invasions, Species distribution modelling, Fungal biogeography

come

Introduction

Climate is one of the most important factors explaining species distributions [START_REF] Moura | Disentangling the role of climate topography and vegetation in species richness gradients[END_REF], and global warming is thereby a significant threat to biodiversity, as it may affect animal [START_REF] Radchuk | Adaptive responses of animals to climate change are most likely insufficient[END_REF], plant [START_REF] Suggitt | Widespread effects of climate change on local plant diversity[END_REF], and fungal communities [START_REF] Větrovský | A meta-analysis of global fungal distribution reveals climate-driven patterns[END_REF] as well as entire ecosystems [START_REF] Walther | Community and ecosystem responses to recent climate change[END_REF]. Climate changes strongly influence the abundances and distribution of numerous organisms, and these changes can cause latitudinal or altitudinal range shifts. Global warming can also enhance the spread of some organisms, including exotic ones [START_REF] Hellmann | Five potential consequences of climate change for invasive species[END_REF][START_REF] Dyderski | How much does climate change threaten European forest tree species distributions[END_REF].

As defined by the International Union for Conservation of Nature (IUCN), an invasive species (invader) is an alien organism introduced beyond the area where it occurs naturally. The introduction of invasive species is an agent of change and threatens native biodiversity (http://www.issg.org 1 ). Invasive species can be animals, plants, or fungi that are introduced accidentally or deliberately into a place where they are not normally found. Studies on anthropogenic introductions focus mainly on animals and plants, with too few investigations into microbial species, especially fungi (but see [START_REF] Desprez-Loustau | The fungal dimension of biological invasions[END_REF]. In contrast to macroorganisms, the occurrence of fungi is often difficult to detect without concerted efforts or costly molecular tools.

Apart from fungal pathogens, which cause major ecological impacts on plants and animals, the biogeography and influence of other invasive fungi often remain poorly recognized. According to another definition [START_REF] Blackburn | A proposed unified framework for biological invasions[END_REF], the term 'invasive' can be used for organisms (including fungi) with self-sustaining populations at significant distances from their places of original introduction. The crucial difference between these two definitions is the ability of the particular species to threaten biological diversity in the invaded areas.

The few existing biogeographic studies of non-pathogenic fungi have focused mainly on mycorrhizal fungi [START_REF] Hao | Using species distribution models for fungi[END_REF]. Suilloid fungi, which were cointroduced with their pine tree mycorrhizal partners, are among the well-described examples showing the biogeography and probable pathways of introduction [START_REF] Vellinga | Global patterns of ectomycorrhizal introductions[END_REF][START_REF] Pietras | Current and potential distribution of the ectomycorrhizal fungus Suillus lakei ((Murrill) A.H. Sm. & Thiers) in its invasion range[END_REF][START_REF] Pietras | First record of North American fungus Rhizopogon pseudoroseolus in Australia and prediction of its occurrence based on climatic niche and symbiotic partner preferences[END_REF][START_REF] Pietras | Predicted potential occurrence of the North American false truffle Rhizopogon salebrosus (A.H. Sm.) in Europe[END_REF][START_REF] Policelli | Suilloid fungi as global drivers of pine invasions[END_REF]. Descriptions of the introduction of wood-inhabiting Favolaschia calocera to Italy [START_REF] Vizzini | Alien fungal species distribution: the study case of Favolaschia calocera[END_REF], saprobic Amanita inopinata [START_REF] Fraiture | Amanita inopinata its ecology and expansion in Europe[END_REF] to Belgium, and ectomycorrhizal Amanita phalloides [START_REF] Pringle | Last chance to know? Using literature to explore the biogeography and invasion biology of the death cap mushroom Amanita phalloides (Vaill. Ex fr. : Fr.)[END_REF][START_REF] Pringle | The ectomycorrhizal fungus Amanita phalloides was introduced and is expanding its range on the west coast of North America[END_REF] to the USA are among the few examples of the well-described introduction pathways of fungi.

Clathrus archeri (Berk.) Dring is a saprobic fungus native to Australia, but it is also frequently recorded in New Zealand. In its native range, the fungus has been detected mainly in mountainous regions and across seasides, where it forms sporocarps on the wood debris of eucalypt and nothofagus trees. Currently, it is also frequently recorded in gardens, parks, and other urban areas as a result of the transport and use of wood chips and mulch for gardening. Since the beginning of the twentieth century, it has been gradually spreading throughout the world. Currently, sporocarps of C. archeri are found in North and South America, but the highest concentration of introduced records are from Europe, where it grows mainly in broadleaf and mixed forests and beyond forest areas in gardens, orchards, old cemeteries and meadows [START_REF] Stengl-Rejthar | Expansion of the fungus Clathrus archeri (Berk.) Dring. (Gasteromycetes) in Europe and Poland. -Zeszyty Naukowe Uniwersytetu Jagiellońskiego 752[END_REF]. According to [START_REF] Bãžrsan | Distribution and ecology of Clathrus archeri in Romania[END_REF], C. archeri has no special habitat preferences and occurs frequently on different soil types, commonly forming sporocarps on wood debris, mulch, and forest litter. However, the fungus avoids calcareous soils [START_REF] Parent | Nouvelles données sur la répartition de Clathrus archeri en particulier dans l'Ouest et le Sud-Ouest de l'Europe[END_REF][START_REF] Bãžrsan | Distribution and ecology of Clathrus archeri in Romania[END_REF]. Reports on the distribution of C. archeri in Europe are mostly embedded as short notes or snippets in regional lists since 1920 [START_REF] Parent | État actuel de l'extension de Clathrus archeri (Berkeley) Dring (syn. : Anthurus archeri (Berk.) Ed. Fischer) en Europe et particulièrement en France et au BÉnÉlux[END_REF]. Our previous research conducted in Poland indicated that even at the scale of a single country, latitudinal shifts in the geographical distribution of C. archeri could be detected in a relatively short time period [START_REF] Pietras | Distribution and molecular characterization of an alien fungus Clathrus archeri in Poland[END_REF]. Outside its native range, the impact of C. archeri on native ecosystems is unknown. Negative effects of the fungus on native ecosystems have not been proven so far regarding either its interactions with native mycobiota or its influence on soil biochemistry.

Therefore, based on the IUCN definition, C. archeri cannot yet be termed "invasive" but is rather considered potentially invasive.

The aim of the present study was to present the climaterelated ecology of C. archeri. We chose this fungus because it is a highly suitable model species for invasion ecology research due to the availability and reliability of associated data (both in its native and introduced ranges: GBIF.org, 26 October 2018a, 2018b, respectively; including disjunct occurrences: GBIF.org, 26 October 2018c) and its set of distinct morphological features, which make it easy to identify under field conditions. Considering all these aspects, C. archeri can be regarded as a model organism in studies describing the biogeography and invasion ecology of nonpathogenic fungi [START_REF] Dickie | Towards management of invasive ectomycorrhizal fungi[END_REF]. Additionally, we evaluated the similarities between the bioclimatic niches occupied by native and introduced populations of C. archeri. The impact of global warming on the future distribution of this species was further assessed using a species distribution modelling approach. Our research focused on estimations of the locations of suitable C. archeri habitats during the last glacial maximum (LGM; approximately 22,000 years ago) and on resolving its status in New Zealand, where the fungus is suspected to be nonnative.

Methods

List of localities

Data on the current distribution of C. archeri were gathered from the literature, herbaria, the Global Biodiversity Information Facility (GBIF2 ) and our own extensive field survey. Out of 2222 records submitted to GBIF (31 December 2019) we select only localities of preserved specimens. This approach give possibility to eliminate observations, that could be easily misidentified and confused with other closely-related taxa such as Aseroe rubra-a fungus that occurs in the same area in Australia and New Zealand. Direct use of all GBIF records could be also biased because of autocorrelation of several closely located (in the same 21.62 km 2 equator) records. Locations for which the exact geographic coordinates were available were used exclusively in further analyses. From a total of 150 locations of C. archeri assembled during this study, the multiple species records from each grid cell were removed. A final database included 146 occurrence data (Appendix).

Species distribution modelling

The species distribution modelling was conducted using the maximum entropy method implemented in MaxEnt version 3.3.2 [START_REF] Phillips | A maximum entropy approach to species distribution modeling[END_REF][START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF] based on the species presence-only observations. We use MaxEnt application because its prediction accuracy is stable and reliable, even with limited data and small sample sizes [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF][START_REF] Guo | Prediction of the potential geographic distribution of the ectomycorrhizal mushroom Tricholoma matsutake under multiple climate change scenarios[END_REF]. The dataset of 146 records, which is more than the minimum number of records (10 observed localities) required to obtain reliable predictions in the MaxEnt application [START_REF] Pearson | Predicting species distributions from small numbers of occurrence records: a test case using cryptic geckos in Madagascar[END_REF][START_REF] Wisz | Effects of sample size on the performance of species distribution models[END_REF], was used as input data. Because niche shifts were observed in some invasive species (e.g., [START_REF] González-Moreno | Beyond climate: disturbance niche shifts in invasive species[END_REF], we decided to conduct separated assessment for populations of C. archeri in its native range. For that reason two different approaches were used to conduct analysis. In the first, locations from Australia and New Zealand (58 records) were exclusively used. In the second, all available records of the studied species (58 from Australia and New Zealand, 88 from Europe and South Africa) were included in the modelling. From 19 climatic variables ("bioclims", Table 1) in 2.5 arc minutes (± 21.62 km 2 at the equator) as developed by [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF] and provided by WorldClim (version 1.4 release 3, www.worldclim.org3 ), seven variables were removed due to their significant correlation (above 0.9) as evaluated by the Pearsons' correlation coefficient calculated using ENMTools v1.3. The following variables were excluded from the dataset: bio6, bio7, bio9, bio10, bio11, bio16, and bio17 (Table 1).

In all analyses, the maximum number of iterations was set to 10,000 and the convergence threshold to 0.00001. The "random seed" option, which provided random test partition and background subset for each run, was applied. The run was performed as a bootstrap with 1000 replicates, and the output was set to logistic. All operations on GIS data were carried out on ArcGis 10.6 (ESRI) and QGIS 2.18.20.

Niche similarity

The differences between the niches occupied by the native and introduced populations of C. archeri (as calculated in ENM analyses) were evaluated using the niche identity test indexes: Schoener's D and I statistic as available in ENMTools v1.3 [START_REF] Schoener | The anolis lizards of bimini: resource partitioning in a complex fauna[END_REF][START_REF] Warren | Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution[END_REF][START_REF] Warren | ENMTools: a toolbox for comparative studies of environmental niche models[END_REF]. To visualize the difference between the climatic preferences of both groups of populations, the predicted niche occupancy profiles (PNOs) were created using the Phyloclim package [START_REF] Heibl | Phyloclim: integrating phylogenetics and climatic niche modeling[END_REF].

History and future range changes

The contemporary species-climate relations were projected for the past, present and three future climate change scenarios using the settings described above. The same bioclimatic variables downloaded from the WorldClim database 3 were used in all analyses. The variables based on the coupled global climate models CGCM2 and CGCM3 were used4 [START_REF] Kim | A coupled climate model simulation of the last glacial maximum part 2: approach to equilibrium[END_REF] to estimate the future changes in the distribution of suitable climatic niches for C. archeri. Three various emission scenarios for 2080 were analyzed: A1b (balance fossil and nonfossil energy sources-CCCMA CGCM3 simulation), A2a (high energy requirements-CCCMA CGCM2 simulation), and B2a (low energy requirements-CCCMACGCM2 simulation, IPCC Special Report Emissions Scenarios. Summary for Policymakers 2000, to see detailed information about differences between analyzed scenarios see paper previously published by [START_REF] Wiatrowska | Current occurrence and potential future climatic niche distribution of the invasive shrub Spiraea tomentosa L. in its native and nonnative ranges[END_REF]. This approach was used in numerous recent studies focused on climate change impact on the distribution of various organisms (e.g., [START_REF] Coops | Estimating the vulnerability of fifteen tree species under changing climate in Northwest North America[END_REF]Ferro et al. 2014, Kukwa and[START_REF] Kukwa | Glacial refugia and the prediction of future habitat coverage of the south American lichen species Ochrolechia austroamericana[END_REF]. The bioclimatic data used to evaluate the possible distribution of the studied species during the last glacial maximum (LGM, about 22,000 years ago) were developed and mapped by Paleoclimate Modelling Intercomparison Project Phase II (PMIP2, CCSM; [START_REF] Braconnot | Results of PMIP2 coupled simulations of the mid-Holocene and last glacial maximum -part 1: experiments and large-scale features[END_REF]).

SDMtoolbox 2.3 for ArcGIS was used to visualize changes in the distribution of suitable climatic niches of studied species caused by global warming [START_REF] Kremen | Aligning conservation priorities across taxa in Madagascar with high-resolution planning tools[END_REF]. To compare the distribution model created for current climatic conditions with future models, all SDMs were converted into binary rasters and mapped using the Goode homolosine projection. The presence threshold value (0.5) was estimated based on the values of grids in which C. archeri occurred in models created using present-time data.

Results

The assessment of the C. archeri distribution was made based on the literature, specimens preserved in herbaria, our data collected during field studies and records gathered from open databases. Altogether, a dataset including 146 records of the fungus was created (Fig. 1, Appendix). To show the worldwide predictions of the distribution of the fungus in the past, present and future, we used occurrence information originating from three continents (Australia, Europe and Africa) and New Zealand.

Model evaluation, niche similarity and limiting factors

The average training area under the ROC (receiver operating characteristic) curve (AUC; [START_REF] Delong | Comparing the areas under two or more correlated receiveroperating characteristic curves: a nonparametric approach[END_REF] for the replicate runs received high scores of 0.990-0.994, indicating that the model performed reliably (Table 2). The AUC has been shown to be prevalence-independent [START_REF] Mcpherson | The effects of species' range sizes on the accuracy of distribution models: ecological phenomenon or statistical artefact?[END_REF], and it is considered to be a highly effective indicator of model performance [START_REF] Shabani | A comparison of absolute performance of different correlative and mechanistic species distribution models in an independent area[END_REF]). Among the bioclimatic variables used in this study, precipitation in the driest month was the most important factor affecting the distribution of C. archeri in the model based on native records only (25.8%, 2) with the optimum volume oscillating near 50 mm per month (Table 2, Supplementary Material S1). Other climate parameters that influenced this model were the mean temperature of the wettest quarter (up to 9 °C) and the precipitation of the warmest quarter (between 100 and 300 mm, Supplementary Material S1). The model generated based on all available records of C. archeri indicated that the precipitation in the driest month (optimal value 45 mm) was of lowest importance (between 14.7% and 16%), followed by the annual mean temperature (assessed as the 9 °C maximum) and the maximum temperature of the warmest month (approaching 22 °C, Supplementary Material S2). The contributions of all variables to the MaxEnt models, as described above, are presented in Table 2.

The climatic niches of the native and introduced populations are rather similar according to the climatic niche identity test results (I = 0.957, D = 0.797). The created PNOs (Supplementary Material S1 and S2) showed relatively small differences in the most suitable values of annual mean temperature (bio1), isothermality (bio3), mean temperature of wettest quarter (bio8), precipitation seasonality (bio15), and precipitation of coldest quarter (bio19) between the native and introduced populations, while significant variations were detected for the temperature seasonality (bio4: 2100-8200 for invasive vs 2800-5800 in native populations). Overall, the ecological tolerances estimated here based on the values of bioclimatic variables in areas suitable for the occurrence of C. archeri seem to be broader in populations occurring outside the native range of the fungus (see Supplementary Material S1 vs S2).

Past and present distributions

The estimated range of C. archeri during the LGM was broader than that in the present time (Fig. 1b, Fig. 2a). During the glacial period, refugia of the studied fungus were common in Southeastern Australia (including in Tasmania, which was then an integral part of Australia) and most of New Zealand. However, the highest concentration of climatic niches suitable for C. archeri in the past was observed on several small islands located near the east coast of Australia, which are now underwater. The current observed potential range of C. archeri is smaller by ca. 18% than it was during the LGM (Table 3). Currently, the distribution of suitable climatic niches of the fungus considerably overlaps with the actual acreage in which the fungus is found in its native range (Fig. 1b). According to the created model, climatic niches of the fungus are distributed in southern Australia and both New Zealand islands, with the highest concentrations observed in the Sydney Basin, southern New South Wales, and the East Coastal Plain in Australia, as well as in Tasmania and the North Island of New Zealand. Interestingly, the analysis conducted based on both native and introduced records indicated some potentially available climatic niches in Madagascar and coastal regions of South Africa. Our model predicted that the climatic conditions of the south coasts of Brazil, Uruguay, and Argentina, as well as some disjunctive locations at the foot of the Andes in Chile, are favorable for C. archeri (Fig. 1a).

The model also showed the existence of numerous climatic niches suitable for the fungus in Europe, where the areas potentially available for the fungus were found in the western, central, and southern portions of the continent. Highly suitable regions that are potentially habitable for C. archeri extend from Northern Spain throughout Western and Central Europe to the eastern Polish border and the Carpathian Mountains in the east. Additional suitable habitats were recorded in the Southern Scandinavian Peninsula and the British Isles as well as from the mountainous regions of the Iberian Peninsula, the Dinar Mountains and the Balkans to the southeastern shores of the Black Sea (Fig. 1A).

Future changes in the coverage of suitable habitats

In our model, over 3.5 million km 2 worldwide are currently climatically suitable for C. archeri (Table 2). The total area of the potential occurrence of the fungus is generally expected to decline due to climate change. On a global scale, only one model (A1b scenario; Fig. 3a) projected a small increase (of 6.2%) in suitable climatic niches for C. archeri by 2080 (Table 4). Under both scenarios A2a and B2a, the area of preferred habitats is expected to decline, and the highest contraction was anticipated in the A2a scenario (Fig. 3b). This future climatic niche loss was also predicted in the native range of C. archeri. All generated models showed significant reductions in suitable habitats in Southeastern Australia, with the highest contraction expected under the A1b scenario (Fig. 2). In the pessimistic projection (scenario A1b), the range of C. archeri is expected to be narrowed to the East Coastal Plain. Highly suitable habitats for the fungus will still be available in Tasmania and New Zealand, and these regions are predicted to be future refugia for the studied species. Notably, an expansion of suitable climatic niches in the native range of C. archeri was predicted in the B2a (11.1%) climate change scenario (Fig. 2b, Table 3). According to our analyses, the total loss of suitable habitats for C. archeri outside its native range will be observed in Madagascar, South Africa, and eastern South America. In contrast, the future potential range of the fungus in Chile is projected to grow (Fig. 3). All created models revealed the existence of highly suitable areas across the western coast of South America and throughout Patagonia to Cape Horn. Additionally, according to the A1b model (Fig. 3a), preferred habitats will become available in North America (in Newfoundland and Nova Scotia, as well in the northwest region of the continent, from the seaside of British Colombia to the southern coast of Alaska). However, the greatest overlap between the current and future suitable climatic niche aggregates will be observed in Europe (Fig. 3). All future projections showed that the southern part of the European range of C. archeri will become unsuitable for the fungus. In general, losses of available habitats are expected in the Iberian, Apennine, and Anatolian Peninsulas, as well as in Central Europe and the Balkans. The different scenarios predict different degrees of habitat coverage. The most significant reduction in the C. archeri range is expected under the A2a scenario (Table 4, Fig. 3b). In terms of total available habitat, the studied species is expected to fare better under the A1b and B2a scenarios than under the other scenarios. The suitable areas estimated by these scenarios are projected to be substantially larger than those predicted by the A2a scenario in the southern range of the fungus, where the model outcomes contradict each other to the greatest extent. On a continental scale, the area of suitable C. archeri habitat is generally expected to increase over time as habitat is gained in new regions, especially in northeastern Europe (in the eastern part of the continent and the central and southern parts of the Scandinavian Peninsula), as well in the entire British Isles (Fig. 3). However, our model predicts a notable decline in the global suitable area for the fungus by 2080 under the A2a scenario (Table 4, Fig. 3b).

Discussion

Species distribution models are rather infrequently used in research on nonpathogenic fungi. Thus far, analyses of the future distributions of these organisms have been implemented in biogeographical studies on Boletales in North America [START_REF] Pietras | Current and potential distribution of the ectomycorrhizal fungus Suillus lakei ((Murrill) A.H. Sm. & Thiers) in its invasion range[END_REF][START_REF] Banasiak | Aureoboletus projectellus (Fungi Boletales)an American bolete rapidly spreading in Europe as a new model species for studying expansion of macrofungi[END_REF][START_REF] Pietras | First record of North American fungus Rhizopogon pseudoroseolus in Australia and prediction of its occurrence based on climatic niche and symbiotic partner preferences[END_REF][START_REF] Pietras | Predicted potential occurrence of the North American false truffle Rhizopogon salebrosus (A.H. Sm.) in Europe[END_REF] and on rare wood-inhabiting saprotrophs restricted to European boreal forests [START_REF] Kujawska | Distribution and ecological traits of a rare and threatened fungus Hericium flagellum in Poland with the prediction of its potential occurrence in Europe[END_REF]. Therefore, studies of C. archeri and its introduction outside of its native range contribute to the discussion on the occurrence of nonpathogenic fungi outside their natural habitats. In this study, we attempted to answer three basic questions regarding the historical occurrence, current biogeography and future distribution changes of C. archeri.

Where did C. archeri originate?

In contrast to plants, fungal micro-and macrofossils are recorded extremely rarely in paleoecological studies. The earliest uncontested fungal fossils are multicellular organic-walled microfossils dating to 1010-890 million years ago, discovered in Arctic Canada [START_REF] Loron | Early fungi from the Proterozoic era in Arctic Canada[END_REF]. In general, the lack of paleomycological data confirms that fungi rarely fossilize due to the structures of sporocarps and delicate spores [START_REF] Ascaso | Fossil protists and fungi in amber and rock substrates[END_REF][START_REF] Beimforde | Estimating the Phanerozoic history of the Ascomycota lineages: combining fossil and molecular data[END_REF]. The environmental DNA of fungi dating back to the LGM can sometimes be detected in glacier sheets using molecular tools [START_REF] Gavin | Climate refugia: joint inference from fossil records species distribution models and phylogeography[END_REF]. Species distribution modelling based on current species records is a useful tool for reconstructing the history of fungi when fossil material is absent. In Australia and Tasmania, species distribution models have been successfully used to reconstruct the occurrence patterns of past refugia, where suitable habitats for endemic plants occurred [START_REF] Mokany | Past present and future refugia for Tasmania's palaeoendemic flora[END_REF][START_REF] Das | Identifying climate refugia for 30 Australian rainforest plant species from the last glacial maximum to 2070[END_REF]. Our study showed for the first time the potential distribution of a saprotrophic fungus, C. archeri, 22,000 years ago during the LGM. At that time, the climate throughout eastern Australia was cooler and drier than it is at present [START_REF] Das | Identifying climate refugia for 30 Australian rainforest plant species from the last glacial maximum to 2070[END_REF]. Such conditions were favorable for the broader occurrence of this fungus during the LGM than that at present. A favorable past environment is also confirmed in our results, showing the precipitation of the driest month and the mean temperature of the wettest quarter, together contributing more than 50%, as main predictors of C. archeri occurrence (Table 2). When the climate becomes warmer and wetter, the range of the fungus decreases to 18.7 % of Australian hot-spot areas (Southeastern Queensland and the North Coast of New South Wales, [START_REF] Das | Identifying climate refugia for 30 Australian rainforest plant species from the last glacial maximum to 2070[END_REF] as well to the southern states of New South Wales and Victoria (Fig. 2b). In contrast, our model indicated that the highly suitable climatic niches of C. archeri are located in Tasmania and New Zealand in both the past and present. The wide occurrence of suitable climatic niches during the LGM in New Zealand suggests that C. archeri may have occurred in this area during that period (Fig. 2A). Therefore, we hypothesize C. archeri is a likely natural component of the New Zealand mycobiota; however, the natural divergence of the Australian and New Zealand populations could have occurred before the glacial period. This hypothesis needs to be investigated using genetic tools to highlight the differences between the populations and species histories. Population genetic analysis would also reveal demographic events, such as bottlenecks and founder effects, within populations of alien/invasive fungi, which may underlie the invasion success of these fungi [START_REF] Gladieux | The population biology of fungal invasions[END_REF][START_REF] Séne | A pantropically introduced tree is followed by specific ectomycorrhizal symbionts due to pseudo-vertical transmission[END_REF]. This approach may also reveal the possibility of the occurrence of cryptic relatives of the fungus, especially in distinct areas where fungal sporocarps can differ in size and morphological features.

How has the fungus found new homes?

C. archeri was introduced incidentally to Europe in approximately 1920 [START_REF] Parent | État actuel de l'extension de Clathrus archeri (Berkeley) Dring (syn. : Anthurus archeri (Berk.) Ed. Fischer) en Europe et particulièrement en France et au BÉnÉlux[END_REF], probably with Australian wool [START_REF] Parent | Nouvelles données sur la répartition de Clathrus archeri en particulier dans l'Ouest et le Sud-Ouest de l'Europe[END_REF] or with Australian soldiers during the First World War (Desprez-Loustau and Rizzo 2011).

The probable pathways of C. archeri introductions to other parts of the world are unknown. As revealed in the current study, the climatic conditions preferred by the fungus in Europe are comparable to those observed in its native range.

The similarity of available habitats allows quick expansion to new areas, as no new adaptation is required, at least for the studied traits of the niches. The broader ecological tolerance of the studied species in its introduced range compared to that in its native range revealed in the PNO analysis also suggests that either the species does not use its fundamental climatic niche completely or that it was able to shift its preferences outside the conditions of its native geographical region. Perhaps the enlargement of its niche was allowed by losses of mortality sources, i.e., the absence of pathogens and parasites of the species [START_REF] Liu | Testing the enemy release hypothesis: a review and meta-analysis[END_REF]. Considering the analyzed variables, populations outside the native range can occur in areas with isothermality (bio3) values of ca. 22-38% and 45-55% (native populations are localized within areas with values of ca. 43-53% only), temperature seasonality (bio4) values of ca. 2100-8200 (vs 2800-5800), and precipitation seasonality (bio15) values up to 80 mm (vs up to 40 mm).

The current potential range of C. archeri in Europe is much broader than its actual occurrence as revealed by the occurrence data (Table 5), further confirming our previous studies showing that the realized niches of nonnative fungi are smaller than the predicted niches [START_REF] Pietras | Current and potential distribution of the ectomycorrhizal fungus Suillus lakei ((Murrill) A.H. Sm. & Thiers) in its invasion range[END_REF][START_REF] Banasiak | Aureoboletus projectellus (Fungi Boletales)an American bolete rapidly spreading in Europe as a new model species for studying expansion of macrofungi[END_REF][START_REF] Pietras | First record of North American fungus Rhizopogon pseudoroseolus in Australia and prediction of its occurrence based on climatic niche and symbiotic partner preferences[END_REF][START_REF] Pietras | Predicted potential occurrence of the North American false truffle Rhizopogon salebrosus (A.H. Sm.) in Europe[END_REF]. Our model indicated that numerous suitable habitats are located in Europe, e.g., in the British Isles, Apennine Peninsula, Carpathians, Balkans, and Eastern Europe, where no GBIF records of the fungus are currently available (GBIF.org, 26 October 2018b). Clathrus archeri, similar to other fungi, is a cryptic organism most of the year and is found mostly during field surveys. Ephemeral fruiting bodies make it challenging to detect specimens in field surveys. Moreover, the production of sporocarps strongly depends on climatic conditions, mainly on precipitation during the driest month [START_REF] Pietras | Distribution and molecular characterization of an alien fungus Clathrus archeri in Poland[END_REF]. In contrast to regional scales, where soil conditions become more important, the climate can be considered the most important factor influencing the species distribution at the continental scale [START_REF] Pearson | Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful?[END_REF]. Precipitation values during the driest or coldest quarters were the most decisive factors influencing the occurrence of suilloid fungi outside their native ranges [START_REF] Pietras | Current and potential distribution of the ectomycorrhizal fungus Suillus lakei ((Murrill) A.H. Sm. & Thiers) in its invasion range[END_REF][START_REF] Pietras | First record of North American fungus Rhizopogon pseudoroseolus in Australia and prediction of its occurrence based on climatic niche and symbiotic partner preferences[END_REF][START_REF] Pietras | Predicted potential occurrence of the North American false truffle Rhizopogon salebrosus (A.H. Sm.) in Europe[END_REF]. Changes in precipitation during the vegetation season are regarded as crucial drivers affecting the sporocarp production of edible fungi [START_REF] Taye | Meteorological conditions and site characteristics driving edible mushroom production in Pinus pinaster forests of Central Spain[END_REF][START_REF] Salerni | Linking climate variables with Tuber borchii sporocarps production[END_REF]. Precipitation in the coldest month determines the seasonality of sporocarp production. Therefore, this study has shown that the assessed range of C. archeri in Europe is underestimated compared to that determined by specimens observed to date. This assumption can be supported by our latest findings in northern Poland, where, using next-generation sequencing analysis of soil DNA, the presence of C. archeri was recorded in sites where no sporocarps had ever been found (Pietras unpublished). Based on the results obtained in this study and our best knowledge about the biology and ecology of C. archeri, we conclude that the fungus could be much more prevalent in Europe than previous evidence would indicate.

In the present study, we used an extended dataset of C. archeri records (Appendix), including occurrence data from South Africa and Madagascar (GBIF.org, 26 October 2018c). In these regions, the created models of C. archeri occurrence overlapped with the observed range of the fungus. African records were classified as nonnative because they constitute a distant and separated group from the main range of C. archeri in Australia and New Zealand. Moreover, sporocarps recorded in Africa indicate different morphological features, mostly consisting of sporocarps colored from red to yellow. There is also no molecular evidence clarifying the phylogenetic relationship between C. archeri in its native range and those found in South Africa. Therefore, further molecular studies need to be performed to clarify the taxonomy and systematics of the fungus. Our model also predicted climatic conditions favorable for the fungus in regions where it has never been found (Brazil, Uruguay, Argentina, and Chile). However, sporocarps of C. archeri were recently found in other regions of South America (Pinzón-Osorio and Pinzón-Osorio 2020). Citizen science observations revealed the presence of C. archeri in Bogota Province in Colombia5 ; however, these isolated records could not be included for modelling purposes in this study [START_REF] Pearson | Predicting species distributions from small numbers of occurrence records: a test case using cryptic geckos in Madagascar[END_REF]. Similarly, single isolated sites were noted in California in the USA, with the first record of C. archeri occurring in the early 1980s. (voucher no. NY01938977). Apart from the observations above, no other specimens have been preserved in North American fungaria.

Quo vadis?

It is notable that within its native range, the coverage of suitable climatic niches of C. archeri has decreased since the LGM, but anthropogenic climate changes have accelerated habitat loss. This process of habitat loss is especially intensive in Australia, where the highest rate of range contraction is expected to occur (Fig. 2c-e). Based on the most pessimistic climate change scenario (A1b), the acreage of optimal habitats of the fungus will decrease to the limited area of the Australian Alps and Coastal Plains in southeast Australia (Fig. 2c). Therefore, C. archeri would become a threatened species in all of Australia, even though the IUCN attributes the vulnerable status for saprotrophic fungi whose occupied area tends to decline by more than 30% over 15-30 years. In contrast, no habitat loss or even range extension is observed in Tasmania and New Zealand, where suitable habitats will still be available in the areas where the fungus potentially occurs; the distribution of the fungus will even enlarge into new regions. Generally, islands are among the most vulnerable to habitats to climate changes regarding shifts in rainfall, temperature, and sea level, including the expected negative influences of these changes on biodiversity [START_REF] Veron | Vulnerability to climate change of islands worldwide and its impact on the tree of life[END_REF], Office of The Prime Minister's Science Advisory Committee 2013). However, based on the obtained predictions, no negative impacts due to climate change on the distribution of C. archeri will be observed in either Tasmania or New Zealand in the future, where a wide range of areas will constitute highly valuable refugia. Different responses to climate change are predicted to occur in the native and introduced ranges of C. archeri. In Australia and New Zealand, habitat losses were predicted in the A1b and A2a climate change scenarios, the declines in suitable climatic niche coverages in other parts of the world are expected to occur in the A2a and B2a scenarios. In the most extreme climate simulations, the species will lose 7.0% (native range)-28.5% (introduced range) of its preferred habitats. The estimated increase in suitable climatic niche coverage will not be significant within the invasive range (6.2%), but it is expected to be twice as high within Australia and New Zealand (11.1%).

Generally, the most common responses of organisms to climate change in Europe include latitudinal or altitudinal shifts in the geographical distributions of the organisms [START_REF] Lenoir | Climate-related range shiftsa global multidimensional synthesis and new research directions[END_REF][START_REF] Gange | Trait-dependent distributional shifts in fruiting of common British fungi[END_REF]. For C. archeri, the effect of ongoing climate change is predicted to involve migration from Southwest and Central Europe to areas farther northeast. Additional suitable climatic niches will appear at higher elevations of the Alps and Carpathians (Table 3). The spread of non-pathogenic fungi seems to be unstoppable. Therefore, potential invasion of C. archeri should be considered as a possible threat to native mycobiota. Moreover, in the case of the studied fungus, any management strategies proposed by [START_REF] Dickie | Towards management of invasive ectomycorrhizal fungi[END_REF], such as blocking animal vectors, cannot be used to slow down the expansion of the fungus. The ratio between the areas of contraction and expansion depends on the analyzed climate change scenario. The highest intersection between the current potential and future suitable climatic niche coverages is observed in the B2a scenario, where areas in Central Europe that are suitable for C. archeri will be maintained. In contrast, the most pessimistic climate change scenario (A2a, Table 3B) reveals significant habitat loss in the central part of Europe, making the occurrence of the fungus limited mainly to Eastern Europe and several disjointed areas in the central region of the continent.

Conclusion

Fungi are among the most poorly studied groups of organisms regarding their biogeography and invasion ecology. Based on a worldwide dataset originating from open databases, fungaria, and field observations, we have prepared a comprehensive description of the past, current and future distributions of C. archeri. Our study indicated that (1) the area of habitats suitable for the fungus has been decreasing since the LGM and (2) anthropogenic climate change may accelerate the process of habitat loss. Different responses to climate change are predicted to occur in different parts of the world. In the most pessimistic scenario, the highest risks of regional fungal extinction were predicted in Australia and South Africa. As a native species, C. archeri should be regarded as a threatened fungus in Australia. In other regions within its native range, climate change will not be entirely detrimental for C. archeri, and these areas can be regarded as highly valuable future refugia. The expansion of the fungus in Europe is a result of the similarity between the climatic conditions preferred by the fungus in its native range and those in Europe. This similarity allowed a quick expansion of the fungus to new areas, as no adaptation was required. This study also revealed the presence of suitable climatic niches during the LGM in both Australia and New Zealand, supporting the possibility that C. archeri is a native species in both countries.

Fig. 1

 1 Fig. 1 Present potential distribution of C. archeri: a Global range, b native range. Dots represent localities used in the modelling

Fig. 2

 2 Fig. 2 Potential distribution of C. archeri in its native range during LGM (a) and the changes in area suitable for the occurrence of C. between the present and LGM (b), future distribution according to A1b (CCCMA CGCM3 simulation) (c), A2a (CCCMA CGCM2 simulation)

Fig. 3

 3 Fig. 3 Changes in the area suitable for the occurrence of C. archeri between presence and future according to A1b (CCCMA CGCM3 simulation) (a), A2a (CCCMA CGCM2 simulation) (b) and B2a

Table 1

 1 Codes of climatic variables developed by[START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF], variables used in this study in bold

	Code Description
	bio1 Annual mean temperature
	bio2 Mean diurnal range=mean of monthly (max temp-min temp)
	bio3 Isothermality (bio2/bio7) * 100
	bio4 Temperature seasonality (standard deviation * 100)
	bio5 Max temperature of warmest month
	bio6 Min temperature of coldest month
	bio7 Temperature annual range (bio5 -bio6)
	bio8 Mean temperature of wettest quarter
	bio9 Mean temperature of driest quarter
	bio10 Mean temperature of warmest quarter
	bio11 Mean temperature of coldest quarter
	bio12 Annual precipitation
	bio13 Precipitation of wettest month
	bio14 Precipitation of driest month
	bio15 Precipitation seasonality (coefficient of variation)
	bio16 Precipitation of wettest quarter
	bio17 Precipitation of driest quarter
	bio18 Precipitation of warmest quarter
	bio19 Precipitation of coldest quarter

Table 2

 2 The AUC indexes, coverage in native and outside native range and the estimates of relative contributions of the most important environmental variables to the MaxEnt models. Standard deviation values are given in brackets

	Model	LGM	Present	A1b	A2a	B2a
	Dataset	Native	Native and nonnative
	AUC (SD)	0.995 (0.002)	0.994 (0.002)	0.993 (0.002)	0.990 (0.002)	0.994 (0.002)
	Coverage [km 2 ]	574,229.3	3,587,973.29	3,824,713.22	2,793,158.65	3,361,766.02
	Coverage AUS and NZ [km 2 ]	574,229.37	451,213.60	421,856.36	443,098.93	507,378.98
	Coverage Beyond native range [km 2 ]	-	3,136,759.69	3,402,856.87	2,350,059.72	2,854,397.04
	Var_1	bio14 (25.8)	bio14 (16.0)	bio14 (14.7)	bio14 (15.8)	bio14 (15.4)
	Var_2	bio8 (25.4)	bio1 (13.1)	bio1 (14.2)	bio1 (13.3)	bio1 (15.3)
	Var_3	bio18 (8.7)	bio5 (12.3)	bio5 (12.8)	bio5 (12.8)	bio5 (13.8)

Table 3

 3 Changes in area suitable for the occurrence of C. archeri within its native range [km 2 ]

	Range expansion	No occupancy	No change	Range contraction	Range change
	A	B	C	D	(C+D)/(C+A)
	LGM vs present [km 2 ] (Fig. 2b)
	160,338.98	928,5491.98	413,890.39	69,675.66	Contraction 18.7%
	Present vs A1b [km 2 ] (Fig. 2c)
	93,403.62	912,1499.27	328,452.74	122,760.86	Contraction 7.0%
	Present vs A2a [km 2 ] (Fig. 2d)
	97,269.76	9,117,633.12	345,829.17	10,5384.43	Contraction 1.8%
	Present vs B2a [km 2 ] (Fig. 2e)
	125,182.51	9,089,720.38	382,196.46	69,017.14	Expansion 11.1%

Table 4

 4 Changes in area suitable for the global occurrence of C. archeri [km 2 ]

	Range expansion	No occupancy	No change	Range contraction	Range change
	A	B	C	D	(C+D)/(C+A)
	Present vs A1b [km 2 ] (Fig. 3a)
	2,098,826.09	12,9546,234.07	1,725,887.13	1,862,086.15	expansion 6.2%
	Present vs A2a [km 2 ] (Fig. 3b)
	1,461,942.69	130,183,117.47	1,331,215.96	2,256,757.33	contraction 28.5%
	Present vs B2a [km 2 ] (Fig. 3c)
	1,353,516.98	130,291,543.19	2,008,249.04	1,579,724.25	contraction 6.7%

Table 5

 5 The occurrence of Clathrus archeri in Europe according to all available data submitted to Global Biodiversity Information Facility (GBIF, 31 December 2019) and localities used in this study (records of preserved specimens)

	Country	Number of all records	Number of records
		submitted to GBIF	used in this study
	Poland	27	23
	Germany	370	13
	Spain	138	13
	France	143	2
	Netherlands	831	2
	Slovenia	203	1
	Austria	50	1
	Sweden	32	1
	Italy	13	1
	Luxembourg	11	1
	Czech Republic	9	1
	Switzerland	8	1
	Croatia	4	1
	Portugal	2	1
	Bosna and Herzegovina 1	1

www.gbif.org

www.worldclim.org

http://ccafs-climate.org

Mycol Progress (2021) 20:299-311

https://mushroomobserver.org/266860

Acknowledgments We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the National Center for Atmospheric Research for producing and making available their model output. For CMIP, the U.S. Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and leads the development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. We thank Dr. Margarita Dueñas (Curator of Herbarium of the MA-Fungi), Dr. Ulrike Damm (Curator of Herbarium of Senckenberg Museum of Natural History, Görlitz), Dr. D. Christine Cargill (Curator of Cryptogam Herbarium Australian National Herbarium, CANB), Pina Milne (Curator of Royal Botanic Garden Victoria Herbarium), Dr. Matthew Pace (Assistant Curator of Herbarium of The New York Botanical Garden). We are grateful to Patric Alber, Bronisław Burza, Anne-Marie Grenier, M. Mayer, Maurice Pelissier, Laura Schneider-Maunoury, Zuzanna Michalska, Sebastian Piskorski and Rozalia Ziółkowska for help in specimen gathering. Special thanks to Prof Ian Dickie, who assisted us in collecting material during our trip to New Zealand.

Funding This study was supported by the Institute of Dendrology Polish Academy of Sciences within statutory activity and the Research Fund for Young Scientists in 2013, by the Ministry of Education, Youth and Sports of CR within the CzeCOS program [grant number LM2018123] and by the Polish National Science Center under the funding of an internship after obtaining a doctoral degree [decision no. DEC-2015/16/S/NZ9/ 00370].

Data availability

The authors declare that all data and materials as well as software application support their published claims and comply with field standards.

Code availability Not applicable

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1007/s11557-021-01669-w.

Authors' contribution All authors contributed to the study conception and design. Material preparation, data collection, and analysis were performed by Marcin Pietras, Marta Kolanowska, and Marc-André Selosse. The first draft of the manuscript was written by Marcin Pietras and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Compliance with ethical standards

Competing interests The authors declare that they have no competing interests.

Ethics approval Not applicable.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


C Métais

C Petuya

S Espert

J Ollivier

L Martin-Gondre

A Desmedt
email: arnaud.desmedt@u-bordeaux.fr

Nitrogen Hydrate Cage Occupancy and Bulk Modulus Inferred from DFT-derived Cell Parameters

Gas clathrate hydrates solid materials, ubiquitous in the nature as found either on the ocean floor, permafrost on earth or in extraterrestrial planets and comets, are also technologically relevant, e.g., for energy storage or carbon dioxide sequestration. Nitrogen hydrate, in particular, is of great interest as a promoter of the kinetics of methane replacement reaction by carbon dioxide in natural gas hydrates. This hydrate may also appear in the chemistry of planets wherever the nitrogen constitutes the majority of the atmosphere. A fine understanding of the stability of this hydrate under various thermodynamic conditions is thus of utmost importance to assess its role in the many fields where it occurs. In the present work we have investigated the structural properties of the nitrogen hydrate by means of DFT calculations. We show that the lattice parameters strongly depend onto the cage occupancy and that sI structure has higher bulk elasticity than sII structure.

An energy analysis reveals the key role played by the cage occupancy onto the type of hydrate structure formed, that could be used for experimentally estimating the cage occupancy through the lattice parameter measurement.

I. Introduction

Gas clathrate hydrates are natural porous solid materials. Since few decades, they became of particular interest for their applications in various research fields, such as energy, technology and environment. [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF][START_REF]Gas Hydrates 1: Fundamentals, Characterization and Modeling[END_REF][START_REF]Gas Hydrates 2: Geosciences and Applications[END_REF] As examples, gas hydrates are involved in flow assurance, in gas separation technologies or for desalination and they are considered as a substantial energy reservoir. [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF][START_REF]Gas Hydrates 1: Fundamentals, Characterization and Modeling[END_REF][START_REF]Gas Hydrates 2: Geosciences and Applications[END_REF] These systems are also studied from a fundamental point of view in the areas of physical-chemistry and astrophysics. [START_REF]Gas Hydrates 1: Fundamentals, Characterization and Modeling[END_REF][START_REF]Gas Hydrates 2: Geosciences and Applications[END_REF] Nitrogen hydrate has been widely studied due to its particular relevance. For example, a new technology has been developed for the replacement of methane by carbon dioxide in natural gas hydrates, using N2 as a promoter of the replacement kinetics. [START_REF] Qin | Calibration of Raman Quantification Factors of Guest Molecules in Gas Hydrates and Their Application to Gas Exchange Processes Involving N2[END_REF] The investigation of nitrogen hydrates is also relevant in the fields of atmospheric chemistry and in astrophysics, since nitrogen is one of the most abundant constituents of the atmosphere of various planets in our solar system. Researches conducted in that fields rely on clathrate hydrates to explain and reconstruct the evolution of the atmospheric composition, either on Earth [START_REF] Pauer | Raman Spectroscopic Study on the Nitrogen Ratio in Natural Ice Clathrates in the GRIP Ice Core[END_REF][START_REF] Pauer | Raman Spectroscopic Study on the Spatial Distribution of Nitrogen and Oxygen in Natural Ice Clathrates and their Decomposition to Air Bubbles[END_REF][START_REF] Champagnon | Nitrogen and Oxygen Guest Molecules in Clathrate Hydrates: Different Sites Revealed by Raman Spectroscopy[END_REF][START_REF] Kuhs | Cage Occupancy and Compressibility of Deuterated N2-Clathrate Hydrate by Neutron Diffraction[END_REF][START_REF] Chazallon | Micro-Raman Analysis of Synthetic Air Clathrates[END_REF][START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF] and on planets or comets, such as 67P/Churyumov-Gerasimenko [START_REF] Lectez | A ∼32-70 K Formation Temperature Range for The Ice Grains Agglomerated by Comet 67P/Churyumov-Gerasimenko[END_REF] .

Clathrate hydrates are composed of water molecules, arranged in cages forming the host network. Small gas molecules trapped in these cages are referred to as guest molecules. Nitrogen hydrate is known to crystallize into the so-called type I and type II cubic structures (called sI and sII), lattices with unit cell parameter of approximately 12 Å and 17 Å, respectively. The sI and sII unit cells consist of 46 and 136 water molecules, respectively. Both structures are composed of two types of cages, commonly called small cages and large cages (labeled SC and LC, respectively). The sI contains 2 pentagonal dodecahedra (noted 5 [START_REF] Lundgaard | Calculation of Phase Diagram of Gas Hydrates[END_REF] ) SCs and 6 tetradecahedron (5 [START_REF] Lundgaard | Calculation of Phase Diagram of Gas Hydrates[END_REF] 6 2 ) LCs while the sII is made of 16 SCs (5 [START_REF] Lundgaard | Calculation of Phase Diagram of Gas Hydrates[END_REF] ) and 8 hexakaidecahedral (5 [START_REF] Lundgaard | Calculation of Phase Diagram of Gas Hydrates[END_REF] 6 4 ) LCs. The small cages are the same in both structures, but the large cages in sII are slightly bigger than in sI.

Thermodynamic modeling of nitrogen hydrate predicted that sI would be the most stable form at high pressures (>1000 bar) while sII would be more stable at lower pressures (~100 bar) [START_REF] Lundgaard | Calculation of Phase Diagram of Gas Hydrates[END_REF] , as shown in Fig. 1.

Experimentally, it has been proved that, at moderate pressure, nitrogen hydrate initially forms a type I structure and transforms into its thermodynamically stable form, the type II structure, after few days. [START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF] Similar structural metastability has been observed for carbon monoxide clathrate hydrate, which also forms in sI and converts into sII after 17 weeks. [START_REF] Zhu | Encapsulation Kinetics and Dynamics of Carbon Monoxide in Clathrate Hydrate[END_REF][START_REF] Petuya | Guest Partitioning in Carbon Monoxide by Raman Spectroscopy[END_REF] Small molecules, like N2 or CO, preferentially form the type II structure due to their ability to fill the larger SCs proportion in sII when compared to sI. It becomes then of great interest to investigate how the cages occupancy by the guest molecules influence the gas hydrate structural stability. corresponds to the melting equilibrium of water. Black symbols correspond to experimental data. [START_REF] Sugahara | Thermodynamic Stability and Structure of Nitrogen Hydrate Crystal[END_REF][START_REF] Marshall | Hydrates at High Pressures: Part I. Methane-Water, Argon-Water, and Nitrogen-Water Systems[END_REF][START_REF] Van Cleeff | Gas Hydrates of Nitrogen and Oxygen[END_REF] The red square is the dissociation point of the nitrogen hydrate measured experimentally at 1 bar. [START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF] The cage occupancy corresponds to the average number of guest molecules in small or large cages. For instance, empty sI SCs correspond to θSC = 0 (all SC empty) and one guest molecule per sI SC correspond to θSC = 1 (all SCs filled). Experimental studies with various techniques have been conducted to investigate the occupancy qSC and qLC of small and large cages, respectively, in nitrogen clathrate hydrates. [START_REF] Qin | Calibration of Raman Quantification Factors of Guest Molecules in Gas Hydrates and Their Application to Gas Exchange Processes Involving N2[END_REF][START_REF] Kuhs | Cage Occupancy and Compressibility of Deuterated N2-Clathrate Hydrate by Neutron Diffraction[END_REF][START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF][START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF][START_REF] Hansen | Lattice Constants and Expansivity of Gas Hydrates From 10 K up to the Stability Limit[END_REF][START_REF] Petuya | Ageing and Langmuir Behavior of the Cage Occupancy in the Nitrogen Gas Hydrate[END_REF] These studies revealed that N2 molecules were localized in both types of cages, regardless of the structure formed. Notably, neutron diffraction studies showed that a sII nitrogen hydrate formed at 150 bar and 258 K for 20 days had θSC = 0.87 and θLC = 1.2 (meaning that 20% of the sII LCs are occupied by two nitrogen molecules) [START_REF] Hansen | Lattice Constants and Expansivity of Gas Hydrates From 10 K up to the Stability Limit[END_REF] and that a sI nitrogen hydrate formed at 1093 bar and 268 K for 11 days had θSC = 0.98 and θLC = 1.12 4 . Overall, SC occupancy ranges from 82 to 100 % and LC occupancy ranges from 97 to 125 %. The LC double occupancy has almost always been observed for samples studied by neutron diffraction. A Raman scattering study of the nitrogen hydrate recently evidenced that the LCs have the ability to catch and release guest molecules by simply varying T-P conditions. [START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF] This ability was also observed in the case of the carbon monoxide hydrate. [START_REF] Petuya | Guest Partitioning in Carbon Monoxide by Raman Spectroscopy[END_REF] The effect of both the temperature and pressure on cage occupancy and lattice parameter of nitrogen hydrates has been described. In general, for a given temperature, the cage occupancy increases with increasing pressure [START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF][START_REF] Kuhs | Filling-Isotherms in Clathrate Hydrates[END_REF] and, at fixed pressure, the cage occupancy increases when the temperature increases [START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF][START_REF] Sasaki | Microscopic Observation and In-Situ Raman Scattering Studies on High-Pressure Phase Transformations of a Synthetic Nitrogen Hydrate[END_REF] . However, at a given pressure, the lattice parameter increases with temperature [START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF][START_REF] Hansen | Lattice Constants and Expansivity of Gas Hydrates From 10 K up to the Stability Limit[END_REF] , whereas decreases with the pressure (for P>500bar) at constant temperature [START_REF] Kuhs | Cage Occupancy and Compressibility of Deuterated N2-Clathrate Hydrate by Neutron Diffraction[END_REF][START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF] . Studying these effects of pressure on lattice parameters has also led to the determination of the compressibility of the nitrogen hydrate. Experimentally, the structure type I has a lower bulk modulus BT than structure type II, [START_REF] Kuhs | Cage Occupancy and Compressibility of Deuterated N2-Clathrate Hydrate by Neutron Diffraction[END_REF][START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF] meaning that sI may undergo more structural deformations than sII.

Several theoretical studies of the nitrogen hydrate have also been conducted. Molecular dynamics calculations were performed on sII nitrogen hydrate to investigate the possibility of LCs double occupancy. [START_REF] Van Klaveren | Stability of Doubly Occupied N 2 Clathrate Hydrates Investigated by Molecular Dynamics Simulations[END_REF][START_REF] Van Klaveren | Molecular Dynamics Simulation Study of The Properties of Doubly Occupied N 2 Clathrate Hydrates[END_REF] These calculations showed that the hydrate structure remained stable for all tested pressures, temperatures and compositions. The volume and energies were found to be linearly dependent on the cage occupancy, and there was no marked difference in the structure between singly and doubly occupied systems. Grand-canonical Monte-Carlo studies were performed on sI and sII nitrogen clathrate hydrates. [START_REF] Patt | A Grand Canonical Monte-Carlo Study of the N2, CO and Mixed N2-CO Clathrate Hydrates[END_REF][START_REF] Ballenegger | Cage Occupancies in Nitrogen Clathrate Hydrates from Monte Carlo Simulations[END_REF] These studies showed that at low temperature (at 50K), the SCs are filled before the LCs, but at higher temperature (at 150K) the opposite happens, i.e. the LCs are filled before the SCs. Carbon monoxide hydrate has also been simulated and, in the latter case, whatever the structure formed or the temperature, the SCs were always occupied before the LCs, although at 150 K the occupancy of both types of cages becomes simultaneous.

All these results highlight the key role of cage occupancy, and in particular the importance of water-guest interactions on the structural stability of the hydrate. [START_REF] Desmedt | Energy Landscape of Clathrate Hydrates[END_REF] However, as important as they are, these water-guest interactions cannot be considered alone to fully describe the structural and energetic properties of the system. [START_REF] Cox | Benchmarking the Performance of Density Functional Theory and Point Charge Force Fields in Their Description of SI Methane Hydrate Against Diffusion Monte Carlo[END_REF] Indeed, it is important to consider the host-host interactions contained in the H-bonded solid network formed by the water molecules that are the main constituent of a gas hydrate. [START_REF] Chakraborty | Hydrogen-Bond Vibrational and Energetic Dynamical Properties in SI and SII Clathrate Hydrates and in Ice Ih: Molecular Dynamic Insights[END_REF] As an example, the encapsulation of strong acids within the water cages leads to generate H-bond defects (breaking the socalled ice rule) within the water network modifying the H-bonded water dynamics [START_REF] Desmedt | Proton Dynamics in the Perchloric Acid Clathrate Hydrate HClO4-5.5H2O[END_REF][START_REF] Desmedt | Hydronium Dynamics in the Perchloric Acid Clathrate Hydrate[END_REF][START_REF] Bedouret | Proton Diffusion in the Hexafluoro-phosphoric Acid Clathrate Hydrate[END_REF] and physicochemical properties such as hydrate stability [START_REF] Desmedt | Modifying the Flexibility of Water Cages by Co-Including Acidic Species Within Clathrate Hydrates[END_REF] or formation kinetics [START_REF] Nguyen | Promoting the Insertion of Molecular Hydrogen in Tetrahydrofuran Hydrate with The Help of Acidic Additives[END_REF] . The structural stability of the nitrogen hydrate must therefore be investigated considering both guest-host and host-host interactions.

Nowadays, to study molecular interactions in condensed phase systems, it is appropriated to use periodic density functional theory (DFT) calculations. [START_REF] Burke | Perspective on Density Functional Theory[END_REF] Two types of interactions are involved in gas hydrates: the host-host interactions between water molecules, which correspond to hydrogen bonding, and the guest-host interactions between the guest molecule and the host cage, which is mainly due to van der Waals (vdW) dispersion forces. It is important that DFT correctly describes these two kinds of interactions and in particular the vdW interactions that contribute to the structural stabilization of the hydrate with respect to ice. Unfortunately, the most common exchange-correlation functionals in DFT, such as local density approximation (LDA) or generalized gradient approximation (GGA), are not able to describe the long-range electron correlations that are responsible for the dispersion forces. Although many dispersion-based DFT methods have been developed over the recent years 36 , there is not yet a method fully capable to provide an accurate description of clathrate hydrates [START_REF] Cox | Benchmarking the Performance of Density Functional Theory and Point Charge Force Fields in Their Description of SI Methane Hydrate Against Diffusion Monte Carlo[END_REF][START_REF] Gillan | Perspective: How Good is DFT for Water?[END_REF] .

DFT-based clathrate hydrates structural stability studies have been reported recently. [START_REF] Vlasic | Atomistic Modeling of Structure II Gas Hydrate Mechanics: Compressibility and Equations of State[END_REF][START_REF] Petuya | Unraveling the Metastability Of the SI and SII Carbon Monoxide Hydrate with A Combined DFT-Neutron Diffraction Investigation[END_REF] In the first study [START_REF] Vlasic | Atomistic Modeling of Structure II Gas Hydrate Mechanics: Compressibility and Equations of State[END_REF] , DFT has been used in order to determine the equilibrium lattice volume and bulk modulus of sII hydrocarbons hydrates using several functionals and equations of state with a comparison to experimental data. This investigation concluded that the experimental lattice parameters are better reproduced by using the revPBE functional (revised Perdew-Burke-Ernzerhof), and the bulk modulus of sII hydrate was found to increase with increasing atomic volume. In the second study 39 , DFT calculations have been performed to investigate the structural stability and energy properties of the carbon monoxide hydrate. The simulations were done on complete sI and sII unit cells, with varying cage occupancy, and testing two functionals (PBE [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] and vdW-DF 41 ) to study the impact of the dispersion forces on the gas hydrates properties. Optimized lattice parameters were compared to neutron diffraction and a good agreement was found between theoretical and experimental results. The hydrate structural stability was also investigated through a careful potential energy analysis and results showed that sII is the thermodynamically stable CO hydrate structure. Furthermore, increasing the CO content in the hydrate, i.e. doubly occupying the LCs, was found to destabilize the sI and to have a stabilizing effect on the sII.

The aim of the present study is to investigate the structural and energetic properties of nitrogen clathrate hydrate. Both sI and sII unit cells are considered, and DFT-derived structural properties are compared to experimental data. An analysis of the DFT potential energies is presented and aims at providing a better understanding of the hydrate structural stability. Two different exchange-correlation functionals are used, namely the semilocal functional Perdew-Burke-Ernzerhof (denoted PBE) and the nonlocal functional vdW-DF. The PBE functional does not include dispersion effects but still provide a correct description of the Hbonding network. Therefore, the comparison with the vdW-DF functional, which accounts for dispersion interactions, allows to investigate the effect of the dispersion forces on the calculations. Since the hydrate structural stability depends on the presence of the guest molecules (empty hydrate is not stable), the influence of the cage occupancy on structural and energetic properties is also evaluated.

II. Calculations Details

The Vienna ab initio Simulation Package (VASP 5.3.5), a periodic plane wave basis set code [START_REF] Kresse | Ab Initio Molecular Dynamics for Liquid Metal[END_REF][START_REF] Kresse | Ab Initio Molecular Dynamics for Open-Shell Transition Metals[END_REF][START_REF] Kresse | Efficiency of ab Initio Total Energy Calculations for Metals and Semi-Conductors Using Plane-Wave Basis Set[END_REF][START_REF] Kresse | Efficient Iterative Scheme for ab Initio Total-Energy Calculations Using Plane-Wave Basis Set[END_REF] , has been used to perform DFT calculations. In order to calculate the exchange-correlation energy, two functionals were used: the semilocal PBE functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] and the nonlocal van der Waals Density Functional (vdW-DF) [START_REF] Dion | Van der Waals Density Functional for General Geometries[END_REF] .

Using both functionals made it possible to investigate how the dispersion forces influence the calculated properties. The projector augmented waves (PAW) method [START_REF] Blöchl | Projector Augmented-Wave Method[END_REF][START_REF] Kresse | From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method[END_REF] has been used to describe the electron core interaction. The energy cutoff in the plane wave expansion was 520 eV corresponding to high precision VASP calculations. All calculations are made at the Γ point of the Brillouin zone because of the large size of unit cells. where qSC and qLC are the small cage and large cage occupancies respectively. A summary of all the studied cases is presented in Table 1. The first calculations were performed on singly occupied systems, i.e. with all cages occupied by one guest only, which was positioned at the center of the cages (Fig. 2). These calculations are labelled [1/1] according to the above convention. Calculations have also been performed on systems with some doubly occupied large cages (i.e.

qLC > 1), keeping a single occupation of the small cages (i.e. qSC = 1). The double occupancy of the large cages has been varied from one large cage doubly occupied (i.e. 1 < qLC < 2) to all the large cages doubly occupied (i.e. qLC = 2). In the case of a doubly occupied large cage, the guest nitrogen molecules were positioned as follows: the molecular axis were positioned orthogonally at a distance of 2.4Å (corresponding to 1/3 of the LC diameter) from each other (Fig. 2) with the center of mass of the two molecules corresponding to the center of the cage. As an example, Fig. 2 illustrates singly and doubly occupied large

Single occupancy

Double occupancy

cages of sI and sII. Finally, various incompletely filled systems with empty cages were tested (i.e. qSC ≤1 or qLC ≤1].

Once the simulation box was built, a first relaxation of all atomic positions was performed to minimize the energy of the structure. Then a volume relaxation step was carried out in order to determine the optimum lattice parameter of each simulated system. A final global relaxation of the complete unit cell was done until forces on atoms converged to within 0.01 eV/Å. For each case studied, both sI and sII have been considered with the two functionals (PBE and vdW-DF). In order to verify our initial H-bonded water arrangement, a set of calculations has been performed with the initial structure arrangement provided by Takeuchi et al. [START_REF] Takeuchi | Water Proton Configurations in Structures I, II, and H Clathrate Hydrate Unit Cells[END_REF] The test was realized on simply occupied nitrogen hydrates (sI and sII) with both functionals used to run the calculations. All energy values obtained were very close to the ones calculated with the present built simulation box.

[0/1] 0 1 [0/1] 0 1 [0.5/1] 0.5 1 [0.5/1] 0.5 1 [1/0] 1 0 [1/0] 1 0 [1/0.5] 1 0.5 [1/0.5] 1 0.5 [1/1] 1 1 [1/1] 1 1 [1/1.17] 1 1.17 [1/1.125] 1 1.125 [1/1.33] 1 1.33 [1/1.25] 1 1.25 [1/1.5] 1 1.5 [1/1.375] 1 1.375 [1/1.67] 1 1.67 [1/1.5] 1 1.5 [1/1.83] 1 1.83 [1/1.625] 1 1.625 [1/2] 1 2 [1/1.75] 1 1.75 [1/1.875] 1 1.875 [1/2] 1 2

III. Results and Discussion

1) DFT-derived structural properties

Unit cell volume optimizations were performed for all the simulated systems. These optimizations were done by varying lattice parameters using a 0.2 Å step and relaxing all atomic positions. In Fig. 3 examples of calculated vdW-DF potential energy curves are presented as a function of the lattice parameter for both type I and type II structures and for different cage occupancies (PBE functional results are provided in Fig. S1). These potential energies were successfully fitted using the Birch-Murnaghan equation of state (denoted

EOS) 51 𝐸(𝑉) -𝐸 ! = "# ! $ " %& ()* # ! # + ' ( ) -1-
( 𝐵 * + + )* # ! # + ' ( ) -1- ' )6 -4 * # ! # + ' ( ) -2 (1)
where E0 corresponds to the minimum energy associated with the volume V0, BT is the bulk modulus and 𝐵 ! "

is the derivative of the bulk modulus with respect to pressure. The optimal lattice parameter of the cubic structure deduced from these fits was then used to perform the final relaxations and extract properties discussed in the following. The dependence of the unit cell parameter on the large cage and small cages occupancies is shown in Fig. 4 for sI and sII structures. The optimal lattice parameters obtained with the vdW-DF functional (empty symbols) are larger than the one obtained with the PBE functional (solid symbols), for both structures. Such a behavior has been reported in the case of methane [START_REF] Cox | Benchmarking the Performance of Density Functional Theory and Point Charge Force Fields in Their Description of SI Methane Hydrate Against Diffusion Monte Carlo[END_REF] and carbon monoxide 39 hydrates: the vdW-DF functional lead to overestimate the interatomic distances (and thus overestimate the lattice parameter) compared to PBE functional due to the repulsive nature of the exchange-correlation contribution of the vdW-DF functional [START_REF] Petuya | Unraveling the Metastability Of the SI and SII Carbon Monoxide Hydrate with A Combined DFT-Neutron Diffraction Investigation[END_REF] . The PBE-derived unit cell parameter increases with the cage occupancy, regardless of the considered structure and interestingly, a discontinuity in the lattice parameters slope is observed for qLC > 1 (i.e. for double occupancy of the large cages). With the vdW-DF functional, the same trend is observed for the sI structure while the sII does not exhibit significant variations with the cage occupancy.

The lattice parameters presented in Fig. 4, resulting from DFT calculations, correspond to values obtained at zero temperature (no kinetics energy contribution is considered in the calculations) and minimized pressure as relaxed lattice parameters were considered (i.e. P ~ 0 bar). For a more appropriated comparison of the DFT-derived cell parameter with the experimental ones recorded at higher pressures (see Fig. 5), the parameters of the Birch-Murnaghan EOS are used to include the pressure effects in the calculation of the DFT-derived cell parameters. The Birch-Murnaghan EOS (see equation 1) -used to fit the potential energy dependence with volume of Fig. 3 -gives access to the bulk modulus BT , i.e. the inverse of the relative volume variation of the material submitted to a pressure, expressed as:

𝐵 * = -𝑉 * ,- ,# + -.! (2)
The DFT-derived bulk moduli BT of the N2 hydrates (reported in Table 2) exhibit a general trend: PBEcalculated BT are greater than the ones derived from vdW-DF, except for multiple LC occupancy of the sI structure. The variation range of bulk modulus values is also larger for the PBE functional. The overestimated distances for the vdW-DF functional lead to larger dimension of the cages compared with the PBE functional. Therefore, a smaller impact on the cage occupancy is observed for the vdW-DF functional thus giving a smoother evolution of the bulk moduli.

In order to compare with the DFT-derived bulk moduli, the experimental values are reported in the Table 3.

The DFT-derived values are of the same order of magnitude as those measured experimentally with a better agreement for the most recent study. [START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF] Moreover, the sI bulk moduli are lower than those of sII for experimental results and for most cases in theoretical results. The lower BT value measured in the case of structure sI reflects the fact that this structure has a higher elasticity (easier deformation with pressure) than the structure sII. In other words, the sII requires more energy to be deformed. This behavior agrees with the present DFT-derived results shown in Fig. 3: curvatures of the potential energies variation with the lattice constant are less pronounced for sI than for sII -a flatter curvature meaning a lower bulk modulus BT. In addition, the volume variation, illustrated through the cage occupancy dependence of the cell parameter (Fig. 4), is more important for sI than for sII, with a pronounced deviation in the case of the PBE functional.

This supports the observation that the structure sI is more impacted than the structure sII by the cage occupancy variation. In some limited cases, the calculated BT values are lower for sII than for sI. This happens mostly when the small cages are not fully occupied such as in the [0/1] and [0.5/1] cases. Since there are more empty SCs in sII than in sI, the sII structure becomes more deformable than the sI structure. There are only few experimental studies reporting an isothermal variation of the cage occupancy and unit cell parameter with pressure. [START_REF] Kuhs | Cage Occupancy and Compressibility of Deuterated N2-Clathrate Hydrate by Neutron Diffraction[END_REF][START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF][START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF] From Rietveld refinement 8,10 the small cage occupancy qSC has been reported to vary from 85 to 100% while the large cage occupancy qLC is almost always greater than or equal to 100%. The isothermal compressibility coefficient kT (inverse of the bulk modulus) has been used to calculate the DFT-derived cell parameters at pressures matching the experimental values and at cage occupancies (indicated with stars in Table 2) close to the experimentally reported values, for which it was assumed a SC occupancy of 1 in all cases. The Fig. 5 shows the DFT calculated cell parameters versus the large cage occupancy for a full SC occupancy (qSC = 1) of the sI and sII structures. The experimental cell

LC SC

parameters determined from neutron diffraction are also reported on the same figure. As expected, the pressure-corrected theoretical data leads to lower lattice parameter for both functionals than the initially calculated at P = 0 bar. For pressures above 200 bar, the lattice parameter decreases with increasing pressure for both theoretical and experimental data, as expected from the compressibility of the material. It should be mentioned that the experimental values reported at 1 bar has been acquired at 80 K 13 while the other experimental data at higher pressures 8,10 were measured at 273 K. The pressure-corrected sI lattice parameter obtained with the vdW-DF functional is systematically larger than the experimental values while the PBE functional pressure-corrected lattice parameter is in excellent agreement. The same agreement has been reported in the case of carbon monoxide 39 and methane [START_REF] Cox | Benchmarking the Performance of Density Functional Theory and Point Charge Force Fields in Their Description of SI Methane Hydrate Against Diffusion Monte Carlo[END_REF] hydrates, comforting the idea that the hydrogenbonding network is correctly described with the PBE functional. In the case of the sII structure, both functionals lead to identical relative deviation (lower than 2.5%) from the experimental data. 2) Potential energy analysis and hydrate stability.

As reported from time-dependent neutron diffraction analysis, [START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF] the nitrogen hydrate forms in the sI clathrate metastable structure and transforms into the thermodynamically stable sII clathrate structure for which the cage occupancy plays a role [START_REF] Petuya | Ageing and Langmuir Behavior of the Cage Occupancy in the Nitrogen Gas Hydrate[END_REF] . Such structural metastability involves different potential energies for the sI and sII structures. The nitrogen hydrate stability has been analyzed in terms of DFT-calculated sI sII potential energies by following the methodology detailed in the work of Petuya et al. [START_REF] Petuya | Unraveling the Metastability Of the SI and SII Carbon Monoxide Hydrate with A Combined DFT-Neutron Diffraction Investigation[END_REF] The potential energy considered for the present analysis corresponds to intermolecular energies whose reference states are defined on the basis of the intramolecular energies of the isolated molecular species (water and nitrogen). The intermolecular cohesive energy provides information regarding the clathrate formation energy and corresponds to the non-bonding energy E NB . It is decomposed into three main components: the host-host (HH) interactions E HH , the guest-host (GH) interactions E GH and the guest-guest (GG) interactions E GG . E HH is associated with the cohesive energy of the empty clathrate structure and corresponds to the H-bonding interactions between water molecules constituting the host substructure. The E GH contribution, also known as binding energy, mainly corresponds to interactions between the guest molecules trapped inside the cages and the water molecules forming the cages. The guest-guest interaction E GG -whose the contribution should be negligible at least for single occupation -is included into the calculated E GH . These intermolecular energies have been computed for all sI and sII structures and cage occupancies reported in Table 1 using both PBE and vdW-DF functionals. To assess the influence of the cage occupancy variation on the calculated stabilization energies, representative results are presented as a function of the large cage occupancy in Fig. 6 for the vdW-DF functional that takes into account the dispersion forces. All calculated energies are reported in Table S1 and S2 for vdW-DF and PBE functionals, respectively. the binding energy becomes lower for sI than for sII. Above qLC ~ 1, the large cages start to be doubly occupied and the binding energy increases for both structures with a more pronounced evolution for the sI structure than for the sII structure. Since the reported E GH contains not only the water-nitrogen molecular interactions but also the guest-guest interactions, a destabilizing effect might be associated with the repulsive interactions becoming less negligible between two guests encapsulated within the same LC. The LC diameter being larger in sII than in sI, 1 the sII LCs better accommodate two guest molecules so that the sII structure becomes more favorable than the sI structure in terms of binding energy for qLC > 1.

As previously described, it is necessary to consider the water-water interactions (E HH ) in addition to nitrogen-water interactions (E GH ) to fully analyze the hydrate structural stability: the water sub-structure is the major constituent in a hydrate and thus represent the major contribution to the stabilization energy. The nonbonding energy E NB presented in Fig. 6 (right) shows again a difference on both sides of qLC ~ 1. Below this value, E NB is smaller for the sI structure than for the sII structure. For qLC ranging from ca.1 to ca. 1.5, sI and sII exhibit similar stability in terms of non-bonding energy and for qLC > 1.5, the sII structure is clearly stabilized compared to the sI structure. The sI structure becomes less stable while increasing the LCs double occupancy: the non-bonding energy increases whereas the opposite happens for the sII structure.

IV. Conclusion

The nitrogen hydrate structural stability has been investigated by means of DFT calculations, performed using two different functionals: the semi-local PBE functional and the non-local vdW-DF functional in order to investigate the role of the dispersion forces on the structural stability of the clathrate. As experimentally shown, [START_REF] Kuhs | Cage Occupancy and Compressibility of Deuterated N2-Clathrate Hydrate by Neutron Diffraction[END_REF][START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF][START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF] the nitrogen hydrate can be form in the sI and sII structures depending on the thermodynamics conditions. In the present investigation, both structures were considered by simulating complete sI and sII unit cells with periodic boundary conditions. Experimentally, the cage occupancy has been proved to vary with the pressure and temperature. [START_REF] Kuhs | Cage Occupancy and Compressibility of Deuterated N2-Clathrate Hydrate by Neutron Diffraction[END_REF][START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF][START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF] The present DFT approach has thus included the influence of the cage occupancy onto the nitrogen hydrate structural stability and shown its key role played onto the formed hydrate structure. For both sI and sII structures, the DFT pressure-corrected lattice parameters are in good agreement (relative deviation lower than 2.5%) with experimental data. [START_REF] Kuhs | Cage Occupancy and Compressibility of Deuterated N2-Clathrate Hydrate by Neutron Diffraction[END_REF][START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF][START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF] The significant dependence of the calculated lattice parameters with the large cage occupancy thus shows that the lattice parameter measurement could constitute an experimental probe of the cage occupancy in nitrogen hydrate. The sI structure exhibits a more pronounced cage occupancy dependence of the lattice parameter than the sII structure. By analyzing the potential energy dependence of the DFT-derived lattice parameters with the Birch-Murnaghan equation of state, the bulk modulii of both structures have been calculated for various cage occupancies and found to be in excellent agreement with those measured experimentally. [START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF] Moreover, the structural stability of the N2 hydrate has been studied through a detailed analysis of the non-bonding potential energy. This analysis reveals that the large cage occupancy plays a key role in the type of formed clathrate structure (sI or sII). The nitrogen hydrate is preferentially formed in the sI structure for partial large cage occupancy, while the sI structure becomes less stable compared to the sII structure while increasing the double occupancy of large cages. For large cage occupancy greater than ca. 1.5, the sII structure becomes more stable than the sI structure. Such a result outlines the difference of elastic properties between the sI and sII structures: the bulk modulus is larger for sII than for sI, according to the measurements and the present DFT calculations.

Such theoretical energy analysis is in agreement with the structural metastability of gas hydrates recently evidenced, as well as with previous DFT calculations. Nitrogen hydrate initially forms in the sI structure that transforms into the sII structure after few days. [START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF][START_REF] Petuya | Ageing and Langmuir Behavior of the Cage Occupancy in the Nitrogen Gas Hydrate[END_REF] Carbon monoxide hydrate exhibits similar structural metastability as nitrogen hydrate, although its transformation from sI to sII is much longer and takes several weeks. [START_REF] Zhu | Encapsulation Kinetics and Dynamics of Carbon Monoxide in Clathrate Hydrate[END_REF][START_REF] Petuya | Guest Partitioning in Carbon Monoxide by Raman Spectroscopy[END_REF] Experimentally, it has been shown that the LCs in both hydrates have the ability to catch and release guest molecules as the pressure and the temperature vary. [START_REF] Petuya | Guest Partitioning and Metastability of the Nitrogen Gas Hydrate[END_REF][START_REF] Petuya | Guest Partitioning in Carbon Monoxide by Raman Spectroscopy[END_REF][START_REF] Petuya | Ageing and Langmuir Behavior of the Cage Occupancy in the Nitrogen Gas Hydrate[END_REF] Together with previous DFT calculations performed on carbon monoxide hydrate, [START_REF] Zhu | Encapsulation Kinetics and Dynamics of Carbon Monoxide in Clathrate Hydrate[END_REF][START_REF] Petuya | Unraveling the Metastability Of the SI and SII Carbon Monoxide Hydrate with A Combined DFT-Neutron Diffraction Investigation[END_REF] the present investigation clearly underlines that the cage occupancy is the factor driving the stability, i.e. leading to the lower potential energy of the formed hydrate structure.

From partial occupancy to single occupancy of the LCs, the sI structure is the most stable structure, probably because the sI LCs are smaller than the sII LCs. [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF] As soon as the large cage occupancy has reached a threshold and starts involving double occupancy, the sII structure is energetically more favorable than the sI structure.

Both nitrogen and carbon monoxide sI hydrate are significantly destabilized with increasing the double occupancy of the large cages. The larger size of the sII LCs 1 facilitates the encapsulation of two guest molecules per cage because of higher repulsive forces for two guest molecules in sI-LC when compared to sII-LC. Finally, the calculations have evidenced the need to consider the water-water interactions to describe the structural stability of the hydrate. Indeed, the aqueous network may undergo some deformations when hosting more than one guest molecule within the cages. The more deformed the water structure, the higher the host-host energy. Since water constitutes the major part of the hydrate structure, its contribution to the stabilization energy of the system is of prime importance.

The present study allowed to evidence the influence of the cage occupancy on the structural and energy properties of the nitrogen hydrate. The importance of considering both water-water and water-guest interactions to evaluate the structural stability of the hydrate has been highlighted.

Figure 1 .

 1 Figure 1. Theoretical (left) and experimental (right) nitrogen hydrate phase diagrams. Black solid lines correspond to calculated nitrogen hydrate equilibrium lines adapted from Lundgaard, et al 12 . The blue line

Figure 2 .

 2 Figure 2. Representation of sI (Left) and SII(Right) LC for single N2 occupancy (Top) and double (Bottom) N2 occupancy.

Figure 3 .

 3 Figure 3. DFT potential energy difference E-E0, as a function of the lattice parameter from DFT calculations obtained with the vdW-DF functional for both sI (empty squares) and sII (filled circles) structures. The results are shown for all cages single-occupancied (left, [1/1]), for half of the LCs double-occupied (center, [1/1.5]) and for all LCs double-occupied (right, [1/2]). The continuous lines represent the fitted curves using the Birch-Murnaghan equation of state (see Equation (1)).

Figure 4 .

 4 Figure 4. DFT-derived lattice parameter as a function of the large (left) and small (right) cage occupancy, for sI (at the top) and sII (at the bottom) structures. The solid and empty symbols correspond to the values obtained with the PBE and vdW-DF functionals, respectively.

Figure 5 .

 5 Figure 5. DFT-derived (black symbols) and compressibility-corrected (colored symbols) lattice parameter as a function of the large cage occupancy for sI (Left) and sII (Right). The solid and empty symbols correspond to the values obtained with the PBE and vdW-DF functionals, respectively. The crosses and asterisks correspond to experimental values measured at 273 K at several pressures (colored-indicated on the figure) 8,10 and the plus signs, to experimental values 13 measured at 80 K and 1 bar.

Figure 6 .

 6 Figure 6. Guest-host energy E GH (left) and nonbonding energy E NB (right) as a function of the large cage occupancy qLC for sI (squares) and sII (circles) nitrogen hydrates, calculated with the vdW-DF functional.In all cases, the small cages are all simply filled.

Table 1 .

 1 Summary of the calculations performed. In each case, both PBE and vdW-DF functionals have been used.

	sI clathrate structure		sII clathrate structure
	Labelling	Occupancy	Labelling	Occupancy
	[𝛉SC/𝛉LC]	𝛉SC	𝛉 LC	[𝛉SC/𝛉LC]	𝛉SC	𝛉 LC

Table 2 .

 2 Theoretical bulk moduli extracted from DFT simulations data fitted with the Birch-Murnaghan equation of state. Values are given for sI and sII hydrates with both PBE and vdW-DF functionals. Standard deviation errors associated to the fitting procedure are indicated in parentheses. The values marked with stars are the ones used to correct the data shown in Fig. 5.

			PBE			vdW-DF
	Labeling	sI		sII		sI		sII
	[𝛉SC/𝛉LC]	BT (kbar)	BT (kbar)	BT (kbar)	BT (kbar)
	[0/1]	106.0 (0.6) 107.9 (1.3) 76.6 (0.7) 64.2 (3.3)
	[1/0]	102.9 (4.0) 102.6 (0.5) 60.9 (9.6) 74.6 (0.8)
	[0.5/1]	140.6 (10.1) 109.8 (1.6) 80.5 (1.4) 70.4 (1.1)
	[1/0.5]* 100.0 (1.1) 100.6* (0.1) 67.6 (3.3) 74.7* (1.9)
	[1/1]*	106.5* (0.4) 100.3* (0.2) 78.4* (4.9) 80.4* (1.4)
	[1/1.125]*	-	-	99.3* (0.3)	-	-80.8* (1.5)
	[1/1.17]* 79.6* (0.7)	-	-74.7* (1.9)	-	-
	[1/1.25]*	-	-	98.5* (0.3)			76.9* (0.7)
	[1/1.33]	84.0 (1.6)	-	-	75.5 (2.0)	-	-
	[1/1.375]	-	-	100.4 (0.6) -	-	74.0 (1.6)
	[1/1.5]	64.8 (1.7) 98.1 (1.6) 86.1 (0.6) 75.5 (1.2)
	[1/1.625]	-	-	86.2 (4.4) -	-	70.7 (1.5)
	[1/1.67]	82.1 (1.8)	-	-	83.0 (0.9)	-	-
	[1/1.75]	-	-	101.2 (0.5)			86.5 (0.7)
	[1/1.83]	68.1 (2.8)	-	-	82.6 (1.4)	-	-
	[1/1.875]	-	-	93.0 (3.3)	-	-	84.9 (0.9)
	[1/2]	58.9 (1.4) 98.1 (1.6) 79.3 (0.4) 86.6 (1.1)

Table 3 .

 3 Experimental values of the bulk modulus BT measured for N2 hydrates of type I and type II structures[START_REF] Kuhs | Cage Occupancy and Compressibility of Deuterated N2-Clathrate Hydrate by Neutron Diffraction[END_REF][START_REF] Chazallon | In-Situ Structural Properties of N2-, O2-and Air-Clathrates by Neutron Diffraction[END_REF] .

	Reference	Hydrate	T (K)	P (bar)	BT (kbar)
		N2 D2O SI	273	500 -2500	25 -42
	8
		N2 D2O SII	273	500 -2500	~115
		N2 D2O SI	273	750 -1000	79
	10
		N2 D2O SII	273	150 -1000	93

FUNDING SOURCES

This paper falls in the frame of the MI2C project funded by the French ANR "Agence Nationale de la Recherche" (No. ANR-15CE29-0016)..

ACKNOWLEDGMENT

Calculations have been performed thanks to the computing center "Mésocentre de Calcul" at the "Université de Franche-Comté, France".

ASSOCIATED CONTENT

Supporting Information.

The following files are available free of charge.

Additional DFT-derived data regarding the lattice parameters and the calculated energies. (PDF)

AUTHOR INFORMATION

The authors declare no competing financial interests.