Conservation Challenges Imposed by Evolutionary History and Habitat Suitability Shifts of Endangered Freshwater Mussels under a Global Climate Change Scenario

Abstract

Climate change and associated shifts in temperature and precipitation patterns have become an increasing concern as drivers of ongoing biodiversity loss. The Mediterranean region is particularly vulnerable, being both a biodiversity hotspot and a region very prone to desertification. Freshwater mussels are amongst the most threatened invertebrate taxa worldwide. Unio tumidiformis is an endemic and endangered species restricted to the southern Iberian Peninsula, living in temporary Mediterranean-type streams. Freshwater mussels need a fish host for successful larval transformation, meaning U. tumidiformis must belong to the genus Squalius. The main objective of this study was to evaluate the vulnerability of U. tumidiformis to climate change, by studying its population genetics and evolutionary history, its current and future habitat suitability, and that of its hosts. Genetic population structure and diversity were assessed using Single-Nucleotide Polymorphisms through Genotyping by Sequencing and used to infer species evolutionary history. The species potential distribution was modeled using an ensemble forecasting approach, and future shifts in habitat suitability were assessed with the projected climate data layers from Worldclim. Most populations showed extreme genetic differentiation (Fst up to 0.745), even from close neighboring ones. Upper Guadiana populations were more diverse and less differentiated. We hypothesize that U. tumidiformis originated in Upper Guadiana and followed the same colonization routes as their hosts with numerous founder effects and bottlenecks. Our results also predicted a reduction of 99% of climatically suitable areas for U. tumidiformis in the Iberian Peninsula until 2040. For the fish hosts, a maximum 42% reduction in suitable areas was estimated throughout the century, with remaining adequate habitats in the north. Our results suggest that difficult conservation options are necessary, prioritizing the preservation of populations, translocations to the northern area of its historical range and stream engineering to increase resilience to droughts.

1. Introduction

Climate change is now recognized as a major driver of biodiversity loss on the planet [1,2] and has led to widespread local extinctions [3]. In response to climate change and expected shifts in temperature and precipitation patterns, some species may shift their distributions to more poleward locations, a pattern enhanced in northern latitudes [4]. However, for some aquatic species with reduced dispersal mechanisms or living in low-connectivity environments, such as freshwater bivalves, these mechanisms to cope with climate change might not be readily available [5].

The Mediterranean region is prone to water scarcity and overall increased temperature due to ongoing climate change [6]. In Mediterranean-type streams, flow conditions show high natural seasonal variability, being reduced to a few pools in summer and flooding in winter [7]. The frequency, timing, and intensity of these seasonal events vary greatly over the years [8]. Climate change is expected to lead to an increase in drought frequency and intensity in this area, which might lead to pools drying out, further increases in temperature, and overall decreases in oxygen levels [9]. This may exacerbate already difficult conditions for freshwater organisms, particularly imperiled invertebrates and fish, leading to local extirpations [10].

Species distribution models (SDMs) that consider several types of variables (e.g., bioclimatic, topographic, anthropogenic) are important tools for predicting the geographic distribution of species. In recent years, SDMs have been used to predict the impacts of climate change on species distribution, with multiple future scenarios being considered [2]. By providing information on habitat suitability, SDMs are important tools in ecological, biological, and resource management, especially in the context of global climate change [11], and are essential to improve the definition of conservation strategies [12]. However, SDMs alone cannot predict the impact of distribution changes and shifts in the population structure of species and the corresponding impacts on their viability. In fact, the genetic population structure is recognized as an important element for evaluating the extinction risk of individual populations, particularly if climate change is a threat [1,4]. Traditionally, higher genetic diversity is associated with lower extinction risk [13,14]. However, this general assumption is not always applicable, and it has been shown that higher genetic diversity in small populations may be associated with higher extinction risk [14], even though small and fragmented populations are especially prone to be affected by genetic drift and inbreeding depression [15]. Climate change-induced changes in species range can impact the genetic structure itself [4], so the two factors are interconnected and should be considered together.

