Environmental DNA for detecting Bulinus truncatus: a new environmental surveillance tool for schistosomiasis emergence risk assessment

Under ongoing climate changes, the development of large-scale monitoring tools for assessing the risk of disease emergence constitutes an urging challenge. This is particularly the case for snail-borne diseases such as the urogenital bilharziasis that emerged in Corsica and threat European countries. The expansion of this tropical disease mainly relies on the local presence of competent snail hosts such as Bulinus truncatus . Unfortunately, very little is known about the actual repartition of freshwater snails worldwide which makes new emergences difficult to predict. In this study, we developed two ready-to-use environmental DNA-based methods for assessing the distribution of B. truncatus from water samples collected in the field. We used two approaches, a quantitative PCR (qPCR) and a droplet digital PCR (ddPCR) approach. We successfully detected B. truncatus in natural environments where the snail was previously visually reported. Our environmental DNA diagnostic methods showed a high sensitivity (≈60 DNA copy per mL of filtered water) and a high specificity to B. truncatus . Results obtained in qPCR and ddPCR were very similar. This study dem-onstrates that environmental DNA diagnostics tools enable a sensitive large-scale monitoring of snail-borne diseases hence allowing the delimitation of areas potentially threatened by urogenital schistosomiasis.

Among the emerging diseases documented so far, 17% are induced by pathogens transmitted by vectors and/or intermediate hosts (i.e., vector-borne diseases; VBD) (WHO, 2017). Predicting the emergence and spread of VBD is particularly challenging because it requires knowledge on the distribution and the spatio-temporal dynamics of pathogens and of their intermediate hosts and/or vectors that are also potentially influenced by current global change (Blum & Hotez, 2018;Richards, Lord, Pesko, & Tabachnick, 2009). In this regard, recent studies together with the development of international monitoring networks (e.g., VBORNET in Europe; ECDC, 2011) have made it possible to establish real-time risk maps of emergence and spread of some arthropod-borne-diseases such as malaria, chikungunya, and dengue (ECDC, 2014;Schaffner et al., 2013).
The neglected aspect of most SBD diseases contrast with their significant impact on human health. For instance, more than 250 million peoples are infected by Schistosoma parasites worldwide (Stensgaard, Vounatsou, Sengupta, & Utzinger, 2019). This infectious disease is also responsible for more than 200,000 deaths per year (WHO, 2002). The SBDs concern several vertebrates including humans, livestock, and wild animals. Consequently, these diseases also have a significant economic impact especially in the livestock industry (e.g., 2.5€ billion of annual losses linked to the fascioliasis worldwide; Jones et al., 2018). They are induced by some parasites that depend on a snail host to complete their lifecycle. Our current knowledge on the distribution and dynamics of SBD mainly relies on epidemiological studies based on the prevalence and intensities of adult parasites within their definitive hosts (either humans or animals), while the distribution of snail species involved in the transmission of SBD is generally poorly documented (Kincaid-Smith et al., 2017). This lack of knowledge can be explained by at least four main reasons: first, until very recently, freshwater ecosystems received less ecological attention than terrestrial systems (Collen et al., 2014;Puth & Post, 2005). Second, the distribution of snail populations is highly heterogeneous at both the spatial and temporal scale which makes their monitoring difficult (Lamy et al., 2012).
Third, the available classical methods used to collect and identify snail species are laborious and time-consuming and important sampling efforts are generally necessary for assessing, with certainty, the presence/absence of a targeted species in the field (Lamy et al., 2012). Thus, it is currently inconceivable to monitor snails' populations at large spatial and temporal scales. Finally, once the snails collected from a given site, the taxonomic identification of specimens at the species level requires thorough malacological expertise and generally these identifications still need to be completed by DNA analyses.
Methods based on environmental DNA (eDNA) constitute a promising approach for monitoring biodiversity (Bohmann et al., 2014). These methods were recently applied in a parasitological context to detect (cryptic) parasite species, their hosts and/or vectors (Bass et al., 2015). In this respect, eDNA has been used for detecting targeted parasites species (Huver, Koprivnikar, Johnson, & Whyard, 2015;Richey et al., 2018;Rusch et al., 2018), for characterizing parasite communities (Hartikainen et al., 2016), and to a lower extent for the detection of parasites' vectors (e.g., the invasive mosquito vectors Aedes albopictus; Schneider et al., 2016).
Regarding SBDs, studies have mainly focused on the detection of free-living parasite infectious stages present in the water including Opisthorchis viverrini (Hashizume et al., 2017); Fasciola hepatica (Jones et al., 2018); Schistosoma japonicum (Worrell et al., 2011); and Schistosoma mansoni . This approach targeted toward the detection of parasites is particularly relevant to identify active transmission sites, that is, where the parasite is currently present and susceptible to infect definitive hosts such as humans or livestock species. Surprisingly, very few studies have developed eDNA-based methods to detect the presence of snail species that constitute intermediate hosts for the parasites and that are thus responsible for the transmission of the associated diseases. To the best of our knowledge, eDNA-based methods were developed only for the detection Galba truncatula and Austropeplea tomentosa (Jones et al., 2018;Rathinasamy et al., 2018), both species being involved in the transmission of the liver fluke (Fasciola hepatica), and more recently for the detection of Oncomelania hupensis involved in the transmission of S. japonicum (Calata et al., 2019). Such approach targeted toward the detection of snail hosts is of prime interest to identify sites where the associated disease could emerge (Kincaid-Smith et al., 2017).
In this study, we develop an eDNA-based method to detect the presence of B. truncatus from freshwater water samples. This snail species hosts two flatworm parasites S. bovis, a cattle parasite, and S. haematobium, the latter being responsible for the urogenital form of schistosomiasis. The urogenital schistosomiasis is one of the most severe and prevalent forms of schistosomiasis in Sub-Saharan Africa threatening 436 million individuals and affecting 112 million people throughout the world (WHO, 2018). Moreover, infection with S. haematobium also constitutes an extra risk factor for secondary infections with other infectious agents such as HIV, especially in women (WHO, 2002). The morbidity associated with urogenital schistosomiasis fortunately globally decreased in the last decades although its prevalence has increased throughout Africa (Herricks et al., 2017).
Moreover, the future distribution of this disease in response to climate change is difficult to predict at the continental scale (Blum & Hotez, 2018) and outside the endemic area (Kincaid-Smith et al., 2017). In this regard, an outbreak of urogenital schistosomiasis has recently occurred in Southern Europe in Corsica (Boissier et al., 2016;Holtfreter, Mone, Muller-Stover, Mouahid, & Richter, 2014). This situation has called for an urgent need to develop environmental diagnosis tools to monitor the distribution of compatible freshwater mollusks and hence predict possible new transmission sites in Europe (Kincaid-Smith et al., 2017). Importantly, no vaccine exists against schistosomiasis and mass chemotherapy by the administration of praziquantel, an antihelminth drug, remains the main strategy to control the disease in endemic areas. Today, the control of snail populations under a certain threshold to limit transmission is acknowledged as a complementary strategy necessary to eradicate schistosomiasis (Bergquist et al., 2017). Better characterizing the distribution of host snails is thus fundamental for targeting snail control strategies. So far, eDNA-based methods were developed to detect the presence of free-living stages (cercariae) of S. japonicum (Worrell et al., 2011) and S. mansoni  present in the water and susceptible to infect humans, or within their snail hosts (Amarir et al., 2014), but no methods were developed to detect host snail species except for O. hupensis (Calata et al., 2019).
We here propose a ready-to-use eDNA protocol to detect the presence of B. truncatus established in natural river streams from water samples. We also compared a classical optimized qPCR approach to a protocol using the droplet digital PCR technology (ddPCR), also called "third PCR generation" that provides interesting properties for its application in eDNA .

