Quentin C K Hazard

Jérémy S P Froidevaux

Natalie Yoh

Jonathan Moore

Juliana Senawi

Luke Gibson

Ana Filipa Palmeirim

Foraging guild modulates insectivorous bat responses to habitat loss and insular fragmentation in peninsular Malaysia

Keywords: Passive acoustic monitoring, Habitat fragmentation, Hydroelectric dams, Island Biogeography Theory, Land-bridge islands, Tropical forests AICc: Akaike Information Criterion Constant Frequency CI: Confidence Interval FM: Frequency Modulated FMqCF: Frequency Modulated quasi-Constant Frequency LF: Low Frequency LM: Linear Model NDVI: Normalised Difference Vegetation Index NMDS: Non-Metric Multi-Dimensional Scaling QCF: Quasi-constant Frequency

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

Introduction

Humanity currently faces a need to reconcile human population growth, increasing energy demands, and the decarbonization of that energy. In this context, hydropower is an increasingly appealing option, representing 73% of the renewable energy produced in the world (The World Bank, 2016). Yet, river damming is also a major driver of habitat loss and insular fragmentation across lowland forests [START_REF] Gibson | How Green is 'Green' Energy?[END_REF]. By flooding the lowland areas, dam construction often creates insular forest fragments matching the previous hilltops that become isolated within an inhospitable aquatic matrix [START_REF] Jones | Extinction debt on reservoir land-bridge islands[END_REF]. Although recent efforts have been made to understand the ecological consequences of hydropower [START_REF] Palmeirim | Emergent properties of species-habitat networks in an insular forest landscape[END_REF][START_REF] Terborgh | Maintenance of tree diversity in tropical forests[END_REF], few studies have targeted the Southeast Asian forests [START_REF] Jones | Extinction debt on reservoir land-bridge islands[END_REF]. Such understanding is therefore considered a priority for biodiversity conservation in the region [START_REF] Coleman | Top 100 research questions for biodiversity conservation in Southeast Asia[END_REF]. Species diversity persisting in insular forest fragments is typically affected by fragment size and isolation, which limit species population size and colonisation rates respectively [START_REF] Mcarthur | The theory of island biogeography[END_REF]. Edge effects arising from insularisation, namely increased exposure to windthrows and floristic transition towards light-wooded fast-growing pioneer assemblages, affect habitat quality within insular fragments, thereby influencing remaining species diversity (Benchimol & Peres, 2015a). Likewise, further human disturbances, including logging and fires, affect forest structure and thus habitat quality. Altogether, habitat quality can be assessed using the Normalized Difference Vegetation Index (NDVI), distinguishing between forest with denser canopy (higher NDVI values) and sparser canopy (lower NDVI values) [START_REF] Tovar | NDVI as indicator of degradation[END_REF].