The genetic population structure can be accessed using different markers, with Single-Nucleotide Polymorphisms (SNPs) becoming increasingly used to evaluate evolutionary and conservation issues in non-model organisms [16]. The use of SNPs presents several advantages when used to study population genetic variations, including a more extensive screening of the genome if compared to the use of microsatellites [17]. Even though microsatellites may be more polymorphic, SNPs seem to lead to more reliable inferences of patterns of genetic structure and diversity in different aquatic organisms [18]. Additionally, using more molecular markers based on the whole-genome SNPs of a few individuals, rather than many individuals at few markers, allows us to overcome the restrictions imposed by very low population densities where only a few specimens are available. Therefore, we used SNPs as genetic markers in this study as they are an efficient genomic tool increasingly used to study evolutionary patterns and make informed decisions regarding conservation priorities for endangered species [19].

Freshwater mussels are amongst the most threatened invertebrate taxa worldwide [19], currently with ten extant species known in the Iberian Peninsula. Unio tumidiformis Castro, 1885, is a narrow endemic with a distribution range restricted to the southwestern Iberian Peninsula [20], living in Mediterranean-type streams mainly in the Guadiana basin, with isolated populations in the Mira and Sado basins in Portugal and the Guadalquivir basin in Spain [7]. There are also historic records further north, dating from the late XIX century and early XX century, particularly from the Tagus basin [20]. The genetic structure and diversity within the extant range of U. tumidiformis and its respective genetic diversity have not yet been assessed, and there is little information about the demography of its populations [7]. This species is endangered and classified as Vulnerable by the IUCN red list as well as protected under the Habitats Directive (Annex II and IV) as it was previously identified in the Iberian Peninsula as U. crassus [7,21]. The main recognized threat for U. tumidiformis is the modification of the hydrological regime of rivers and streams caused by overexploitation of water resources coupled with climate change impacts [22,23]. The destruction of riverine habitats (with consequent increases in water temperature and evaporation due to reduced shadow) adds to this threat, while the introduction of invasive species is also a known and growing problem.

Almost all freshwater mussels (order Unionida) need a host fish to successfully reproduce [24]. In the case of U. tumidiformis, it can only metamorphose successfully on fish that belong to the genus Squalius [25]. Within the mussel’s native range, three species have been identified as suitable hosts: Squalius alburnoides, S. pyrenaicus, and S. torgalensis. However, two geographically neighboring although allopatric species, S. carolitertii and S. aradensis, have been proven successful hosts in a laboratory setting [25]. The dispersal capacity of freshwater mussels is directly related to their hosts [26], which in turn determines the gene flow within and between populations and consequently their genetic population structure [27]. The dispersal using fish hosts also determines the potential displacement of mussels to more adequate habitats under climate change shifts, which in the case of riverine habitats implies that movement is limited to the dendritic nature of the watersheds [28]. In fact, aquatic organisms such as freshwater fish are usually unable to disperse across river basins, effectively preventing potential range shifts promoted by climate change to more adequate habitats [29]. For these reasons, fish hosts are essential components to consider in conservation plans for freshwater mussels [30,31], especially when assessing climate change impacts on future species distribution [32] and on patterns of genetic population structure [19,27,33].

The main objective of our study was to evaluate the vulnerability of U. tumidiformis to climate change. Specific objectives included (1) evaluating the genetic structure and diversity of the species within its range by using genome-wide markers, inferring its evolutionary history; (2) establishing the current and future habitat suitability of mussels and their hosts based on climate change projections; (3) recommending conservation actions to minimize the species extinction risk.

2. Methods

2.1. Sampling

A total of 60 individuals of U. tumidiformis were collected from 14 locations belonging to ten hydrological sub-basins and three hydrological basins (Figure 1 and Table S1). Specimens were visually located by informal sampling (as defined by [34]) in rivers with previously known populations, using batiscopes, snorkeling, or scuba-diving. For genomic analyses, a small piece of tissue from each individual was collected in situ and preserved in 95% ethanol until DNA extraction. Sampled specimens were then placed back in the collecting site. For some populations, previously collected tissues preserved at −80 °C at the Museo Nacional de Ciencias Naturales in Madrid were used (Table S1). Up to 5 specimens were used from each location, as rather than sampling many individuals with few molecular markers, we decided to sample more markers in fewer individuals. More markers provide more information about the evolutionary history, and this was the only adequate sampling strategy considering only a few mussels were located from most populations.