| Sampling sites
We focused our sampling on five sites along the Cavu and the Solenzara Rivers (Corsica), two geographically isolated rivers where a urogenital schistosomiasis outbreak occurred very recently (Boissier et al., 2016;Ramalli et al., 2018) (Table 1) (Table 1): at the Mulinu Bridge, hereafter denoted "MB"; at the recreational A Tyroliana Park, hereafter named "AT," and at the "3 Pool" site, hereafter denoted "3P" these sites respectively correspond to site 5, 8, and 9 in Boissier et al. (2016).
We sampled water from a fourth site upstream the Cavu River, hereafter denoted "WI" (i.e., Water Intake; corresponding to site 10 in Boissier et al., 2016) where no B. truncatus has ever been observed. This site was sampled to validate or invalidate the absence of B. truncatus at this site so far determined through visual inspection only. Along the Solenzara River, we sampled water from one site (the U Rosumarinu Snack, hereafter denoted "UR") where local infection with schistosomiasis occurred recently (Noel et al., 2017).
Additionally to these sampling sites in Corsica, we also sampled water from an artificial channel in Perpignan city (Southern France) where B. truncatus never occurred. This site was considered as a field negative control. All samples were collected during the first week of September 2018. This period succeeds to the highest density period of B. truncatus in Corsican rivers (i.e., august; Kincaid-Smith et al., 2017) and before high-water period (i.e., October-November).
At each site, we also measured environmental conditions: including water temperature, pH, flow at the water surface (measured in the river bed), and human activity (i.e., presence/absence of anthropic structures/individuals).