Species responses to habitat loss and insular fragmentation may further vary between [START_REF] Palmeirim | Emergent properties of species-habitat networks in an insular forest landscape[END_REF] and within biological groups [START_REF] Brändel | Consequences of fragmentation for Neotropical bats: The importance of the matrix[END_REF]Meyer & Kalko, 2008), as influenced by particular species traits (Meyer et al., 2008;[START_REF] Palmeirim | Forest area predicts all dimensions of small mammal and lizard diversity in Amazonian insular forest fragments[END_REF]. For instance, the persistence of mid-sized mammal species on islands can be related to their swimming capacity (Benchimol & Peres, 2015b), while that of lizard species is dictated by their thermoregulation mode (Palmeirim et al., 2017). Yet, even persisting taxa show general decreases in abundance in insular forest fragments compared to surrounding mainland (e.g., bats: [START_REF] Gorresen | Landscape Responses of Bats to Habitat Fragmentation in Atlantic Forest of Paraguay[END_REF], birds: [START_REF] Yong | Do insectivorous bird communities decline on land-bridge forest islands in Peninsular Malaysia?[END_REF] and terrestrial and arboreal mammals: Benchimol & Peres (2021)). As such, while certain species able to use non-forest habitats might become overabundant [START_REF] Moore | Invasive rat drives complete collapse of native small mammal communities in insular forest fragments[END_REF], forest-dependent species become rarer or locally extinct [START_REF] Palmeirim | Small mammal responses to Amazonian forest islands are modulated by their forest dependence[END_REF]. Understanding the drivers of species response to forest insularisationconsidering both environmental and intrinsic species characteristicsallows more efficient management actions to be proposed, which is not trivial given the expansion of the hydropower sector across lowland tropical forests [START_REF] Couto | Global proliferation of small hydropower plantsscience and policy[END_REF].

Although habitat loss and insular fragmentation have been reported as important drivers of bat species' local extinction in the Neotropics [START_REF] Colombo | Functional diversity and trait filtering of insectivorous bats on forest islands created by an Amazonian mega dam[END_REF]Meyer & Kalko, 2008), little is known for Asia. In fact, the only study targeting such effects on insectivorous bats in this region highlighted the importance of island area especially for forest-dependent species in East China [START_REF] López-Bosch | Passive acoustic monitoring reveals the role of habitat affinity in sensitivity of sub-tropical East Asian bats to fragmentation[END_REF]. In addition, in a non-insular matrix setting in peninsular Malaysia, the diversity of insectivorous bats was impacted by forest area, with species-specific responses being modulated by their habitat affinity [START_REF] Struebig | Conservation value of forest fragments to Palaeotropical bats[END_REF]. Insectivorous bats emit echolocation calls to navigate their surroundings and locate food [START_REF] Schnitzler | From spatial orientation to food acquisition in echolocating bats[END_REF]. The characteristics of the calls produced, e.g., call shape, are adapted to a species' foraging preferences [START_REF] Denzinger | Bat guilds, a concept to classify the highly diverse foraging and echolocation behaviors of microchiropteran bats[END_REF][START_REF] Schnitzler | Echolocation by Insect-Eating Bats[END_REF]. For example, species adapted to foraging in the forest interior (forest foragers) use long constant frequency (CF) calls or very short, broadband, frequency modulated (FM) calls, adapted to particularly cluttered environments. Edge foragers use quasi-constant frequency (QCF) calls, or medium frequency calls (FMqCF) composed of an FM component followed by a short and quasi constant element (qCF), allowing them to locate and navigate between background features (medium frequency FM component), and to locate prey at an intermediate distance (qCF component). Open-space foragers use low frequency (LF) (<30 kHz) FMqCF calls with a narrow FM component and a long QCF component, enabling prey detection in vast empty spaces [START_REF] Denzinger | Bat guilds, a concept to classify the highly diverse foraging and echolocation behaviors of microchiropteran bats[END_REF][START_REF] Schnitzler | Echolocation by Insect-Eating Bats[END_REF]. Within foraging guilds, bat calls are not species specific: the calls of several species have evolved in a convergent way to respond to analogous environmental pressures, resulting in very similar calls in species facing analogous ecological conditions and therefore preventing species-specific identification in certain instances [START_REF] Gibb | Emerging opportunities and challenges for passive acoustics in ecological assessment and monitoring[END_REF][START_REF] Russo | Bats are still not birds in the digital era: echolocation call variation and why it matters for bat species identification1[END_REF]. To overcome this issue, bat calls are commonly classified into sonotypes, i.e., calls of similar shape and peak frequency [START_REF] Roemer | An automatic classifier of bat sonotypes around the world[END_REF].

Owing to the role of their foraging habitat in modulating their calls, bats of different foraging guilds are expected to respond antagonistically to habitat loss and fragmentation [START_REF] Denzinger | Bat guilds, a concept to classify the highly diverse foraging and echolocation behaviors of microchiropteran bats[END_REF][START_REF] Schnitzler | Echolocation by Insect-Eating Bats[END_REF]. In the aftermath of damming, once continuous forests are lost to smaller insular forest patches harbouring a gradient of vegetation structures: while cluttered forest represent optimal conditions for supporting diverse assemblages of forest bats, the water matrix can pose a serious obstacle to this guild's mobility, thus acting as a morpho-ecological filter [START_REF] Colombo | Functional diversity and trait filtering of insectivorous bats on forest islands created by an Amazonian mega dam[END_REF]. On the other hand, the newly created edges and adjacent open water surface might constitute high-quality habitat for edge and open-space foragers, as these vegetation structures are adapted to these guilds' morphological and acoustical adaptations [START_REF] Denzinger | Bat guilds, a concept to classify the highly diverse foraging and echolocation behaviors of microchiropteran bats[END_REF][START_REF] Schnitzler | Echolocation by Insect-Eating Bats[END_REF]. In this sense, forest-dependent bats, which have been identified as being of conservation priority in SE Asia due to the general decline of forested areas in this region [START_REF] Kingston | Research priorities for bat conservation in Southeast Asia: A consensus approach[END_REF], are expected to be restricted to continuous mainland forest and large, well-connected islands. Differentiating the response of these three guilds is therefore central to the implementation of adapted conservation measures.

Notwithstanding, the study of insectivorous bats in Southeast Asia, especially forest and edge foragers, has been largely impaired by the ability of most insectivorous species to avoid live trapping methods, such as harp traps and mist nets [START_REF] Kingston | Response of bat diversity to forest disturbance in Southeast Asia: Insights from long-term research in Malaysia[END_REF]. With the increasing affordability of low-cost acoustic devices, insectivorous bat surveys are becoming more accessible and reliable, further allowing for high replication [START_REF] Gibb | Emerging opportunities and challenges for passive acoustics in ecological assessment and monitoring[END_REF][START_REF] Hill | AudioMoth: Evaluation of a smart open acoustic device for monitoring biodiversity and the environment[END_REF].

Utilising such technological advances, here we provide the first assessment of the effects of habitat loss and insular fragmentation on insectivorous bats in a daminduced landscape in peninsular Malaysia. Using passive acoustic monitoring, we surveyed insectivorous bats in 26 forest islands and two mainland continuous forest sites. Across a gradient of habitat loss and insular fragmentation, we tested the effect of island size, isolation, and habitat quality (island shape and NDVI), at the following levels of bat diversity (1) overall assemblage, considering sonotype richness, activity and assemblage composition, (2) foraging guild, separately considering the activity of forest, edge and open-space foragers, and (3) sonotype, given each sonotype activity. We hypothesise that increased forest size and canopy density (i.e., higher NDVI values), and decreased isolation positively influence sonotype richness and activity. Yet, the different foraging guilds are hypothesised to show contrasting responses. In particular, forest foragers are expected to respond positively to forest size, canopy density, and connectivity to the mainland and surrounding forest patches, while open-space and edge foragers are expected to show opposite responses.

Material and methods

Study area

This study was conducted within the insular fragmented landscape of the Kenyir Lake and its surroundings in peninsular Malaysia. This artificial freshwater reservoir was formed in 1986 by the damming of the Kenyir river. The novel insular landscape occupies 260 000 ha and is composed of >340 islands ranging in size between 0.6 and 1428 ha embedded in the water matrix (Figure 1). Tropical humid forest on the islands and the adjacent mainland continuous forest are characterised by lowland and mid-elevation dipterocarp vegetation. The wide reservoir landscape, including the islands and surrounding mainland continuous forest, was subject to selective logging prior to damming [START_REF] Muhammad Yusuf | Environmental Issues in a Federation : The Case of Malaysia[END_REF][START_REF] Qie | Dung beetle assemblages on tropical land-bridge islands: Small island effect and vulnerable species[END_REF]. This practice is still allowed in some parts of the lake's catchment area, but not in Taman Negara National Park, which borders its south-eastern side [START_REF] Mariapan | Ecotourism Potentials in Malaysia[END_REF]. This region experiences a wet season between November and March, and a dry season between May and October. Annual precipitation varies between 2700 and 4000 mm annually [START_REF] Qie | Dung beetle assemblages on tropical land-bridge islands: Small island effect and vulnerable species[END_REF].

Study design and data collection

We selected 26 islands covering a range of sizes (min-max: 0.45 -167.3 ha) and distances from the mainland (135 -2748 m), in addition to two mainland sampling sites. This sampling strategy was set up to study the effects of forest size and isolation independently, i.e. maintaining a low correlation between these variables (r = -0.36 when considering island size and distance to mainland, r = -0.49 when considering log-transformed island size and distance to the mainland). Bat acoustic surveys were carried out between September 8 th and October 13 th 2019, using Audiomoth recorders [START_REF] Hill | AudioMoth: Evaluation of a smart open acoustic device for monitoring biodiversity and the environment[END_REF] set on a sampling rate of 384 kHz and the gain to the second setting ("med"). This sampling period has been selected to avoid the severe monsoons typical of this region, and thus scaping any rain-induced deterioration of the acoustic recording. We deployed one recorder in each sampling site and recorded for six hours, a sampling effort deemed adequate when aiming to detect bat populations' and assemblages' response to habitat modification [START_REF] López-Baucells | Optimizing bat bioacoustic surveys in humanmodified Neotropical landscapes[END_REF]. Recordings were divided into two time periods starting 30 minutes before sunset and ending 30 minutes after sunrise: from 18:00 to 22:00, and from 04:00 to 06:00 [START_REF] Hayes | Temporal Variation in Activity of Bats and the Design of Echolocation-Monitoring Studies[END_REF], covering the two peaks of bat activity at dusk and dawn [START_REF] Fenton | A technique for monitoring bat activity with results obtained from different environments in southern Ontario[END_REF]. Each recorder was attached to a tree, positioned 2 metres above the ground and, to minimise any uncontrolled impact from edge effects, placed as inland as possible relative to island size, i.e., between 14 and 123 m from the edge (median: 50 m).

Acoustic analysis

Using the software Kaleidoscope Version 5.4.7 (Wildlife Acoustics, 2019), we split the recorded sequences into 5-second recordings [START_REF] Kh Torrent | The importance of lakes for bat conservation in Amazonian rainforests: an assessment using autonomous recorders[END_REF]. The same software was used to filter the sequences containing sounds with a minimum frequency of 10 kHz and a maximum frequency of 250 kHz, and a pulse length between 2 and 500 ms. Among these sequences, only those containing one bat pass, i.e., at least two pulses of the same sonotype were kept for subsequent analysis [START_REF] Kh Torrent | The importance of lakes for bat conservation in Amazonian rainforests: an assessment using autonomous recorders[END_REF].

Prior to the acoustic analysis, we first compiled a list of all species of insectivorous bat known to occur in peninsular Malaysia [START_REF] Lim | Diversity of Malaysian insectivorous bat assemblages revisited[END_REF][START_REF] Nor Zalipah | Checklist of small mammals of Hulu Terengganu, Terengganu[END_REF]. Secondly, we collated reference calls for these species. We did so by conducting a literature survey using the Web of Science platform, between September and November 2021. We searched for publications by each species' name followed by the country name: we favoured reference calls obtained in our study area in order to avoid any potential geographical variation in the call parameters. For those species we could not find any reference calls, we used the reference calls available in the bat call library Chirovox [START_REF] Görföl | ChiroVox: A public library of bat calls[END_REF] (Supporting file 1). We therefore matched the call type of the species present in peninsular Malaysia to one of the sonotypes described in [START_REF] Yoh | A machine learning framework to classify Southeast Asian echolocating bats[END_REF] for bat species in Malaysian Borneo, namely CF, FM, FMqCF1, FMqCF2, FMqCF3, FMqCF4, FMqCF5, QCF and LF sonotypes. Using start and end frequency, frequency of maximum energy, duration and interpulse interval as defined in [START_REF] Russo | Identification of twenty-two bat species (Mammalia: Chiroptera) from Italy by analysis of time-expanded recordings of echolocation calls[END_REF], we classified the calls into one of these nine sonotypes. Given their very distinct echolocation parameters, the CF calls belonging to the genera Rhinolophus and Hipposideros could be identified to the species level.

As the shape of the echolocation calls reflects the physical constraints encountered by the bats, we were able to classify bat sonotypes into three foraging guilds: (1) the constant-frequency and FM calls represent forest foragers, [START_REF] Senawi | Kelawar Malaysia[END_REF] FMqCF4, FMqCF5 and QCF represent edge foragers, and (3) LF, FMqCF2 and

FMqCF3 calls represent open-space foragers [START_REF] Yoh | A machine learning framework to classify Southeast Asian echolocating bats[END_REF]. Social calls could not be identified to the sonotype level and were treated as assemblage-level activity.

Patch variables

Patch variables were obtained from a georeferenced LANDSAT 5 image which was transformed into a land/water matrix using an unsupervised classification on the software ArcGIS (ESRI, 2011). We then used the "landscapemetrics" R package [START_REF] Hesselbarth | landscapemetrics: an open-source R tool to calculate landscape metrics[END_REF] to extract: (1) island size (area; ha), ( 2) shortest Euclidean distance to the mainland (dist.main; m), [START_REF] Palmeirim | Small mammal responses to Amazonian forest islands are modulated by their forest dependence[END_REF] shortest Euclidean distance to the nearest neighbour island or mainland (dist.neigh; m), (4) island shape (shape), defined as the ratio between the patch perimeter and the hypothetical minimum perimeter of this patch, i.e. the perimeter of a maximally compact patch [START_REF] Mcgarigal | Comparative Evaluation of Experimental Approaches to the Study of Habitat Fragmentation Effects[END_REF], [START_REF] Palmeirim | Emergent properties of species-habitat networks in an insular forest landscape[END_REF] the normalised difference vegetation index (NDVI), and [START_REF] Palmeirim | Non-random lizard extinctions in land-bridge Amazonian forest islands after 28 years of isolation[END_REF] distance between the recorder and the forest edge (dist.edge; m). Given that the mainland continuous forest sites are characterised by an extensive forest coverage non-isolated area, we attributed these sites with the closest possible values to 'reality'. This included area values of one order of magnitude higher than the largest island (1670 ha) and zero distances to either the mainland or the nearest neighbour.

To streamline, we refer to the area of both islands and mainland size as 'forest area'.

NDVI and shape were calculated as for the remaining sampling sites but considering a 1000 m buffer centred in the sampling site and excluding water.

Data analysis

Assemblage-level metrics include sonotype richness, activity, and assemblage composition. Sonotype richness was defined as the number of sonotypes: this measure is representative of the diversity of call traits present at a site. Activity i.e., the number of bat passes, was used as a proxy for abundance. Although activity is widely used as a surrogate for abundance when studying echolocating bats (see for instance [START_REF] Charbonnier | Bat and bird diversity along independent gradients of latitude and tree composition in European forests[END_REF][START_REF] Ellerbrok | Activity of forest specialist bats decreases towards wind turbines at forest sites[END_REF], this method may overestimate abundance due to the possibility of multiple detection of a single individual [START_REF] Gibb | Emerging opportunities and challenges for passive acoustics in ecological assessment and monitoring[END_REF][START_REF] Kunz | Methods for assessing colony size, population size, and relative abundance of bats[END_REF]). Yet, it offers an effective tool for the detection of damming effects (e.g., [START_REF] López-Bosch | Passive acoustic monitoring reveals the role of habitat affinity in sensitivity of sub-tropical East Asian bats to fragmentation[END_REF][START_REF] Colombo | Functional diversity and trait filtering of insectivorous bats on forest islands created by an Amazonian mega dam[END_REF], as well as other human-induced disturbances on bat communities [START_REF] Mena | The effect of artificial light on bat richness and nocturnal soundscapes along an urbanization gradient in an arid landscape of central Peru[END_REF][START_REF] Williams-Guillén | Ensemble Composition and Activity Levels of Insectivorous Bats in Response to Management Intensification in Coffee Agroforestry Systems[END_REF]. Assemblage composition was summarised as a single variable using a Non-Metric Multi-Dimensional Scaling (NMDS) ordination. This analysis was performed considering sonotype activity and using a Bray-Curtis similarity matrix (stress = 0.130). The scores of the first axis of the NMDS composed the assemblage composition metric. Guild activity was calculated by summing the activity of the individual sonotypes respectively belonging to the forest, open-space and edge guilds (Table 1). Sonotype-level responses were examined for the sonotypes recorded in more than 10 sites and which had more than 50 bat passes. This threshold was intended to ensure a normal distribution of the residuals, as well as homoscedasticity. Social and unidentified calls were only included in the assemblage-level analysis. Among the sonotypes that met the threshold to be analysed, namely FMqCF2, FMqCF3, FMqCF4, FMqCF5, LF, R. trifoliatus, H. diadema and QCF, H. diadema and FMqCF2 had unequal error variances and were therefore excluded from the analysis.

We first accounted for spatial autocorrelation by applying Mantel tests using the R package "ade4" [START_REF] Dray | The ade4 package: Implementing the duality diagram for ecologists[END_REF]. These tests correlate geographic distance between sampling sites and each response variable as well as the residuals of each model introduced in the subsequent section. We found no spatial autocorrelation (p > 0.05) in all instances. We also examined the pairwise correlation between patch variables using Pearson correlation coefficients. Shape and area (log 10 x) (r = 0.720), as well as area (log 10 x) and NDVI were highly correlated (r = 0.800). Given the overall importance of area explaining biodiversity patterns in insular forest fragments [START_REF] Jones | Using relict species-area relationships to estimate the conservation value of reservoir islands to improve environmental impact assessments of dams[END_REF], we preferred to keep this metric to enable comparisons with other studies, whereas shape was excluded from subsequent analysis. Due to the lack of knowledge on the effects of canopy closeness on bats in this region, we chose to also keep NDVI in subsequent analyses. However, area (log 10 x) and NDVI were not included together in a model.

Collinearity between predictor variables was also examined using Variance Inflation Factors (VIFs), with no variable showing substantial collinearity (VIF>5) [START_REF] Dormann | Collinearity: a review of methods to deal with it and a simulation study evaluating their performance[END_REF].

We then analysed the combined effects of patch variablesarea, dist.main, dist.neigh and NDVIon (1) sonotype richness, activity and assemblage composition; (2) activity of forest, edge, and open-space sonotypes, and (3) the individual activity of eight sonotypes. To do so, we applied Linear Models (LMs) to each of these response variables, whose distribution was scrutinised prior to the analysis. The response variables regarding overall assemblage, guild and sonotype level activity, as well as forest area were log-transformed. Although all models were run with a gaussian error distribution, we initially considered a negative binomial distribution for the overall, guild-level and sonotype-level activity responses. Given that none of the models addressing individual sonotypes activity nor forest guild activity converged with a negative binomial error distribution, and that the distribution of these variables' residuals was closer to a normal distribution when using a logtransformation with a gaussian error structure, we chose to retain that transformation and error structure in the models. We further considered dist.edge as a covariate in each model, aiming to control for any eventual effect of distance to the forest edge.

A candidate model set including all possible combinations of patch variables (including the covariate dist.edge), except combinations involving area (log 10 x) and NDVI in the same model, was generated using the dredge function of the "MuMIn" R package [START_REF] Barton | Package "MuMIn" Title Multi-Model Inference[END_REF]. All models were ranked by Akaike Information Criteria corrected for small sample sizes (AICc: [START_REF] Burnham | Model Selection and Multimodel Inference[END_REF]. To account for model uncertainty in multi-model inference, we used a model-averaging approach considering the most parsimonious models, i.e. those having the lowest AICc within a ΔAICc <2 (ΔAICc = AICc i -AICc min , i being the i th model derived from the dredge) [START_REF] Froidevaux | Tree size, microhabitat diversity and landscape structure determine the value of isolated trees for bats in farmland[END_REF]. We report model average estimates along with their 95% confidence intervals (CIs) which were considered significant if not overlapping zero [START_REF] Nakagawa | Effect size, confidence interval and statistical significance: A practical guide for biologists[END_REF]. Assumptions about the normal distribution of the variables and their residuals were verified using the R packages "performance" [START_REF] Lüdecke | performance: An R Package for Assessment, Comparison and Testing of Statistical Models[END_REF] and "Dharma" [START_REF] Hartig | DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models[END_REF]. All data analyses were performed using R (R Core Team, 2022).

Results

In total, we recorded 21 197 bat passes from 16 different sonotypes: 10 forest, three edge and three open-space foragers (Table 1). Sonotype richness varied between 4 and 13 sonotypes per site, activity varied between 43 and 3351 bat passes per six hours recording. Activity varied greatly across sampling sites (43 -3351, 757.03 1 andS1). According to the NMDS, low values in the first axis were mostly associated with larger forest sites and forest foragers (H. cervinus, H. kunzi, H. bicolor, FM), while high values were associated with smaller islands, as well as with edge (FMqCF4, FMqCF5), open-space (FMqCF3), and forest foragers (R. refulgens, R. affinis) (Figure 2). Overall, 988 bat passes could not be identified to either the guild or to the sonotype level, including 981 bat passes corresponding to social calls (Table 1). R. trifoliatus was likely greatly influential on the response of the forest guild: being present in less than half the forest sites, it accounted for nearly 80% of the forest guild's activity (Table 1). Most of the FM sonotype activity (71.42%) was recorded on the mainland sites.

Overall assemblage responses

Sonotype richness increased with NDVI (β = 0.819 ± 0.409, p = 0.045, CI min = 0.017, CI max = 1.620), while total bat activity was unaffected by the patch variables considered (Table S2). Assemblage composition varied among sites based on their size (β = -0.536 ± 0.137, CI min = -0.817, CI max = -0.254), p < 0.001) (Table S2, Figure 3D, I and Figure 4).

Guild-level responses

Forest sonotypes were more active at sites with higher NDVI (β = 1.476 ± 0.439, p = 0.002, CI min = 0.616, CI max = 2.335). The activity of forest sonotypes also increased with decreasing distance to the closest neighbouring forest site (β = -1.468 ± 0.439, p = 0.002, CI min = -2.328, CI max = -0.608), while edge sonotypes activity decreased with increasing both island size (β = -1.050 ± 0.366, p = 0.004, CI min = -1.768, CI max = -0.332) and distances to edge (β = -1.045 ± 0.355, p = 0.003, CI min = -1.741, CI max = -0.350). None of the tested variables had a significant effect on open-space sonotypes (Table S2, Figure 3O -Q and Figure 4). Unlike all other response variables, only one model was selected for the activity of forest foragers (Dist.neigh, NDVI, AICc = 132.129) (Table S3).

Sonotype-level responses

Among all the six individual sonotypes, only FMqCF4 and R. trifoliatus showed a significant response to the tested patch variables. FMqCF4 sonotype showed higher activity in smaller islands (β = -1.127 ± 0.365, p = 0.002, CI min = -1.841, CI max = -0.412), while R. trifoliatus was more active on larger forest sites (β = 1.487 ± 0.446, p = 0.001, CI min = 0.614, CI max = 2.361) (Table S2, Figure 4 andS1).

Discussion

A number of studies have demonstrated that habitat loss and insular fragmentation cause species local extinctions across lowland tropical forests [START_REF] Gibson | Near-complete extinction of native small mammal fauna 25 years after forest fragmentation[END_REF][START_REF] Moore | Invasive rat drives complete collapse of native small mammal communities in insular forest fragments[END_REF][START_REF] Palmeirim | Emergent properties of species-habitat networks in an insular forest landscape[END_REF][START_REF] Pinto Henriques | Avian extinctions induced by the oldest Amazonian hydropower mega dam: Evidence from museum collections and sighting data spanning 172 years[END_REF]. Here, we contribute to fill an important knowledge gap by accordingly demonstrating overall negative bats response to dam-induced disturbance across an insular fragmented landscape in Southeast Asia. Our results highlight the role of canopy density driving the number of sonotypes, whereas forest area dictated which sonotypes were able to persist. Our guild-level analysis revealed that forest foragers were associated with denser forest structures, likely representing higher habitat quality for this guild, and were negatively affected by increasing isolation from neighbouring landmasses. In contrast, edge foragers seemed to benefit from island shrinkage. Fragmentation effects were not so clearly observed at the sonotype-level, with only two of the six sonotypes analysed responding to patch variables, namely to forest area which had a positive effect on the forest forager R. trifoliatus and a negative effect on the edge sonotype FMqCF4.

Effects of forest area

Forest area did not predict sonotype richness at our study site. Owing to their high correlation, habitat quality and forest area were precluded from being combined in models, and models including habitat quality were more parsimonious than those including area: the absence of area effects is therefore more likely a consequence of the study design rather than a proper lack of area effects. In fact, a large body of evidence reports that area affects bat richness in fragmented insular [START_REF] Luypaert | Extending species-area relationships into the realm of ecoacoustics: The soundscapearea relationship[END_REF], and non-insular systems [START_REF] Rocha | Consequences of a large-scale fragmentation experiment for Neotropical bats: disentangling the relative importance of local and landscape-scale effects[END_REF]. This pattern holds true for several taxa on reservoir islands [START_REF] Palmeirim | Emergent properties of species-habitat networks in an insular forest landscape[END_REF], including dung-beetles [START_REF] Qie | Dung beetle assemblages on tropical land-bridge islands: Small island effect and vulnerable species[END_REF], primates and ungulates [START_REF] Yong | Persistence of primate and ungulate communities on Forested Islands in Lake Kenyir in Northern peninsular Malaysia[END_REF], and birds [START_REF] Yong | Do insectivorous bird communities decline on land-bridge forest islands in Peninsular Malaysia?[END_REF] surveyed in nearly the same islands in Kenyir. Notwithstanding, bat assemblage composition varied along the gradient of forest area, with edge foragers being particularly active on smaller islands. This trend was further reflected at the sonotype level by the edge forager FMqCF4 and the forest forager R. trifoliatus, both of whom displayed contrasting responses to forest area, the former being negative, and the latter being positive. These responses were expected given that small islands tend to be edge-dominated. The extensive variation in species-specific home-range size may have further influenced the ability of some species to use the smaller fragments: while that of most edge foragers (5 -210 ha) and R. trifoliatus (5 ha) are notably small, the majority of forest foragers have a vast home-ranges, reaching up to 2199 ha for some species [START_REF] Wilson | Conserving biodiversity in production landscapes[END_REF]. Given their tendency to avoid flying across the water matrix (Meyer & Kalko, 2008), forest foraging species, mostly belonging to the FM sonotype, were therefore rarely recorded on smaller islands. Yet, as a caveat to this study, in small islands, detectors had to be placed closer to edges given the lack of forest interior. It is therefore possible that the detectors on small islands, being mechanically closer to forest edges, recorded a higher activity of edge foragers. This is further supported by the negative relationship between edge foragers activity and the Dist.edge. In any case, this would still demonstrate the preferential use of edges by this bat guild [START_REF] López-Bosch | Passive acoustic monitoring reveals the role of habitat affinity in sensitivity of sub-tropical East Asian bats to fragmentation[END_REF]. In contrast, our results show that the forest forager R. trifoliatus responds positively to forest area, suggesting that this sonotype requires greater habitat complexity associated with larger areas of forest (Benchimol & Peres, 2015a).

Effects of habitat quality

Despite the absence of clear effects of forest size, canopy closeness, as indicated by the NDVI, promoted an increase in the number of bat sonotypes and forest foragers activity across the Kenyir landscape. A bat species' response to habitat quality is likely influenced by the intrinsic habitat characteristics such as 3D forest structure or canopy ruggedness, ultimately impacting which species are able to use each site [START_REF] Froidevaux | From field surveys to LiDAR: Shining a light on how bats respond to forest structure[END_REF]. Our study landscape has been subject to intensive selective logging prior to the construction of the dam [START_REF] Qie | Dung beetle assemblages on tropical land-bridge islands: Small island effect and vulnerable species[END_REF]. Indeed, evidence for previous logging is still noticeable in the area, with logging trails and canopy holes being observed in-situ, resulting in a low but variable overall NDVI on the islands and surrounding continuous forest sites. Whilst the effects of logging on bat species richness seem to be limited both in the Neotropics [START_REF] Meyer | Responses of tropical bats to habitat fragmentation, logging, and deforestation[END_REF] and in the Paleotropics [START_REF] Struebig | Quantifying the Biodiversity Value of Repeatedly Logged Rainforests. Gradient and Comparative Approaches from Borneo[END_REF], logging appears to strongly influence assemblage composition, edge species being indicative of repeatedly logged sites [START_REF] Peters | Effects of selective logging on bat communities in the southeastern Amazon[END_REF]). Yet, the effects of logging on biodiversity depend on the intensity and extraction methods [START_REF] Burivalova | Thresholds of Logging Intensity to Maintain Tropical Forest Biodiversity[END_REF], and further investigations regarding the effects of logging intensity in the context of insular forest fragments are needed to further our understanding of how logging may drive bat sonotype richness.

Notwithstanding the potential effects of logging, forest insularisation led to the creation of edges, whose deleterious effects on vegetation include increased exposure to wind-throws, culminating in shifts towards disturbance-adapted pioneer trees (Benchimol & Peres, 2015a;[START_REF] Santo-Silva | Phylogenetic homogenization of Amazonian tree assemblages in forest islands after 26 years of isolation[END_REF]. While sites with low canopy closeness can be widely used by edge foragers, only those sites harbouring increased NDVI may represent suitable habitat for manoeuvrable forest dependent species that are further adapted to echolocate in more cluttered environments [START_REF] Froidevaux | From field surveys to LiDAR: Shining a light on how bats respond to forest structure[END_REF][START_REF] Suarez-Rubio | Insectivorous bats respond to vegetation complexity in urban green spaces[END_REF]. By allowing forest foragers to persist, habitat quality contributes to maintain bat diversity, as also observed for other biological groups, e.g., large-sized mammals and reptiles [START_REF] Oliveira | Effects of landscape and patch attributes on the functional diversity of medium and large-sized mammals in the Brazilian Cerrado[END_REF][START_REF] Silva | Habitat Quality, Not Patch Size, Modulates Lizard Responses to Habitat Loss and Fragmentation in the Southwestern Amazon[END_REF]. This is further supported by the increase in forest bat activity we observed in denser canopies. Higher NDVI values may also be associated with higher availability of mature trees that provide roosting sites for species such as R. trifoliatus, R. sedulus, K. papillosa and K. pellucida, all of which depend on these structures to rest and thus to persist. For instance, in Malaysia, the absence of tree cavities due to forest disturbance was associated with the decline of the forest foragers Kerivoula sp. [START_REF] Struebig | Quantifying the Biodiversity Value of Repeatedly Logged Rainforests. Gradient and Comparative Approaches from Borneo[END_REF]. Our findings reiterate the importance of habitat quality as a key driver of species diversity in fragmented landscapes [START_REF] Armstrong | Using experimental reintroductions to resolve the roles of habitat quality and metapopulation dynamics on patch occupancy in fragmented landscapes[END_REF][START_REF] Poniatowski | Patch occupancy of grassland specialists: Habitat quality matters more than habitat connectivity[END_REF].

Effects of isolation

Contrary to our expectations, isolation was not an important variable explaining bat assemblage-level responses. These results contrast with an insular fragmented landscape in Panama, where isolation to the mainland was the main predictor of bat richness (Meyer & Kalko, 2008). However, our results are in agreement with findings from a non-insular Malaysian fragmented landscape, where isolation has also been found to be a poor prediction of bat richness [START_REF] Brändel | Consequences of fragmentation for Neotropical bats: The importance of the matrix[END_REF]. The lack of isolation effects may be related to the overall small distance separating most of the sampling sites, and the overall size of the lake, as the home ranges for most local species exceed the distance separating most of the study sites [START_REF] Wilson | Conserving biodiversity in production landscapes[END_REF].

Nevertheless, forest foragers were more active in sites less isolated from neighbouring landmasses, which might be due to morphological constraints [START_REF] Norberg | Ecological morphology and flight in bats (Mammalia; Chiroptera): wing adaptations, flight performance, foraging strategy and echolocation[END_REF]. Indeed, forest foragers have a wing morphology characterised by a low aspect ratio (wingspan²/wing area) and a low wing loading (body mass/wing area) [START_REF] Norberg | Ecological morphology and flight in bats (Mammalia; Chiroptera): wing adaptations, flight performance, foraging strategy and echolocation[END_REF]. Although this characteristic allows them to have a slow and highly manoeuvrable flight, it also makes flight over open spaces particularly energetically demanding [START_REF] Altringham | Bats: From Evolution to Conservation[END_REF][START_REF] Bader | Mobility explains the response of aerial insectivorous bats to anthropogenic habitat change in the Neotropics[END_REF].

Furthermore, the absence of distance-to-mainland effects in favour of distance-toneighbour effects for forest foragers underlines the value of intermediary islands to act as stepping-stones for forest bats to cross the water matrix and reach more remote islands. This idea is supported by [START_REF] Saura | Stepping stones are crucial for species' long-distance dispersal and range expansion through habitat networks[END_REF] who also stress that these intermediate islands need to be sufficiently large and of high quality in order to act as stepping-stones.

Limitations and potential caveats

Our results emphasise the valuable use of passive acoustic monitoring techniques to survey bat assemblages, further allowing us to examine bat responses at multiple levels. However, the use of sonotype richness instead of species richness likely biased the estimated sonotype richness towards forest species. Indeed, while the CF calls produced by forest foragers could be identified to the species level, other sonotypes including FMqCF, QCF and LF contained multiple species. Likewise, given the similarity of the ecological constraints faced by edge foragers, the calls produced by the species belonging to this guild can only be separated between three sonotypes (FMqCF4, FMqCF5 and QCF). For this reason, our results might underestimate the effects of habitat loss and insular fragmentation on insectivorous bats, which is further enhanced by the fact that larger forest sitesexpected to harbour higher species diversitywere proportionally less sampled. In addition, as species detectability is a function of call intensity [START_REF] Hayes | Assumptions and practical considerations in the design and interpretation of echolocation-monitoring studies[END_REF], forest bats producing low-intensity FM calls such as Vespertilionidae (e.g., Kerivoula and Myotis spp.) tend to be under-detected [START_REF] Waters | Echolocation call structure and intensity in five species of insectivorous bats[END_REF]. This might further explain the relatively weak responses observed at the sonotype-level, which should therefore be interpreted with caution. Likewise, even the most commonly used devices in bat detection, including Audiomoths, tend to lack the sensitivity required to detect some of the high intensity calls emitted by smallest-bodied bats such as H. cervinus, H. larvatus and H. bicolor [START_REF] Kingston | Research priorities for bat conservation in Southeast Asia: A consensus approach[END_REF]: although common in Peninsular Malaysia [START_REF] Lim | Diversity of Malaysian insectivorous bat assemblages revisited[END_REF], these species have probably been under detected in our study. Live trapping remains the most efficient method to monitor these species [START_REF] Kingston | Response of bat diversity to forest disturbance in Southeast Asia: Insights from long-term research in Malaysia[END_REF], and still, studies using these trapping methods have highlighted the high sensitivity of these forest genera to forest disturbance (Huang et al., 2019).

Conservation implications

SE Asia may lose over 74% of its original forest cover by the end of the century, putting forest-dependent species at tremendous risk [START_REF] Sodhi | Southeast Asian biodiversity: An impending disaster[END_REF]. Forest foragers from the Kenyir lake are no exception: being absent from the smallest islands, they showed low activity in isolated and degraded forest patches, which nevertheless make up the majority of the landscape of Lake Kenyir. Additionally, three forest species, namely H. bicolor, H. cervinus and H. kunzi, were only found on the mainland. Our results suggest that conservation efforts should target forest bats which, given their forest-adapted morphology and the rapid vanishing of their foraging habitat, are particularly extinction prone [START_REF] Jones | Biological Correlates of Extinction Risk in Bats[END_REF][START_REF] Safi | A Comparative Analysis of Specialization and Extinction Risk in Temperate-Zone Bats[END_REF]. The FM sonotype, of which most potential representatives such as Kerivoula. intermedia, K. pellucida or Nycteris tragata are listed as near threatened [START_REF] Senawi | Kelawar Malaysia[END_REF], was mostly active in the mainland, and completely absent from 21 islands. We therefore stress that bats highly dependent on forest can only persist in large undisturbed forest tracts. The independent presence of dense canopies and connection to nearby landmasses does not guarantee the use by forest bats. Indeed, only patches that are large, well connected to the mainland, and harbouring a high habitat quality can serve as stepping-stones, and therefore allow less vagile species to commute over the water matrix [START_REF] Saura | Stepping stones are crucial for species' long-distance dispersal and range expansion through habitat networks[END_REF]. In insular forest patches, species are lost in a sustained and delayed manner according to the time elapsed since isolation, a process referred to as an "extinction debt" [START_REF] Jones | Extinction debt on reservoir land-bridge islands[END_REF]. Furthermore, these islands suffer an inevitable vegetation decay (Benchimol & Peres, 2015a). In line with [START_REF] Jones | Extinction debt on reservoir land-bridge islands[END_REF], our results suggest that islands should not be used as conservation units for mitigation purposes due to the imminent extinction debt. Instead, conservation efforts should prioritise maintaining mainland habitat quality, for instance by minimising logging activity in these highly forested areas [START_REF] Poudyal | Evolutionary dynamics of selective logging in the tropics: A systematic review of impact studies and their effectiveness in sustainable forest management[END_REF][START_REF] Harvey | Effects of mechanized careful logging on natural regeneration and vegetation competition in the southeastern Canadian boreal forest[END_REF]. This can be achieved by legally protecting the islands and surrounding mainland continuous forest/wider reservoir landscape as to minimise disturbances. Such measures have successfully helped minimising further human activities in a Brazilian Biological reserve encompassing part of a mega-dam (Benchimol & Peres, 2015a, 2015b). In addition, future hydropower developments should consider how dam placement is likely to affect the creation of different island systems. These developments should aim to reduce the creation of a myriad of small, isolated, and habitat-degraded forest fragments, for instance by targeting craggy locations over flat areas, therefore drastically minimising the flooded area. The biological impacts of damming tend to be largely underestimated, as current environmental impact assessment methods poorly predict the extent and location of dam-induced flooding, possibly leading to a 64.5 % increase of the flooded surface compared to initial predictions [START_REF] Cochrane | Landsat-based analysis of mega dam flooding impacts in the Amazon compared to associated environmental impact assessments: Upper Madeira River example 2006-2015[END_REF]. Working towards a more accurate planning is therefore essential to minimise the flooding area, avoid the flooding of ecologically valuable zones, and adequately locate future dams.

Conclusions

Hydropower development is set to massively expand across Southeast Asian forests, with energy production expected to increase threefold by 2035 [START_REF] Petinrin | Renewable energy for continuous energy sustainability in Malaysia[END_REF][START_REF] Tang | Current and future hydropower development in Southeast Asia countries (Malaysia, Indonesia, Thailand and Myanmar)[END_REF]. In Malaysia alone, at least four additional major dams will soon be constructed (> 34 000 MW) [START_REF] Foo | A vision on the opportunities, policies and coping strategies for the energy security and green energy development in Malaysia[END_REF]. Coupled with a steadily declining share of forest in this region, these threats put insectivorous bat assemblages at risk, especially for forest-dependent species. Here, we showed that canopy density played a central role in promoting both sonotype richness, persistence and activity of forest foragers, while this guild was negatively impacted by isolation to neighbouring forest patches. Forest area further affected the assemblage composition, larger forest patches being associated with forestdominated assemblages, and smaller patches being associated with edge-foragers.

Large, dense, well connected forest fragments still supported a subset of the mainlands' assemblage diversity. Yet, our results suggested that preserving a high canopy density, especially in the mainland, is essential to serve the preservation of forest-dependent species. Keeping forest disturbances such as logging minimum and developing more efficient tools to predict the spatial extent of future dam's effects is therefore a priority for the preservation of bat assemblages in SE Asia.

Figures and Tables:

Colour should be used for figures [START_REF] Yoh | A machine learning framework to classify Southeast Asian echolocating bats[END_REF][START_REF] Senawi | Kelawar Malaysia[END_REF][START_REF] Palmeirim | Small mammal responses to Amazonian forest islands are modulated by their forest dependence[END_REF]4 1). Given that the sonotypes H. kunzi, H. cervinus and H. bicolor are overlapping, for the sake of clarity, the position of these sonotypes is replaced by a green line and the sonotypes labels are separated.

Single column fitting image

Two columns fitting image

Table 1. List of all sonotypes identified during the acoustic bat surveys across the insular fragmented landscape of the Kenyir Lake, peninsular Malaysia. For each sonotype, we indicate the potential species matching that sonotype, corresponding foraging guild, total activity (number of bat passes), and number of sampling sites in which each of the sonotypes was recorded. Bat passes that could not be identified were labelled as "unknown".

Sonotype

  744.18), and among foraging guilds (0 -689, 101.89 174.73 for forest foragers, 2 -2923, 489.39 698.75 for edge foragers, and 11 -641, 130.39 144.01 for openspace foragers). While the edge forager FMqCF4 and the open-space forager LF were present at every site, the following forest foragers were found at only one site: Rhinolophus refulgens (island 13), CF.46 (island 25), Hipposideros cervinus, H. kunzi and H. bicolor (CF A) (Tables

Figure 1 .

 1 Figure 1. Map of the Kenyir Lake (A) Location of the study area and sampling sites in the Kenyir lake, peninsular Malaysia. The solid dots depict each of the 28 sampling sites. (B) Distribution of island sizes in Kenyir lake. (C) Photo of the Kenyir Lake. Two columns fitting image

Figure 2 .

 2 Figure 2. Non-Metric Multi-Dimensional Scaling (NMDS) ordination plot denoting both sampling sites and sonotypes. Sampling sites are represented by circles, matching the islands which are sized proportionally to their size (log 10 x), and triangles correspond to the mainland continuous forest sites. Sonotypes are represented by their name and colour-coded according to the corresponding foraging guild: forest (in green), edge (red) and open-space (blue) (for further details regarding each sonotype, see Table1). Given that the sonotypes H. kunzi, H.

Figure 3 .

 3 Figure 3. Results of the model averaging procedure. Response variables are shown on the y axis, and predictor variables are shown on the x axis. The variables that were retained in the model selection and averaging appear at the intersection between response and predictor variables. The value of each estimate is depicted with a colour gradient, negative estimates being increasingly red, and positive estimates being increasingly blue. The significance of each selected variable is shown with point shape: significant variables are represented with a triangle, and non-significant variables are shown with a circle. The relative importance of each selected variable within the averaged model is shown with point size, smaller points having a lower relative importance than larger points.Two columns fitting image

Figure 4 .

 4 Figure 4. Relationships between bat sonotype richness (A-D), activity (log 10 x ) (E-H), assemblage composition (axis 1 of the NMDS) (I-L), and the activity of forest (log 10 x) (M-P), edge (log 10 x) (Q-T), and open-space foraging bats (log 10 x) (U-X) and Area (log 10 x) (A, E, I, M, Q, U), distance to the mainland (Dist.main) (B, F, J, N, R, V), Distance to neighbour (Dist.neigh) (C, G, K, O, S, W) and NDVI (D, H, L, P, T, X). In each panel, the solid black line is the prediction given by the averaged model obtained from the dredge, and the shaded area represents the 95% confidence interval. The predictions of the selected model were only shown for significant variables.

This species has not been sampled around our study site to our knowledge

Acknowledgements:

We are grateful for the fieldwork assistance of David L´opez-Bosch and Ahmad Faizul Bin Zulkifli during data collection. We thank the Economic Planning Unit, Department of Prime Minister, Malaysia for permission to conduct research and Department of Wildlife and National Parks Peninsular Malaysia for permission to work in Kenyir (JPHL&TN (IP):100-34/1.24 Jld 14(57). AFP was supported by the Outstanding Postdoctoral Fellowship of the Southern University of Science and Technology (SUSTech), and is currently funded by the European Union's


Ruel Cedeno

Romain Grossier
email: romain.grossier@cnrs.fr

Nadine Candoni

Nicolas Levernier

Stéphane Adrian Flood
email: adrian.flood@vistec.ac.th

Stéphane Veesler
email: stephane.veesler@cnrs.fr

time distributions in contracting microdroplets

CNT effective interfacial energy and pre-exponential kinetic factor from measured NaCl crystal nucleation

INTRODUCTION

Nucleation plays a key role in a wide array of applications including nanosynthesis 1 , energy storage 2 , pharmaceutical production, biomineralization, and climate modeling. 3 Sodium chloride, being the most abundant salt on earth 4 , is of particular interest due to its influence on metal corrosion 5 , building material degradation 6 , oil well productivity 7 , atmospheric science 3 and so on. Thus, fundamental understanding of its nucleation kinetics is of paramount importance yet it remains poorly understood from both experimental and theoretical perspective. 8 Indeed, nucleation rates from simulations still differ from experiments by several orders of magnitude. 8,9 With diverse simulation approaches being applied, 8,[10][11][12][13][14] reliable benchmarking of computational results with experiments remains a challenge as there are only very few experimental studies that quantitatively measure the nucleation kinetic parameters of NaCl in water. These include experiments which employ an efflorescence chamber 15 , electrodynamic levitator 16 , and microcapillaries 4 ; all of which treated nucleation deterministically. In this context, deterministic methods calculate nucleation rates directly from the mean nucleation time whereas stochastic methods employ the probability distribution of nucleation times. However, primary nucleation has been shown to be inherently stochastic rather than deterministic. [START_REF] Deck | Crystal Growth & Design[END_REF][START_REF] Candoni | Advances in Organic Crystal Chemistry: Comprehensive Reviews[END_REF] Moreover, because nucleation is a rare event and the critical nucleus is a transient species, it is barely undetectable by classical experiment. This can be addressed by introducing a local bias in the solution so that nucleation is much more probable. [START_REF] Grossier | [END_REF] Recently, with in situ electron microscopy, Nakamuro et al. 20 have captured atomically-resolved images of NaCl nucleation in confined conical carbon nanotubes (with volume confinement as a bias). They observed that a critical cluster must have at least 48 NaCl units and that the nucleation periods follow a normal distribution spanning from 2 to 10 s based on nine nucleation events. Although more data points are needed to generate a reliable statistical distribution, this is a strong evidence for the stochasticity of NaCl nucleation yet surprisingly, there are no existing experimental studies that measure its kinetic parameters using the stochastic view of nucleation. In this communication, we address this by measuring the primary nucleation kinetic parameters of aqueous NaCl in confined microdroplets, in the pL range, with a stochastic model. We demonstrate that by combining the deliquescence-recrystallization cycle for measuring induction times 21 , an appropriate evaporation model 22 , together with an inhomogeneous Poisson probability distribution of nucleation time 23 and Classical Nucleation Theory (CNT), one can obtain a reliable estimate of nucleation parameters, effective interfacial energy 24 and pre-exponential kinetic factor, which are consistent with theoretical and experimental values from literature. Finally, we analyze the impact of confinement by volume 25,26 on the measured kinetic parameters.

II. MATERIALS AND METHODS

Arrays of sessile NaCl microdroplets are generated on PMMA-coated glass immersed in a thin layer of PDMS oil at ambient conditions (1 atm, 25 o C). Supersaturation to achieve nucleation is obtained via the droplet contraction technique 27 (similar to an evaporation process) under controlled humidity (10% RH) using the method described in our previous work. 21 A schematic overview of the setup and illustration of the droplet contraction technique is shown in Figure 1a-b while selected microdroplet images are shown in Figure 1c-d. For each droplet, we calculate the dimensionless nucleation time 21 τ i.e. the time elapsed between saturation and the nucleation point divided by a characteristic time accounting for the evaporation rate. 21 This is a pre-processing procedure to robustly take into account the variability of evaporation rate in the set-up. 28 We also verified that there is no interference of diffusion-mediated interactions between microdroplets. 21 To find the volume and supersaturation ratio as a function of τ, we employed an evaporation model tailored to our specific configuration. 22 We then plot the cumulative probability distribution of supersaturation ratio at nucleation (i.e. fraction of the microdroplets that has nucleated at a given supersaturation ratio) for each set of microdroplet sizes. We then fit distribution of nucleation times with the help of an inhomogeneous Poisson distribution, as proposed by Goh et al. 29 , under conditions of time-varying supersaturation

𝑃(𝑡) = 1 -exp [-∫ 𝐽(𝑡)𝑉(𝑡)𝑑𝑡 𝑡 nuc 𝑡 sat ] (1)
with P(t) as the fraction of microdroplets nucleated after time t, V(t) droplet volume at time t, tsat time at which the microdroplet becomes saturated and tnuc nucleation time.

For the primary nucleation rate J(t), we used the classical nucleation theory (CNT) as

𝐽(𝑡) = 𝐴 exp [ - 16𝜋 3 𝛾 𝑒𝑓𝑓 3 𝜌 𝑠 2 (𝑘 𝑏 𝑇) 3 ln 2 ( 𝛾 ± 𝛾 ± 0 𝑆(𝑡)) ] (2)
where A is the pre-exponential factor, γeff is the effective interfacial energy, ρs is the number density of formula units in the solid (2.27 × 10 28 m -3 for NaCl), kbT is the thermal energy, S(t) is the supersaturation expressed as the ratio of concentrations (concentration at nucleation / concentration at saturation) at nucleation at time t ,and 𝛾 ± and 𝛾 ± 0 are the ionic activities at S(t)

and at saturation respectively (see section S1-S3 of SI). This represents a general equation for primary nucleation with γeff accounting for the "degree" of heterogeneity. While analytical solutions for equations ( 1) and ( 2) exists for special cases, 30 our evolution of supersaturation ratio S(t) with time is governed by a system-specific evaporation model 22 (see Section S4 of SI) whose complexity requires a numerical approach (see Section S5 of SI). Furthermore, thanks to the small droplet size coupled with hindered evaporation rate (due to the oil layer), diffusion dominates over convection (Peclet number < 10 -3 , see Section S4.3 of SI). As a result, the microdroplets maintain a uniform concentration distribution 22 , instead of being concentrated at the interface or contact line (in contrast to µL droplets directly evaporating in air). 31

III. RESULTS AND DISCUSSION

With the nucleation times obtained via image analysis, we computed the normalized nucleation time τ and the results are plotted in Figure 2a. Using our tailored evaporation model, τ can then be converted into distributions of supersaturation ratios at nucleation, as presented in Figure 2b. As Figure 2b shows, the combination of the inhomogeneous Poisson distribution (eq 1) and classical nucleation theory (eq 2) well captures the sigmoidal nature of the distribution. Unlike the use of empirical distributions (such as Weibull, Gompertz, Gumbell etc), whose parameters cannot be interpreted in terms of CNT, this method allows the extraction of the kinetic parameter A (pre-exponential factor) and thermodynamic parameter γeff (effective interfacial energy between crystal and solution). In theory, γeff reflects how the thermodynamic barrier is reduced due to heterogeneous nucleation, that is, if γHOM is the homogeneous interfacial energy, then γeff lies within 0 < γeff < γHOM. The fitted nucleation parameters are listed in Table 1. These parameters can then be used to calculate the nucleation rate J and apparent critical size n* as a function of supersaturation ratio S which are shown in Figure 2c and Figure 2d respectively. Note that the apparent critical size n* is calculated based on an equivalent spherical nucleus with an interfacial energy of γeff. This is an approximation to the actual critical size which theoretically depends on the contact angle between the nucleus and the substrate 24 as well as the shape factor of the nucleus. 32 Moreover, since n* is calculated based on γeff, it also takes into account the effect of heterogeneous nucleation. To explain the observed trends with respect to volume in Figure 2 and Table 1, we consider two relevant phenomena: (1) the interplay between homogeneous and heterogeneous mechanism (2) confinement effects (kinetic and thermodynamic).

Although the median supersaturation 𝑆 ̅ n are essentially identical for the three studied volumes, the effective interfacial energy γeff of the 0.5pL set is significantly lower than that of the rest.

Note that a lower effective thermodynamic barrier is characteristic of the heterogeneous nucleation mechanism. As smaller droplets have higher surface area to volume ratio, the probability of surface nucleation becomes more apparent relative to bulk nucleation. Moreover, lower supersaturation generally favors heterogeneous mechanism while higher supersaturation S generally favors homogeneous nucleation. This is consistent with the plot of nucleation rate J against supersaturation ratio (Figure 2c) where we see that at lower S, the smallest droplet size nucleates faster. However, if the process is governed entirely by the interplay between heterogeneous and homogeneous mechanism, then we would expect that the thermodynamic barrier γeff should obey 0.5 pL < 1.5 pL < 5.5 pL. Interestingly, what we observe for γeff is 0.5 pL < 5.5 pL  1.5 pL. This suggests that other phenomena must be at play. Since the droplet volume is in the picoliter range, kinetic and thermodynamical confinements could have an impact. Note that kinetic confinement stems from the fact that nucleation time scales inversely with the nucleation rate J and system volume V, i.e. tn ∝ 1

𝐽𝑉

. On the other hand, thermodynamic confinement originates from the depletion of the effective supersaturation 24 level during the formation of the pre-critical clusters in finite-sized system. In such finite systems, critical size is determined at a lower "effective" supersaturation, so the critical size would be larger than that of an infinite systems where no depletion occurs. In other words, translated in an infinite system where such confinement through depletion is not taken into account (as in the proposed model herein), it would correspond to a (virtually) higher effective surface energy γeff. However, the exact quantification of such thermodynamic confinement through depletion effects as subcritical population emerge in solution is difficult. This would necessitate, at least, the exact knowledge of the distribution of pre-critical clusters. Nevertheless, we could get some trends from previous works: the higher the solute solubility, the lesser the impact of thermodynamic confinement. For instance, if one compares NaCl and AgCl (known for its low solubility in water) using a "toy model" 24 , we see that the depletion effects in supersaturation needed to form a single critical-cluster has a minor role in the case of NaCl, but a dramatic effect in the case of AgCl (as sketched in Figure S5 of SI). This could explain the slight increase in γeff as the volume decreases from 5.5 pL to 1.5 pL. In other words, the observed trend in γeff could originate from the competition between heterogeneous nucleation and confinement effects, i.e. heterogeneous nucleation tends to decrease γeff while confinement tends to increase γeff. Specific experiments to address this particular problem would need to be made, with a larger range in addressed droplet volumes and comparison with salts/molecules of diverse solubilities.

Regarding the kinetic prefactor, its interpretation is generally more complex as it is related to the mass transfer rate of the monomers to the cluster surface which is a function of the attachment-detachment frequency, viscosity, diffusivity, Zeldovich factor, concentration of nucleation sites etc.

Nevertheless, we remind the reader that the analysis presented here is based on classical nucleation theory (CNT), which has been shown to have inherent limitations particularly at high supersaturations. In such conditions, critical size shrinks to dramatically small sizes where the discontinuity of matter (atoms) should play a role, or at least be the source of large deviations from CNT set of hypothesis (in particular the capillary approximation 33 ). Although nonclassical nucleation theories have been developed 34,35 , their governing equations are more complex and contain multiple fitting parameters. For this reason, we chose CNT to model our experimental data. Moreover, using CNT would allow us to compare our results against the literature (both experiments and simulations). For instance, our fitted γeff ranges from 47.5 to 64.5 mJ/m 2 (Table 1). These values are consistent with that of Hwang et al. 36 who reported γeff = 46.17 mJ/m 2 via electrostatic levitation experiment. Interestingly, a seeded atomistic simulations performed by the group of Peters 10 resulted in an interfacial energy of γ = 47 mJ/m 2 while the molecular dynamics simulation of Bahadur et al 12 yielded a γ of 63 mJ/m. This is a reasonable agreement between theory and experiment 2 . 13 , and seeding simulations (Lamas et al.) 37 .

For further comparison with the literature, we then plotted the nucleation rate J as a function of supersaturation ratio S in Figure 3 together with the experimental and simulated data from multiple research groups. We can see that the magnitude of our measured nucleation rate is very close to that of Gao et. al 15 measured in an efflorescence chamber experiment and Na. et. al 16 who used an electrodynamic levitator trap, a setup that aimed to minimize all possible heterogeneous nucleation sites. Although they reported an interfacial energy between crystal and solution γ = 87 mJ/m 2 , they calculated it from the average induction time while taking A = 10 30 m -3 s -1 as a fixed value (taking induction time as deterministic rather than stochastic). Interestingly, when we used a similar calculation procedure (average induction time and A = 10 30 m -3 s -1 ), we obtained a value of effective interfacial energy γeff of 77, 79, and 80 mJ/m 2 for 0.5 pL, 1.5 pL, and 5.5 pL respectively (while the stochastic approach yields 47.5, 64.5, and 61.9 mJ/m 2 as shown in Table 1 which are closer to the simulation results in the literature).

Thus, the discrepancy in the measured interfacial energy is likely due to two main reasons. First, their approach assumes nucleation as a deterministic process (based on average induction time) while our treatment considers its inherent probabilistic nature (a more realistic view of nucleation). Second, we did not assume any pre-defined value of the pre-exponential factor in the parameter estimation. In the experimental work of Gao et. al 15 where they measured mean efflorescence time, they also fixed the prefactor at a value of 2.8 × 10 38 m -3 s -1 . Furthermore, in the microcapillary experiments of Desarnaud et. al. 4 , they reported J = 0.004 m -3 s -1 at S ≈ 1.6 but they fixed the value of γ at 80 mJ/m 2 . While there are other experimental studies on NaClwater nucleation 31,36,38,39 , values of either nucleation rates or driving forces were not explicitly reported so we are unable to include them in Figure 3. Thus, to the best of our knowledge, our work is the first experimental work that employed a probabilistic approach to measure the interfacial energy between crystal and solution for NaCl-water system without assuming a fixed value of pre-exponential factor. This suggests that the commonly accepted experimental value of A and γ for NaCl crystallization may need to be re-examined. Given that the current theoretical simulations generally overestimate the experimental nucleation rates (Figure 2b), our findings can serve as an additional experimental data for comparison, leading to new insights which could bridge the gap between theory and experiments.

Overall, we highlight that these interesting finite-size effects are clearly observable in our microfluidic experiments which would not be observed in bulk solution experiments. The data treatment of our experiments with CNT model allows us to have a better understanding of nucleation, providing kinetic and thermodynamic information on the system NaCl/water.

IV. CONCLUSION

To summarize, we report a stochastic approach to extract the nucleation kinetic parameters from the induction time distribution of evaporating sessile microdroplets, using NaCl-water as a model system. We showed that by combining a modified Poisson distribution analysis together with an accurate evaporation model, one could obtain reliable nucleation kinetic parameters (both kinetic and thermodynamic). Our results also reveal the competition between the nucleation-enhancing heterogeneous mechanism and the nucleation-inhibiting confinement effects. However, to fully elucidate the underlying mechanisms, it would be interesting to investigate a wider range of droplet sizes together with finer control of evaporation rate, allowing the study of the impact of the rate of change of supersaturation ratio. To investigate quantitatively the impact of thermodynamic confinement, modeling the distributions of the precritical clusters would be essential. The use of non-classical nucleation theories can also be explored.

Given the numerous simulation studies on NaCl nucleation 8, 10-14 , our experimental nucleation parameters based on stochastic approach presented here can serve as additional data points for comparison with theoretical predictions. Moreover, our experimental approach and data-treatment protocol can also be extended to study the nucleation of other salts, biological, and pharmaceutical crystals of interest.

Stéphane Veesler -CNRS, Aix-Marseille University, CINaM (Centre Interdisciplinaire de Nanosciences de Marseille), Campus de Luminy, Case 913, F-13288 Marseille Cedex 09, France; ORCID: https://orcid.org/0000-0001-8362-2531 ; Email: stephane.veesler@cnrs.fr

S4.2 Comparison with constant evaporation rate model

The use of simple linear model has been shown to overestimate the supersaturation ratio at nucleation particularly at the later stages of the evaporation process where the changes in water activity due to the presence of salt becomes significant (demonstrated in our previous work). 4 Here, we show the survival probability plot comparing our evaporation model against that of the simple model. The linear model predicts a maximum S of more than 4 which is highly unlikely for the NaCl-water system. 0.5 pL 1.5 pL 5.5 pL Figure S3. Survival probability plot of supersaturation ratios at nucleation calculated using a system-specific evaporation model (blue) compared against that of a simple linear model (red).

S4.3 Homogeneity of Concentration in Microdroplet

Peclet number Pe is a dimensionless parameter which correlates the relative importance of convective and diffusive transport phenomena.

𝑃𝑒 = 2𝑅𝜅 𝐷 𝑖

where 𝜅 is the evaporation flux (volume loss dV/dt per unit area A), R is the droplet radius and Di is the diffusion coefficient of the solute in the droplet (in this case NaCl ions in water). Generally, Pe << 1 denotes a homogeneous droplet concentration as the diffusion dominates over convection (i.e., no accumulation of solutes in the droplet interface or contact line).

Figure S4.

Peclet number as a function of time.

S4.4. Sensitivity with respect to the empirical parameter of the evaporation model

The evaporation model 4 contains an empirical parameter nx which corresponds to the effective number of neighboring droplets that contribute to the local relative humidity. This parameter nx was adjusted for each droplet size such that the experimental matching time is accurately reproduced.

Here, we show that the resulting distribution of SN (supersaturation at nucleation) is not sensitive to this empirical parameter. Therefore, the estimated CNT parameters are not sensitive to the adjustment of nx. From the microscopy images, we analyzed the standard deviation of the gray-level pixels σ of the region surrounding a microdroplet. The plot of σ vs time allows the extraction of characteristic times tsat (saturation time), tmatch (matching time), and tnuc (nucleation time). From these characteristic times, we can calculate a dimensionless nucleation time 𝜏 = 𝑡 𝑛𝑢𝑐 -𝑡 𝑠𝑎𝑡 𝑡 𝑚𝑎𝑡𝑐ℎ -𝑡 𝑠𝑎𝑡 More details of this procedure (as well as the statistical treatment of outliers) is in Ref. 5 To obtain the supersaturation ratio S at any τ, we employed a tailored evaporation model. 4 The evaporation model contains an empirical parameter nx which can be adjusted such that the theoretical matching time 𝑡 ̃𝑚𝑎𝑡𝑐ℎ reproduces the experimental median matching time 𝑡 ̅ 𝑚𝑎𝑡𝑐ℎ after saturation (results are not sensitive to nx adjustments as shown in S4.2).

The evaporation model then allows the calculation of supersaturation ratio S(t) and volume V(t) as a function of time t. Note that the time variable t in the evaporation model is referenced with respect to tsat (i.e., t = 0 at S = 1). Consequently, the correspondence between the t-scale and τ-scale is given by Note that P(t) is coupled simultaneously with J(t) to find the values of A and γeff that minimize the squared residuals between the experimental and modeled P(t), similar to that of Goh et al 2 . This means that the nucleation rate J is accessible after the fitting. For the curve fitting, we used the least-squares method 'leastsq' as implemented in the LMFIT 6 module.

𝜏 = 𝑡 𝑡 ̃𝑚𝑎𝑡𝑐ℎ

S6. Effect of Solubility on Thermodynamic Confinement

S7. Statistical Analysis of τ-Distributions

S7.1 Reproducibility of Experimental Sets

To assess the reproducibility of the τ-distributions, we analyzed the cumulative probability plots of different experimental sets. Each set corresponds to a group of microdroplets located on the same "line" in the oil bath at a distinct cycle number (i.e. 3 lines with 2 cycles would result in 6 sets). As shown in our previous work, each set can be regarded as "independent" (i.e. negligible patterning and memory effects). 5 We applied a statistical treatment for detecting and eliminating outliers (k-sample Anderson-Darling test) demonstrated in our previous work. 5 The test employs a null hypothesis that the given distributions are statistically identical. Upon eliminating outliers, the results reveal that at 95% confidence level, there is no significant difference between the distribution of each set (i.e. pvalues > 0.05), thus they can be aggregated to form a larger distribution containing 131, 337, and 187 microdroplets for 0.5 pL, 1.5 pL, and 5.5 pL respectively.

S7.2 Comparing τ-Distributions from different volumes

To determine whether the observed difference in τ-distributions obtained from different volumes is statistically significant, we again employed the k-sample Anderson-Darling test (null hypothesis: the given distributions are statistically identical). This resulted in a p-value of less than 0.001 (null hypothesis is rejected). Therefore, the three distributions are statistically different. Moreover, all pairwise comparisons (i.e., 0.5 vs 1.5 pL, 1.5 vs 5.5 pL, 0.5 vs 5.5 pL) also resulted in p-values < 0.001.

Figure S8.

Cumulative probability of aggregated dimensionless nucleation time τ for each microdroplet volume. The p-value < 0.001 suggests a highly significant difference between the three distributions.

Figure 1 (

 1 Figure 1(a) Schematic overview of the experimental setup (adapted from Ref 22 ) (b) illustration of droplet evaporation until nucleation (c) image of a typical 2-D array of microdroplets (d) image of microdroplets at saturation across different sizes at saturation.

Figure 2 .

 2 Figure 2. (a) Cumulative distribution of normalized nucleation Time τ for 0.5 pL (N=131), 1.5pL (N=337), 5.5 pL (N=187) and their corresponding error bands. The median τ is marked by dashed vertical lines and P = 0.5 is marked by grey dotted line. (b) Fitting of the inhomogeneous Poisson distribution (Equation1) with the experimental distribution of supersaturation ratio at nucleation. The median S is marked by dashed vertical lines. (c) nucleation rate (J in m -3 s -1 ) as a function of supersaturation ratio computed using the kinetic parameters in Table1. (d) Apparent critical size (# of NaCl units) as a function of supersaturation ratio computed using the fitted effective interfacial energy γeff.

Figure 3 .

 3 Figure 3. Linearized CNT plot containing the nucleation rate (expressed as log10J) of NaCl in water as a function of supersaturation (1/ln 2 S and S for bottom and top axis respectively). Comparison of our experimental results (dashed blue line is via extrapolation of CNT) to relevant experimental literature data and theoretical simulations. Experiments were based on efflorescence chamber (Gao et. al.)15 , spherical void electrodynamic levitator trap (Na. et. al.)16 , and microcapillaries (Desarnaud et. al.)4 while the simulations were based on seeded atomistic simulations (Zimmerman et. al.)8,10 , forward flux sampling (Jiang et. al.)14 , direct molecular dynamics (Lanaro & Patey)13 , and seeding simulations (Lamas et al.)37 .

Figure S5 .

 S5 Figure S5. Impact of different values of empirical parameter nx on (a) Supersaturation vs time (b) supersaturation ratio SNUC vs normalized time τ (c) probability distribution of SNUC

  This correspondence allows as to calculate S(t) and V(t) at any given τ which are needed in the evaluation of the inhomogeneous Poisson equation coupled with classical nucleation theory as follows:𝑃(𝑡) = 1 -exp [-∫ 𝐽(𝑡)𝑉(𝑡)𝑑𝑡 𝑡 nuc 𝑏 𝑇) 3 ln 2 𝛾 ± 𝑆(𝑡) 𝛾 ± 0 ]

Figure S6 .

 S6 Figure S6. Effective supersaturation upon cluster formation as a function of droplet diameter for model compounds NaCl (high solubility) and AgCl (low solubility). Due to the lower solubility of AgCl, confinement-induced depletion of supersaturation needed to form a single critical-cluster is more dramatic effect in AgCl than NaCl. (Adapted from Ref 7 ).

Figure S7 .

 S7 Figure S7. Cumulative probability of dimensionless nucleation time τ of different experimental sets for each microdroplet volume. The p-values > 0.05 suggest no significant difference between each experimental set within the same volume.

Table 1

 1 Nucleation kinetic parameters and their corresponding standard error obtained from the fit in Figure2b.

	Volume at saturation	Median S at	Effective Interfacial Energy	Kinetic
	Vsat (pL)	nucleation, 𝑆 ̅ n	γeff (mJ/m 2 )	Prefactor
				log10 A(m -3 s -1 )
	0.5	1.89(±3%)	47.5 (±0.4%)	19.8 (±0.2%)
	1.5	1.92(±3%)	64.5 (±0.2%)	24.4 (±0.3%)
	5.5	1.91(±5%)	61.9 (±0.4%)	22.3 (±0.5%)

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

SUPPLEMENTARY MATERIAL

S1. Classical Nucleation Theory for Ionic Systems

Classical Nucleation theory expresses the primary nucleation rate J as the product of the preexponential factor A and an exponential factor containing the free energy cost of forming a critical nucleus ΔG* and thermal energy kbT.

An important difference between the treatment of ionic systems and molecular systems is in the expression of chemical potential difference between solid and liquid. 1 For ionic systems, it is a function of the number of ions forming one formula unit ν (2 for NaCl), and the mean ionic activity coefficient of the solute γ±. These lead to the following expression for ΔG*

with interfacial energy γ between crystal and solution, critical radius Rc, number density of formula units in the solid ρs (2.27 × 10 28 m -3 for NaCl), and S is the supersaturation ratio (c/csat) 1 .

S2. Modified Poisson Distribution Function

In the stochastic view of nucleation, the probability distribution of the nucleation times must be analyzed. In the context of microdroplets, it is normally assumed that the time it takes for a nucleus to grow to detectable size is negligible 2 .

Thus, for constant supersaturation experiments, the cumulative probability of obtaining a droplet with at least one nucleus after time t is a function of nucleation rate J and droplet volume V given as

In the case of evaporating droplet, both the supersaturation and the volume vary with time.

As suggested by Goh et. al., 2 the cumulative probability distribution function becomes

In equation (S4), J(t) can be expressed as a function of supersaturation S(t) by combining with equations (S1) through (S2) . The value of V(t) and S(t) was obtained from an evaporation model (section S4).

S3. Ionic Activity Coefficient vs Supersaturation Ratio

The ionic activity coefficient is necessary to calculate the chemical potential from concentration. To obtain the activity coefficient as a function of supersaturation ratio, we employed the experimental data of Na et al. 3 We then used an empirical function (logistic) to fit the data.

S4. Evaporation Model

S4.1. Evolution of Microdroplet Volume , Supersaturation Ratio, and Evaporate Rates

To determine the supersaturation ratio as a function of nucleation time, we used an evaporation model tailored for our specific system. 4 The evolution of volume and supersaturation with time are plotted as follows. As expected, smaller droplets evaporate faster due to higher surface area to volume ratio. For comparison, the prediction of a simple linear model (constant evaporation rate) is shown (dotted line).