2.2. DNA Extraction and Genotyping by Sequencing

Total genomic DNA was extracted for each individual from the preserved tissue using the commercial kit DNeasy Blood and Tissue (QIAGEN), following the manufacturer’s protocol. The final DNA concentration of all samples was determined using the Qubit 2.0 Fluorometer, and DNA fragmentation was verified in a 1% agarose gel. To obtain genome-wide Single-Nucleotide Polymorphisms (SNPs), a paired-end Genotyping by Sequencing (GBS) approach (adapted from [35]) was used. This was performed at LGC Genomics, GmbH, Germany. Upon arrival, DNA was fragmented using the MslI restriction enzyme, and libraries were constructed and sequenced using Illumina NextSeq with a read length of 150 base pairs (bp).

Sequence quality was assessed using FastQC 0.11 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/, accessed on 1 March 2020) and MultiQC [36]. All sequences passed quality standards to be included in the subsequent analyses. The Process_Radtags program of Stacks v.2.5 [37] was used to truncate all reads to the same length (135 bp) and to discard reads with low quality scores and uncalled bases, using the default settings for the window size (0.15x read length) and the base quality threshold (10 Phred score). The Stacks software only cuts the final end of the sequence reads; therefore, Trimmomatic v.0.36 [38] was used to eliminate the first 5 bases in all sequences. The final length of reads was 130 bp. The de novo pipeline from the Stacks v.2.5 software [37] was used to build a catalog of loci. A subset of individuals covering all sampling sub-basins was used in the construction of the catalog (Table S1). In the de novo Stacks pipeline, first, the identical reads are grouped in stacks and these stacks are merged into loci in each sampled individual. Then, a catalog is created by determining which loci are homologous across all the analyzed samples. We followed the recommendation of [39] to determine the optimal parameters for the construction of the catalog by testing different values of m (minimum number of reads required to form a stack), M (number of mismatches allowed between loci when processing a single individual), and n (number of mismatches allowed between loci during the construction of the catalog). We selected M4m4n5 for the construction of the catalog. Given the possibility that forward and reverse sequences of the same DNA fragment were treated as different loci, similar reads within the catalog were clustered using CD-HIT-EST from CD-HIT v.4.8.1 [40] with a word length of 9 and a sequence identity threshold of 0.98. The paired pre-processed sequences (those truncated to 130 bp) of all individuals were aligned to this catalog using BWA-MEM from BWA v.0.7.17 [41] with default parameters. We sorted the output alignments and removed unmapped reads using Samtools v1.10 [42]. We used Freebayes v.1.3.1 [43] to carry out the SNP calling based on haplotypes for variant detection and considering the cleaned catalog as a reference to genotype SNPs, discarding reads and bases with low quality, and without using Hardy–Weinberg equilibrium priors (-p 2-min-mapping quality 30-min-base-quality, 20-hwe-priors-off). We filtered the obtained vcf file to eliminate indels, sites with an excess of heterozygosity, sites with more than 50% missing data, and alleles with a minimum allele frequency below 5%, using a combination of options from VCFtools v0.1.15 [44] and BCFtools v1.6 [42]. After using these filters, individuals with more than 50% missing data were removed. The final vcf file contained 72736 SNPs, and 48 individuals from the initial 60 remained. In this process, we eliminated the sole specimen from the river Terges e Cobres due to its high amount of missing data. Raw data sequencing has been deposited in GenBank (Bioproject: PRJNA1087543; Biosample accession numbers: SAMN40447825-SAMN40447884; SRA accession numbers: SRR28340993).

2.3. Genetic Diversity and Genetic Structure Analyses

Genetic diversity and structure were investigated at different hierarchical levels (sampling location and sub-basin). Genetic diversity was evaluated by analyzing the expected and observed heterozygosity using custom scripts. The program populations from Stacks v.2.5 [37] were used to estimate different genetic diversity statistics such as number of private alleles, percentage of polymorphic sites and nucleotide diversity, and FIS values.