| Water sampling protocol
Our sampling protocol was designed to quantify possible effects of filtration volume and sampling location relative to the streambed (i.e., shore vs. streambed of the river) on our ability to detect B. truncatus DNA at a given site. To this aim, a total of 7 water samples were collected at each site. Specifically, three water volumes (1, 3, and 5 L) were filtered at two different locations within the river (i.e., shore vs. streambed) at each sampling site. Moreover, 500 ml of commercial spring water was filtered at each site and rigorously following the same protocol as a technical field negative control. Each water filtration was achieved using a filtration unit including a 0.45-µm PES sterile membrane of 90 mm (VWR PES filter Unit 1 L, Model 514-0301) connected to a manual vacuum pump (Mityvac, Model MV8500). This mesh size constitutes a good compromise between membrane clogging rate and eDNA capture capacity for lotic environment (Hinlo, Gleeson, Lintermans, & Furlan, 2017). The resulting filtration membrane was then removed from its initial filtration unit, immediately stored in a 50 ml falcon filled with a Longmire buffer (Longmire, Maltbie, & Baker, 1997), and stored at room temperature until subsequent DNA extraction. We followed classical specific recommendations to avoid DNA contamination overall the sampling process (Taberlet, Coissac, Hajibabaei, & Rieseberg, 2012). In particular, scissors and forceps used to remove the membranes from filtration units were decontaminated before each sampling by successively placing instruments in a 10% bleach bath during 1 min, 1 min in a DNA AWAY™ solution (Thermo Scientific) followed by a flame sterilization step.
To evaluate the repeatability of our results, technical duplicates were achieved at two sampling sites (i.e., MB and AT) for which we respectively observed the lowest B. truncatus density (≈10 mollusks/ m 2 ) and an intermediate density of ≈50 mollusks/m 2 . A total of 13 filtering membranes were thus obtained for each of these two sites ([3 filtrations volumes × 2 river localities × 2 duplicates] + 1 technical negative control).
Once all the water samples were collected, a visual inspection was achieved to qualify the local density of B. truncatus in the sampling area (low: 10-50 mollusks; medium: 50-100 mollusks; and high: more than 100 mollusks). This was achieved by inspecting rocks and the vegetation on a 1 m 2 surface on the shore of the river and following the same capture effort (30 min per site).

| In silico development of qPCR primers
To design specific primers, we focused on the mitochondrial cytochrome c oxidase I (cox1) gene that have been first proposed for barcoding studies on animal DNA (Hebert, Cywinska, Ball, & deWaard, 2003) and classically (although not exclusively) used for the detection of animal eDNA (Tsuji, Takahara, Doi, Shibata, & Yamanaka, 2019). A total of 25 primer pairs were developed based on an alignment of 210 cox1 sequences from B. truncatus from several countries available from the Barcode of Life Database BOLD (Ratnasingham & Hebert, 2007) to account for intraspecific variability. Primer pairs TA B L E 1 Geographical coordinates (using Google Mercator projection), biotic, and abiotic parameters at each prospected sites species. This bioinformatic pipeline is summarized in Figure S1. All primer pairs showing a 100% specificity to B. truncatus (N = 4) were kept for subsequent qPCR specificity tests.