FST values and their p-values after Bonferroni correction [45] were also estimated through Weir and Cockerham’s estimator [46] in Arlequin v.3.5 [47]. In Arlequin v.3.5, population structure was also assessed with an Analysis of Molecular Variance (AMOVA) at different hierarchical levels (sampling locations and sub-basins). Significance was tested with a 10,000-permutation test. Genetic structure was further analyzed using sNMF v.2.0 from the R package LEA [48], testing values of ancestral populations (K) from 1 to 10 (number of sub-basins), and cross-entropy was visually inspected to determine the most plausible value for K.

Within the Guadiana river basin, approximate geographic distances between sampling locations were measured by hand (in km), following the natural meandering of rivers, using the Path tool available in Google Earth Pro v.7.3.3.7786. Insurmountable barriers for fish (dams, weirs, and waterfalls) were mapped along with the drawing of pathways. To evaluate the possible existence of isolation by distance, linear regression analyses were conducted using approximate geographic distances between sampling locations (log transformed) and their respective pairwise FST/(1-FST) estimates (using pairwise FST values at sampling location level—Table S3).

2.4. Ecological Niche Modeling and Habitat Suitability

Georeferenced presence data for each modeled species were obtained from different sources: Global Biodiversity Information Facility (GBIF), Portuguese Nature Conservation Institute (ICNF), and sampling records from the Invertebrate Red List from Portugal [22]. All novel presence data were made available on GBIF. A database was produced with the following number of locations: U. tumidiformis—78 and Squalius (combining records of S. alburnoides—890, S. carolitertii—710, and S. pyrenaicus—1289)—2889. We included S. carolitertii even though it is not sympatric with U. tumidiformis because it is ecologically very close to S. pyrenaicus [49] and proved to be equally suitable as a host under a laboratory setting [25]. The suitability was modeled for the whole Iberian Peninsula, as the original distribution of U. tumidiformis is assumed to be greater than the current.

To describe the environmental conditions that might influence habitat suitability, eight bioclimatic variables were selected from nineteen available on Worldclim (www.worldclim.org): BIO1 = Annual Mean Temperature; BIO5 = Max Temperature of Warmest Month; BIO6 = Min Temperature of Coldest Month; BIO8 = Mean Temperature of Wettest Quarter; BIO9 = Mean Temperature of Driest Quarter; BIO12 = Annual Precipitation; BIO13 = Precipitation of Wettest Month; BIO14 = Precipitation of Driest Month. Those were considered to be ecologically relevant for the species, representing extreme temperature or precipitation and general climatic trend variables. Highly correlated variables (Pearson’s correlation > 0.80) were excluded to avoid collinearity [50].

The projected climate data layers include four different time frames, from 2021 to 2100, with 20-year intervals and five shared socioeconomic pathways (SSPs). Global climate models (GCMs) used were BCC-CSM2-MR, CNRM-CM6-1, CNRM-ESM2-1, CanESM5, IPSL-CM6A-LR, MIROC-ES2L, MIROC6, and MRI-ESM2-00, averaged into one raster for each bioclimatic variable using Arcgis Toolbox Spatial Analyst tools. We decided to select one moderate scenario for modeling—SSP245—and a more extreme one—SSP585. Both are considered “Tier 1” scenarios and the top priority for modeling organizations (scenarios to be modeled first). Both mussels and hosts were modeled for the following combinations: SSP245 for all the periods available: 2021 to 2100; SSP585 for the first period: 2021–2040.

To model the potential distribution, an ensemble forecasting approach was used to create a “mega/meta model” [51]. Modeling was performed in R version 2.14.0 (R Core Team, 2011) using algorithms in the BIOMOD2 package [52]. These algorithms included regression algorithms: GLM (generalized linear model) and GAM (generalized additive model); classification methods: CTA (classification tree analysis) and FDA (flexible discriminant analysis); machine learning methods: ANN (artificial neural network), RF (random forest for classification and regression), and GBM (generalized boosted regression model); a climatic envelope method: SRE (surface range envelope). Modeling outputs were combined in an “ensemble” [53] using the median as a measure of central tendency. Presence records were coupled with randomly generated pseudo-absence data to increase model robustness and avoid bias towards more prevailing responses [51], preventing underestimations due to the lack of species absence data. Models were built with 80% of the data, retaining the remaining 20% for evaluating predictions [52]. Model performance was evaluated using true skill statistics (TSS) as it takes into account omission and commission errors, ranging from −1 to 1, and is not affected by prevalence [54]. For TSS, values ranging from 0.2 to 0.5 are considered poor, from 0.6 to 0.8 are considered useful, and values larger than 0.8 are considered good to excellent (as in [55]). We assumed that the most important variables contributing to the model would be those with a relative importance above the mean of the predictor variables [51]. Ensemble models were produced using a weighted approach based on TSS values, and to reclassify the resulting continuous maps into binary maps, the sensitivity–specificity equality approach (http://r-forge.r-project.org/projects/biomod/, accessed on 6 March 2023) was used, by minimizing the absolute value of the difference between sensitivity and specificity [56].

3. Results

3.1. Genetic Diversity and Genetic Structure

Genetic diversity was assessed using genome-wide SNPs and measured as average observed (Ho) and expected (He) heterozygosity, ranging from 0.073 to 0.229 (Table S2). The highest values of Ho, which reflect the proportion of heterozygotic sites per individual for all the positions analyzed in their genomes, were found in rivers from the Upper Guadiana basin (Bullaque and its tributary Milagro, respectively, 0.185 and 0.155) and the lowest in the Torgal and Marateca sub-basins (0.073 and 0.083), belonging to the Mira and Sado basins, respectively, the two westernmost basins of the distribution area of the species. Regarding the expected heterozygosity (He), which depends on allele frequencies and reflects the expected frequency of heterozygotes in a population after one generation of random mating, the values generally followed the same geographic pattern: lower diversity in the westernmost populations of Sado and Mira river basins and higher in the Guadiana (Table S2). Amongst the Guadiana sub-basins, a decreasing gradient of He values was observed from the Upper to the Lower Guadiana sub-basins (Table S2). At both sampling location (population) and sub-basin levels, the highest number of private sites was found in rivers or sub-basins from the Lower Guadiana (Table S2). However, the highest level of polymorphism was found in Upper Guadiana rivers and sub-basins, while the lowest percentage of polymorphism, although with an intermediate number of private sites, was found in Torgal. All populations showed negative FIS values, except for the Bullaque and Estena sub-basins, suggesting a generalized trend of excess heterozygosity for all populations.

Genetic structure was clearly marked, and pairwise F_ST comparisons based on Weir and Cockerham’s estimator revealed generally high values between the basin and sub-basin comparisons. At the basin level, F_ST was higher between the Sado and Torgal basins (0.726) and ranged from 0.451 to 0.487 between the Guadiana and the other basins. F_ST between the Upper and Lower Guadiana groups of sub-basins was 0.211.

At sub-basin and population levels (

Table 1. Genetic differentiation (Fst pairwise comparisons) for 10 sub-basins based on Weir and Cockerham’s estimator.

	Torgal	Sado	Marateca	Ruidera	Bullaque	Estena	Ardila	Chança	Vascão	Odeleite
Torgal	-	0.630	0.743 *	0.563	0.364 *	0.574	0.565	0.597	0.664 *	0.628
Sado		-	0.252	0.474	0.329	0.450	0.325	0.548	0.594	0.561
Marateca			-	0.643 *	0.500 *	0.628 *	0.647	0.686 *	0.703 *	0.703 *
Ruidera				-	0.214 *	0.383	0.300	0.479	0.532 *	0.522
Bullaque					-	0.140	0.072	0.140	0.355 *	0.393 *
Estena						-	0.297	0.441	0.507 *	0.476
Ardila							-	0.348	0.496	0.402
Chança								-	0.567 *	0.558
Vascão									-	0.581 *
Odeleite										-

* value is significant after Bonferroni’s correction (p-value = 0.001).

and Table S3), the highest F_ST values were found in the comparisons between the Marateca (sub-basin of the Sado) and most of the other populations; nevertheless, values between Marateca and populations from the sub-basins of the Upper Guadiana were lower (0.50–0.64), although geographically more distant, than those between the Marateca sub-basin and the closer populations from the Torgal and Lower Guadiana sub-basins (0.68–0.70) (Table 1). In the same way, F_ST values between Torgal and Upper Guadiana were lower than those between Torgal and other locations geographically closer to this basin.

Although the AMOVA analyses showed the highest percentage of explained genetic variation (71.51%) within individuals, this variation was negligible within locations (Table S4), indicating almost no variation at each site. In contrast, 10.02% of the genetic variation could be explained among locations, and as much as 43.98% among sub-basins, showing strong differentiation between sites and sub-basins, which is consistent with the high F_ST values obtained.

SNMF analysis results supported a strong structure as well, with K = 6 being the most probable number of ancestral populations based on the cross-validation (Figure S1). The main genetic clusters obtained were Torgal, Sado, Upper Guadiana (including the São Pedro population), and three differentiated rivers within the Lower Guadiana (Barranco do Vidigão, Odeleite, and Vascão) (Figure 2). Most populations showed an overwhelming dominance of one ancestry, except for the São Pedro stream (Lower Guadiana) and Estena (Upper Guadiana), which presented a much more balanced contribution of the six ancestries. Overall, F_ST and SNMF analyses show that the São Pedro population seems to be genetically closer to the Upper Guadiana group than to geographically closer populations within the Lower Guadiana group of sub-basins.

3.2. Association between Genetic Differentiation and Geographic Distances

When looking at the association of F_ST values with distances between locations, it becomes evident that in the Guadiana basin, larger geographic distances do not necessarily correspond to higher genetic differentiation (Figure 3). The lowest geographic distances, as expected, do match the lowest F_ST values corresponding to comparisons between sites within the same sub-basin (Bullaque and Vascão, purple in Figure 3). The highest F_ST values, however, correspond to average geographic distances, found between populations within the Lower Guadiana (orange in Figure 3). In contrast, geographic distances between populations within the Upper Guadiana group are similar but correspond to much lower genetic distances (yellow in Figure 3). The highest geographic distances are found between the Upper and Lower Guadiana groups and show F_ST values intermediate between the two previous ones (blue in Figure 3). Finally, the São Pedro population again shows the closest genetic affinity to the Upper Guadiana group, despite being geographically closer to the Lower Guadiana group of sub-basins (green in Figure 3).

3.3. Ecological Niche Modeling and Habitat Suitability

TSS model accuracy (Table S5) varied between 0.77 (average) for BIOCLIM (SRE), the less accurate individual model, and 0.999 for RF (random forest) for both modeled taxa. The rest of the models were all classified as excellent. The overall TSS model accuracy (average = 0.96 for both taxa) was considered excellent in predicting species distribution (TSS > 0.8). Here, seven out of eight models can be considered good to excellent, while the remaining model is classified as useful (

Table 2. Percentage of suitable and unsuitable areas in the Iberian Peninsula for the selected species and the modeled timeframes.

Species	Timeframe	SSP	Suitable (%)	Unsuitable (%)	% Reduction
U. tumidiformis	Present		7.17	92.83
	2021–2040	245	0.11	99.89	99
	2021–2040	585	0.06	99.94	99
Squalius	Present		87.07	12.93
	2021–2040	245	81.37	18.63	19
	2021–2040	585	79.40	20.60	21
	2041–2060	245	67.30	32.70	33
	2061–2080	245	58.54	41.46	42
	2081–2100	245	79.16	20.84	21

).

For U. tumidiformis, relevant climatic variables (Figure S2 are BIO1 (Annual Mean Temperature), BIO12 (Annual Precipitation), and BIO9 (Mean Temperature of Driest Quarter), associated with general climatic trends and precipitation as a proxy of water availability. BIO1 (Annual Mean Temperature), BIO5 (Max Temperature of Warmest Month), and BIO6 (Maximum Temperature of the Warmest Month and Minimum Temperature of the Coldest Month, respectively), were the most relevant variables for assessing habitat suitability for Squalius. Ensemble model performance was classified as excellent, based on the ensemble median TSS score of 0.98. The model correctly predicted 99.69% of Squalius presences (i.e., sensitivity) and 98.47% of its absences (i.e., specificity). For U. tumidiformis, the ensemble model performance was also excellent (TSS score—0.99), and 99.81% of presences and 99.21% of absences were predicted correctly.

For U. tumidiformis, in both SSPs, during 2021–2040, there is a projected reduction of 99% of suitable areas in the Iberian Peninsula, and for later time periods, there is no climatically suitable area up until the end of the century (Figure 4). For the fish hosts, there is a decrease in suitable areas throughout the century that represents a maximum decline in suitable areas of 42% (Figure S2).

4. Discussion

Our results clearly indicate that U. tumidiformis faces extinction in its known remaining range in the Iberian Peninsula. Very few areas will maintain habitat suitability in the next upcoming 20 years, and, because of their size and isolation, these are not likely to sustain viable U. tumidiformis populations. This seems to be the outcome of an already ongoing decline in the past, with a reported 82% reduction in the number of known presence sites in Portugal during the last 20 years [23].

Although we did not perform explicit demographic history modeling, our results on population structure, genetic differentiation, and genetic diversity allow us to propose some potential evolutionary history scenarios. These should be considered as a hypothesis that could be tested in the future with explicit models, provided that whole-genome data become available. Our SNP results seem to corroborate the Betic-Rif Massif origin for U. tumidiformis hypothesized using mitochondrial markers [20], specifically in the Upper Guadiana region, where populations show higher genetic diversity and lower differentiation. Furthermore, populations from all other regions presented very low genetic diversity and were genetically closer to populations from Upper Guadiana than to other geographically closer ones, suggesting founder effects and bottleneck processes. This supports the hypothesis of ancestral asynchronous migratory movements of host fish towards areas away from the Upper Guadiana basin (the proposed radiation center), and posterior isolation of the resulting populations due to paleogeomorphological rearrangements of the Iberian River network, as occurred for its fish hosts [57]. Indeed, the dispersion and differentiation of the genus Squalius in the southwestern Iberian Peninsula is well known, and we argue that it can explain the patterns observed for U. tumidiformis (Figure 5).

While our molecular results clearly support that the evolutionary history of species is the basis of its current population structure, as suggested for other organisms by [58], they do not explain why U. tumidiformis is absent from basins where it is thought to have been present in the past, or its severely fragmented distribution over a few isolated streams and rivers. This is more likely to be explained by environmental pressures and the dispersal capabilities of its host, as suggested by [26,30]). Our modeling results for the current distribution of U. tumidiformis (Figure 4) show that the binary map quite accurately depicts the present distribution of the species; however, the suitability map shows a much more extended potential area of occupation, including the Tagus and Guadalquivir basins and other areas further north. This may suggest that U. tumidiformis can only thrive in very specific optimal environmental conditions, may reflect the higher and older anthropogenic pressure in the Tagus and Guadalquivir basins when compared to the Guadiana basin [59], or point to the need for additional sampling at those basins. In fact, very recent reports have located several previously unknown populations in the Guadalquivir basin (J. Reis, unpublished data), supporting the need for additional sampling in the area. In contrast, our results show most of the Iberian Peninsula is suitable for the fish host (genus Squalius) (Figure 4). This suggests that at this scale, the fish host is not the determining factor in shaping the mussel’s distribution rather than other environmental factors. Fish hosts are widely recognized as fundamental for explaining freshwater mussels’ patterns of occurrence [32]. Nevertheless, in some instances, the distribution and decline of freshwater mussels are not strongly associated with their hosts, especially in the case of host-specific and less mobile species [26,33]. This seems to be the case in U. tumidiformis.

Modeling outputs identify almost no suitable areas for U. tumidiformis in the Iberian Peninsula as soon as in the next 20 years (Figure 4), although their hosts (fish belonging to the genus Squalius) maintain suitable areas, progressively restricted to the north of the peninsula in all our forecasts (Figure 4 and Table S6). The disappearance or northern shift has already been described for other species in the Mediterranean region e.g., [32,60,61] and stresses the increased pressure that this geographic location will suffer according to almost all climatic predictions [6]. Nonetheless, as these results were solely based on climate predictions, they should be interpreted with caution as there can be other factors at play, namely anthropogenic pressure, which might alter suitable areas. Our results for U. tumidiformis do not differ from [61], who used some catchment variables and fish hosts in addition to climatic projections as explanatory factors in SDMs, stressing the idea that climatic factors are the main constrain for the distribution of this species at the Iberian level. Other factors such as increasing pressure from invasive species, especially invasive species of predatory fish (e.g., Micropterus nigricans), crayfish (e.g., Procambarus 14usines), and bivalves (e.g., Corbicula fluminea, Dreissena polymorpha), are a growing problem in the Iberian Peninsula and affect both freshwater mussels and their native hosts [23]. In fact, it is likely that the absence of U. tumidiformis further north is due to a combination of dispersal restriction of their hosts and anthropogenic pressure such as pollution and habitat modification, rather than climatic limitations or host fish availability.

5. Implications for Conservation

Taken together, our results present a challenging scenario for the conservation of U. tumidiformis, and it is unlikely that the species will endure without active human intervention. A rescue program will be essential to prevent local extirpations due to droughts, and captive breeding will be inevitable to restore severely reduced or extinct populations (Figure 6). However, the sustainability of these actions seems unlikely even in the most favorable climate change scenarios, as streams and rivers at the current range for the species will generally not hold suitable habitat conditions. Furthermore, traditional approaches that try to preserve all the genetic lineages would imply actions on most rivers and streams at an individual level, including in captive breeding programs, which is unlikely to be feasible. Alternative approaches that defend mixing lineages from different origins for reintroduction or reinforcement programs may be needed, but no previous experiences with freshwater mussels are known. Populations from the Upper Guadiana seem to be particularly important targets for conservation actions, as they show the highest degree of diversity and, under the traditional viewpoint, would therefore be more prone to adapt to a wider array of conditions [13]. However, caution should be taken because there is evidence that highly diverse specimens may actually increase the extinction risk in small populations [14]. The evolutionary history of the species must therefore be considered in conservation actions as it influences the risk of extinction [58].

In any case, rescue and reinforcement programs must be coupled with habitat restoration aiming to reduce the impacts of climate change (Figure 6), particularly in ensuring that permanent water pools subsist even during severe droughts. In fact, the conservation and management of isolated pools in Mediterranean-type streams are key to ensuring the preservation of biodiversity associated with this habitat [62]. Climate change will challenge the effectiveness of doing so, as isolated pools tend to dry out, leading to local extirpation of species with limited dispersal capacities such as freshwater mussels. Artificial habitats that improve water persistence through drought periods have been proven useful for the conservation of aquatic organisms [63,64] and may need to be considered, but they should be addressed carefully to avoid creating ecological traps [63] (Sousa et al., 2021). Translocations and reintroductions to rivers that do not currently hold U. tumidiformis but are within its historical range and present more favorable conditions further north should be considered (Figure 6). Decisions regarding possible reinforcements and reintroductions must be taken following IUCN guidelines and after careful risk assessment [65]. Introductions outside the species’ historical range should not be considered to avoid unnecessary risks. However, favorable habitats contiguous to the known range may be considered if the absence of physical records of freshwater mussels is believed to be just a consequence of inadequate data collection, constituting a reintroduction rather than an introduction [65]. Nature-based solutions will be necessary, such as creating protected areas with strict surveillance in well-preserved stretches of each basin/sub-basin, simultaneously providing other biodiversity and human well-being benefits, and eventually coupling them with other local actions [3].

The combination of efforts to try to preserve remaining populations either by artificially increasing the resilience of their habitat to climate change or by translocating them and their offspring further north within the species historical range will imply difficult conservation decisions, as it is clear that it may not be possible to preserve all the genetic lineages of the species. Many of the possible options have little to no previous experience and their risks must be evaluated before implementation [65]. Yet, the option of not taking any action will certainly lead to extinction in the short term, effectively forcing the choice of which populations to save. For this and for other Mediterranean-type temporary stream ecosystem species, an outstanding short-term research effort seems necessary to minimize the risks of much-needed biodiversity conservation actions.