| Testing the specificity of the qPCR primers
The specificity of the resulting qPCR primers was assessed in vitro by qPCR on DNA templates from a total of 16 freshwater snail species (see Table 2

| Molecular analyses
All molecular experiments were carried out under a sterile table top hood. Before and after each protocol (i.e., DNA extraction, PCR reaction preparation), the working surface was bleached, washed using a DNA AWAY™ solution (Thermo Scientific™), and then exposed to UV light during one hour.
The DNA positive controls were extracted from 5 B. truncatus individuals that were collected in the field in Corsica. Total genomic DNAs from these mollusks were extracted using the E.Z.N.A.® Tissue DNA kit (OMEGA Bio-Tek, Inc) following the "tissue protocol" for animal tissues, supplied in the kit handbook.
Total DNA from each membrane was extracted using the DNeasy® Blood & Tissue Kit (QIAGEN) which is commonly used in eDNA studies and which is supposed to yield an important amount of extracted DNA (Goldberg et al., 2016;Hinlo et al., 2017). We followed an adapted version of the "tissue" protocol from the kit handbook. Briefly, membranes were subdivided into 4 equal sections to allow the use of microcentrifuge tubes during subsequent extraction steps. Lysis was performed at 65°C during 1 hr using a final solution containing 567 µl of ATL buffers and 63 µl of proteinase K. We next added 630 µl of AL buffers and 630 µl of 100% ethanol in the following steps (steps 2 and 3). The subsequent steps were performed following the manufacturer protocol, and the resulting genomic DNA was finally eluted in 100 µl of AE buffer. At the end of the extraction, we thus obtained four eluates of 100 µl for each membrane.

| Quantitative PCR approach
The four sections from each membrane were analyzed twice meaning that a total of 8 qPCR reactions were ran for each membrane. Finally, we considered that eDNA from B. truncatus was present on a filtering membrane if two qPCR replicates from the same membrane section (i.e., at least 1 section fully positive) displayed a melting curve specific to B. truncatus.
Finally, the eDNA signal of B. truncatus obtained in qPCR was validated on a subset of 20 membranes section. Sequences from these qPCR products were obtained using ABI 3730xl sequencer at the GenoScreen platform (Montpellier, France).

| Droplet digital PCR approach
Because ddPCR is highly sensitive and because it is supposed to be robust even in the potential presence of PCR inhibitors, the ddPCR reactions were ran on pools of DNA extracts from the four sections of each membrane. Noteworthy, this pooling strategy was also tested in our qPCR analyses but failed to provide repeatable results. The ddPCRs were ran using the TaqMan chemistry on a QX200 AutoDG

| Data analysis
Several factors can influence the concentrations of eDNA isolated from environmental samples including the sampling protocols used

| Sampling sites
The environmental conditions (e.g., temperature, pH) were very similar among the sampling sites. Likewise, all sites presented human activities (generally bathing) except the control site that is inaccessible to the public. The water flow varied from one site to another, ranging from 0.1 m/s for the site WI to 1 m/s for the site MB; these data are displayed in Table 1.

| Development of qPCR primers
Among all the primer pairs tested, Btco1F (5′ TYGAAGGAG GGGTTGGAACA 3′) coupled with Btco1R (5′ RKTRATTCC TGGTGCYCGT 3′) provided the best results regarding the qPCR specificity tests (i.e., only the DNA from B. truncatus or pooled DNA containing B. truncatus DNA were positive in qPCR) and were thus used for subsequent qPCR analyses. This primer pair resulted in the amplification of a 179 bp amplicon specific to B. truncatus.

| Quantitative PCR approach
(i.e., University channel) or in the 6 field negative controls (i.e., spring water filtered at each sampling site). Noteworthy, two among the total of 48 field negative control replicates showed a measurable C t value (Site UR: 37.39 and Site MB: 39.17). However, these two samples displayed unspecific melting curves compared to those obtained from positive controls. Likewise, none of our 6 qPCR negative controls showed a measurable C t value. These results indicate that DNA contamination during the overall process from water sampling to qPCR reactions is very unlikely.
The water samples filtered on the river shore and in the streambed provided similar results, that is, 85.7% [18/21] of positive samples at the shore and 81% [17/21] at the streambed location ( Figure 1).
However, we observed a slight tendency in the C t values to be higher for samples filtered in the river shore (35.67; SD ± 1.79 C t ) than those filtered in the river streambed (37.71; SD ± 1.69 C t ).

| Droplet digital PCR approach
The ddPCR provided the same results to those obtained using the qPCR approach. DNA extracts from the filtered membranes ob-

| Statistical analyses
Among all GLMs performed using the qPCR dataset, the best-supported model was the following: C t ~ volume + position + density (Table 3). According to this model, the snail density (p < .001) and the volume of filtered water (p < .01; see Table 3 and Figure 3) significantly affect the C t values obtained by qPCR after molecular process.
Accordingly, a higher water volume and a higher mollusk density decrease the C t value obtained which means that the concentration of B. truncatus eDNA is higher for these samples. A significant positive effect of the sampling position (shore) was also found (p < .05; Table 3), hence indicating that the C t value obtained by qPCR tends to be higher and thus eDNA concentration is lower in water samples collected on the river shore compared to the streambed.
Concerning the ddPCR dataset, the best-supported model explaining c/L predicted a single positive effect of the B. truncatus densities (p < .001; Table 3). This result indicates that a higher mollusk density results in a higher eDNA concentration which is congruent with our expectations (Figure 4).

| D ISCUSS I ON
Under the current global change context, it is essential to develop large-scale monitoring tools to understand and predict the dynamics and epidemiology of emerging diseases (Ellwanger et al., 2019).
In this study, we present the first application of two eDNA-based methods, namely a ddPCR-and a qPCR-based protocols, to de- when using 3 or 5 L of water during filtration but dropped to 91.7% when filtrating 1 L of water using both qPCR and ddPCR. This result indicates that according to our results, at least 3 L of water is necessary to reach an optimal sensitivity whatever the density of mollusks at the sampling site.
More specifically, the sampling procedure that provides the best detection rate consists of a 5-liter filtration preferentially in the river streambed. Whatever the PCR technology used, we can argue that in this current version our protocols designed for water collection, sample preservation, and eDNA extraction are robust. The considerations of these three protocol aspects are critical for establishing a reliable eDNA-based approach (Goldberg et al., 2016) but are generally neglected in most publications. Only 5% of eDNA-based studies provided critical aspects of their sampling protocols (Dickie et al., 2018) and most eDNA studies suffer from a lack of reproducibility especially when their protocols are applied in natural environments for the first time (Tsuji et al., 2019). Finally, no standardized method for eDNA sampling is available so far (Tsuji et al., 2019). This technical lack raises the necessity of further studies benchmarking all existing sampling protocols in an attempt to define standardized sampling pipelines adapted to different environments and environmental samples (i.e., water, sediment, and air).
In this line, several aspects could be developed to improve our eDNA-based monitoring tools (Dickie et al., 2018). First, the robustness of our protocols should be validated in different freshwater environments with different hydrodynamic profiles. In the present study, all of the prospected sites, but the negative University channel site, presented a lotic hydrodynamism (i.e., flowing water) and whether the sampling strategy developed in this study could be applied in lentic environments (i.e., still freshwater) still needs to be assessed. Lentic environments are known to present some specific characteristics that can influence the eDNA decay (e.g., temperature, sunlight, and pH; Klobucar, Torrey, & Phaedra, 2017;Strickler, Fremier, & Goldberg, 2015). They also generally contain significantly more suspended solids (Takahara, Minamoto, Yamanaka, Doi, & Kawabata, 2012), which could increase the PCR inhibitors concentration in eDNA samples. Second, the ergonomic aspect and efficiency of our sampling tools could be further optimized. For instance, filtration units are voluminous and are probably F I G U R E 1 Membranes positiveness obtained from the quantitative PCR approach. The displayed percentage is calculated on the 8 qPCR replicates from the samples filtered on the river shore (a) and at the river streambed (b inappropriate when field conditions are difficult (e.g., isolated prospecting sites; long fieldwork). The size-reduction of filtration devices could also be contemplated, for instance using smaller membrane holders coupled with a vacuum automated pump, although the latter could also be a limitation in some field conditions. Likewise, the use of silicate micro-beads (Wilcox et al., 2018) (Hinlo et al., 2017), these findings could also constitute relevant improvement to consider.

TA B L E 3
Best fitted GLM models that explain (Ct)s and (c/L)s obtained by qPCR and ddPCR based on the AIC criterion and number of parameters (nPars) (a) and the parameter estimates obtained from the chosen models Volume of filtered water (L) and Molluscs densities (individuals/m²) Residuals for cycle threshold 10 molluscs/m² *** 50 molluscs/m² *** 100 molluscs/m² *** Another major finding of this study is the potential of the ddPCR as an alternative quantification platform to the qPCR for environmental samples. In our study, qPCR and ddPCR gave strictly the same results regarding sensitivity since the same samples were found positive whatever the technology used. The only difference between these two PCR technologies is that ddPCR was able to provide reliable results on the pooled samples while qPCR did not. This result hence provides a small advantage of ddPCR compared to qPCR because it requires less replicates to give reliable results. Finally, because we did not conduct comparative analyses using qPCR and ddPCR on the same dilution standards, it is difficult to conclude on the most accurate technology regarding absolute quantification of environmental DNA and ultimately mollusk density in the field. Thus, based on our results and because ddPCR is generally argued as very sensitive and accurate even for the amplification of low quality DNA samples (e.g., high amount of inhibitors and/or untargeted DNA; Hindson et al., 2011;Taylor et al., 2017), we predict that the ddPCR will be increasingly used for detecting eDNA in the next few years.
The financial aspect of these protocols is also an important consideration for the development of large-scale monitoring tools. To this end, we investigated and compared the prices and time cost between the three approaches for monitoring B. truncatus presented in this paper (i.e., qPCR, ddPCR and visual prospecting), following the same method as in . The estimated prices for each site included the reusable materials (e.g., forceps, DNA primers, and probe), meaning that these prices are applicable for a first-time application of the protocol but should be devaluated in subsequent applications (Table S1). As expected, traditional visual survey remains the most economic approach (i.e., ≈20 EUR for all prospected sites; 12 hr/8 sites; Table S1). However, this approach is not adapted for large-scale monitoring as they require extensive workforce proportional to the number of prospected sites and huge malacological knowledge. Regarding the two eDNA-based methods developed in this study, we estimated that the price between the two PCR technologies differed by 5 EUR (i.e., ≈60 EUR/sites; 6 hr/5 L × 8 sites for qPCR and ≈65 EUR/sites; 5h20/5 L × 8 sites for ddPCR; Table S1).
Importantly, however, the estimation for ddPCR included service delivery. When estimating the price of ddPCR if realized in our laboratory facilities, prices are expected to be similar or even less expensive than qPCR analysis in a recurrent use framework. Expenditures related to ddPCR apparatus acquisition, which is much more expensive than those for qPCR, could be the only limiting aspect of this technology . However, this expense is compensated by the ability to obtain an absolute quantification from a single reaction compared with qPCR. Hence as denoted in the previous studies, ddPCR technology would be an interesting alternative to qPCR in quantifying eDNA from environmental samples as this enables absolute quantification from a single sample and due to the inhibitors dilutions during the sample singulation step .
In a schistosomiasis prevention perspective, these two eD- ). In this context, we believe that these developed tools will help better assessing the risk of schistosomiasis emergence worldwide.
Importantly however, other factors need to be accounted to establish efficient risk maps of emergence and spread of schistosomiasis at large scale, which include human population movements, F I G U R E 4 Representation of the estimated B. truncatus DNA copy per liters of filtered water estimated using ddPCR as a function of B. truncatus densities at the sampling sites. *** represent differences with significant levels below 0.001 10*** 50*** 100*** −500 0 500 1,000 1,500 2 ,000 Molluscs densities (individuals/m²) Residuals for DNA copy per litre the spatial fluctuation of schistosomes' intermediate hosts under current climate changes, the invasion ability of African schistosomes species, the effect of hybridization between schistosome species on the invasion ability of the parasite, and finally the susceptibility of other European freshwater snails' species to Schistosoma parasites (Kincaid-Smith et al., 2017). To this end, the combination of eDNA-based tools with spatial modeling approaches such as ecological niche modeling (ENM) would improve our ability to generate "real-time" risk maps of schistosomiasis in Europe and endemic countries (Kincaid-Smith et al., 2017). ENM approaches were already used for spatial monitoring of freshwater snail species involved in SBD transmission such as fascioliasis and schistosomiasis (Cordellier & Pfenninger, 2009;Pedersen et al., 2014). Therefore, the use of ENM modeling approaches coupled with eDNA-diagnosis tools appears to be a next development step for SBD risk assessment endeavors. the ddPCR experiments and subsequent first step analysis. We also wish to thank Sebastien Calvignac-Spencer and another anonymous referee for their help in the revision of previous versions of the manuscript. SM is also greatful to Didier Mulero for his help during field work and also the Occitania region for its involvement in this project.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
The dataset supporting the conclusions of this article is included within the article.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the