Genomic evidence for climatic adaptation in Fejervarya multistriata

Genetic diversity driven by natural selection contributes to population divergence in amphibians, thus facilitating local adaptation to climate change. Understanding the mechanisms of genetic adaptation is one of the important issues in evolutionary biology. This study set out to reveal drivers responsible for intraspecific divergence in Fejervarya multistriata and further investigate the potential involvement of selected genes in responding to climate challenges.


| INTRODUC TI ON
Climate change is directly related to the population dynamics and distribution of animals and plants (Bay et al., 2018;Mousavi-Derazmahalleh et al., 2019).Various climatic factors, including temperature and precipitation, have been related to the adaptive differentiation of populations (Guo et al., 2016;Jin et al., 2022).
Under the pressure of changing climatic conditions, flora and fauna have to evolve adaptive strategies or migrate to more suitable habitats (Chen et al., 2011;Cotto et al., 2017).Thus, genomic variation among species populations may be attributed to their improved ability to cope with rapidly changing climatic conditions.Therefore, mining climate-associated genes is vital for unravelling the adaptation mechanism and aiding efforts to accurate prediction of species responses to future climate changes (Aguirre-Liguori et al., 2021;Waldvogel et al., 2020).Amphibians, with their body temperature regulated through heat exchange with the environment, are recognised as highly susceptible to the impacts of climate change (Corn, 2005;Jiang et al., 2022;Navas et al., 2013).Consequently, considerable progress has been made in documenting the responses of amphibians to climate change.Temperature is regarded as a critical factor affecting amphibian abundance and distributions, and it influences a series of physiological processes, including metamorphosis time, growth rates and energy metabolism (Blaustein et al., 2010;Lowe et al., 2021;Mcwhinnie et al., 2021;Whittaker et al., 2007).Given the dependence of amphibians on water for their reproductive processes, lack of rainfall has been associated with reduced calling activity in anurans, failed reproduction and shrinking body size (Pincheira-Donoso et al., 2019;Semlitsch, 1987;Sheridan & Bickford, 2011;Walls et al., 2013).Overall, climate change plays a broad role in amphibian behaviour, morphology and physiology, with temperature primarily influencing phenological changes and precipitation mainly affecting population dynamics (Daszak et al., 2005;Ficetola & Maiorano, 2016).
Given that climate change has strong selective effects on climate-related traits, it is of interest to explore the genetic basis of adaptation in amphibians to better comprehend the evolutionary mechanisms responsible for climate variation (Funk et al., 2018;Sun et al., 2020).Recent advancements in DNA sequencing technology, especially the high-throughput sequencing in late 2000s, have opened up novel opportunities to investigate amphibian evolution and adaptation through genetic structure and gene function analysis (Jin et al., 2023;Khan et al., 2016;Liedtke et al., 2018;McCartney-Melstad et al., 2016;Nowoshilow et al., 2018;Shaffer et al., 2015).In a genome-scan of 10 geographically diverse Microhyla fissipes populations, researchers identified 69 annual average temperature related single-nucleotide polymorphisms (SNPs) and 248 precipitation related SNPs (Jin et al., 2022).
Another study assembled a high-quality genome of Bufo gargarizans and exhibited a significant evolutionary expansion of gene families associated with sensory perception (Lu et al., 2021).With the help of comparative genomics, the development of webbed feet in Rhacophorus kio responsible for gliding was found to be related to the Wnt signalling and vascular remodelling pathways (Wu et al., 2022).Advances in genomics have enabled understanding of the molecular mechanisms underlying environmental adaptation and organismal phenotypes.However, genomic basis of adaptation to multiple climate factors remains largely unexplored for amphibians.
The rice-paddy frog (Fejervarya multistriata) is one of the most abundant amphibian species in southern China and Southeast Asia (Huang & Tu, 2016;Lin et al., 2020).With its limited mobility and susceptibility to environmental changes, attributed to the ecological diversity of its habitat, F. multistriata is considered as an ideal model organism for studying the impacts of climate change on amphibians (Li et al., 2020;Russell et al., 2005).In this study, we conducted genome RAD-seq of F. multistriata samples collected across a biogeographical space with gradual climatic variation.We aimed to: (1) assess the population structure and genetic diversity; (2) elucidate the evolutionary history of the species; (3) identify genome-wide selective loci and candidate genes associated with climate adaptation.
The findings of our study contribute to a deeper understanding of the genetic mechanisms underlying adaptation to climate change and have potential implications for predicting the future of F. multistriata.

| Restriction-site associated DNA sequencing
Phenol/chloroform protocols were used to extract total DNA from muscle tissues (Dairawan & Shetty, 2020).To ensure DNA quality, 1% agarose gel electrophoresis was performed, and any degraded DNA samples were subjected to re-extraction and re-examination.
Subsequently, non-degraded individual samples of the same gender within each population were pooled and diluted to a concentration of 10 ng/μL.The pooled samples were then quantified and equalised to 10 ng/μL using the NanoDrop® ND-1000 spectrophotometer for the construction of RAD-seq libraries.Genomic libraries with varying insert size lengths (ranging from 300 to 500 bp) were constructed and sequenced with the Illumina HiSeq X platform, generating paired-end reads with a length of 150 base pairs.This sequencing process was performed at BGI Tech Solutions Co., Ltd (Shenzhen, China).We have submitted the raw data to the China National Centre for Bioinformation (https:// ngdc.cncb.ac.cn/ , PRJCA018392).

| Read mapping, SNP calling and filtering
To provide clean data for downstream analysis, the quality control and pre-processing of the generated raw data were performed in fastp (version 0.21) using pre-set parameters.Then, the sequence alignment between clean reads and the reference genome of F. multistriata was conducted in BWA (version 0.7.17) (Li, 2013) with pre-set parameters.
After performing a series of strict filtering steps, we utilised BCFtools (version 1.9) to produce a pileup file using the "mpileup" method, perform SNP/indel variant calling using the "call" method and merged all VCF files into a single file using the merge function.
To further improve variant quality, we kept high-quality sites from the merged VCF file that satisfied the one of the following principles: (1) individuals with a percentage of missing sites below 0.75, (2) bases with a quality score greater than 30, (3) sites with either extremely high average depth below 60 or low coverage depth greater than 4.

| Population analysis
After measuring the genomic differentiation among populations using a sliding window analysis with a 50-kb window sliding in steps of 10-kb, we calculated nucleotide diversity (π) and quantified population genetic differentiation (F ST ) among populations in PIXY (version 1.2.7) (Korunes & Samuk, 2021), and estimated the linkage disequilibrium (LD) levels for every population in PopLDdecay (version 3.4.1)(Yang et al., 2023).Prior to analyses, we generated a consensus sequence that was filtered, aligned and trimmed, and the filtering conditions were bi-allelic sites with high LD (r 2 > .1)and minor allele count (MAC < 2).Phylogenetic relationships based on SNPs were inferred via the Neighborhood-Network Algorithm (NNet) with default parameters using SPLITSTREE (version 4.16.2) (Huson, 1998;Kloepper & Huson, 2008), and the NNet network presents the results of potential evolutionary events (e.g.hybridisation, polyploidisation and recombination) in detail.
To assess the genetic structure of F. multistriata, we conducted the principal component analysis (PCA) and population structure analysis.Before analysis, we filtered the SNP data into unlinked SNP data through PLINK (version 1.90) (Lehmann et al., 2023), and set the indep-pairwise parameter to '50 10 0.1'.PCA was run with the R package LEA (version 3.2.0)(Frichot et al., 2015), and the Tracy-Widom test was conducted to define the significance of the eigenvectors.Additionally, the genomic ancestry of each frog was estimated in ADMIXTURE (version 1.3.0)(Chimusa et al., 2020) with the assumed number of ancestral clusters (K) ranging from 1 to 5.
The above generated data were visualised using self-designed R scripts.

| Demographic analysis using PSMC
The demographic histories of the 15 populations of F. multistriata were reconstructed utilising the pairwise sequentially Markovian coalescent (PSMC) model.A bootstrapping approach involving 100 replicates was employed to assess the variability among the derived Ne trajectories.We set the mutation rate (μ) at 1 × 10 −9 F I G U R E 1 Map showing locations of the 15 sites from where the Fejervarya multistriata samples were collected.The site numbers correspond to Table S1.The colours on the map indicate the two distinct groups identified through phylogenetic analyses and population structure.The red and blue dots represent the SG and HGG group, respectively.per base pair per generation and assumed a generation time (g) of 1 year (Liao et al., 2011).

| Inference of population demographic history
We inferred the population demography using Stairway plot 2. First, we obtained the folded site frequency spectrum (SFS) using the easySFS tool (Liu & Fu, 2020).Then, we used Stairway Plot 2 (version 2.1.1)to simulate the effective population size (Ne) of F. multistriata populations over time.The mutation rate and the generation time were the same as the settings in the PSMC analysis.

| Genes under selective pressure
F ST values and π ratio between the different groups were calculated using a sliding window analysis with a 50-kb window sliding in steps of 10-kb.Based on the maximum 5% of F ST values and maximum and minimum 5% of π ratio, genes under selection in both the SG group and HGG group were identified.In KOBAS (version 3.0) (Hu et al., 2023), gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used to conduct functional enrichment analysis to understand the biological processes and pathways involving the selected genes.The p-value cut-off of .05 was considered to be significantly enriched.

| Genome-climate associations
We downloaded the standard variables with 30 arc-seconds spatial resolution from http:// world clim.org/ : there are 19 bioclimatic variables in total.Through PCA and Pearson correlation analysis, climate factors can be categorised into five distinct groups: A group (BIO1, BIO5, BIO6, BIO8, BIO9, BIO10, BIO12, BIO13, BIO16), B group (BIO3), C group (BIO4, BIO7), D group (BIO2, BIO15) and E group (BIO14, BIO17, BIO18, BIO19) (Figure S1).The cumulative percentage of variance explained by principal component axis 1 and 2 was determined to be 84.99%.To further analyse these climate factors, we selected representative ones from each group, using the contribution values to principal component axes 1 and 2 as weights.Notably, BIO1 (annual mean temperature), BIO3 (isothermality), BIO7 (temperature annual range), BIO15 (precipitation seasonality) and BIO17 (precipitation of the driest quarter) exhibited the highest contribution values within their respective groups (Figure S2).Bioclimatic data corresponding to the five climate factors are added to the supplementary file (Table S9).To further explore the adaptation mechanism, we conducted linear mixed model (LMM) in GEMMA (version 0.98.5) and latent factor linear mixed model (LFMM) in the R package LEA (version 3.2.0) to find the association between the genes under selective pressure and 5 bioclimatic variables (Frichot et al., 2013;Rellstab et al., 2015;Zhou et al., 2013).

| Sampling and genome resequencing
The 300 F. multistriata frogs collected from 15 distinct locations (labelled 1-15, with 10 males and 10 females per location) across a biogeographical range in China were resequenced to understand the genetic basis of climatic adaptation in amphibians (Figure 1; Table S1).We obtained 2.75 to 11.5 Gb of clean data with a target depth of 10× for all samples.Subsequently, the clean data for the 300 individuals and two out-groups (Nanorana parkeri) were aligned to the F. multistriata reference genome and showed an average mapping rate of 98.34% (Table S2), demonstrating high uniformity between the resequencing data and the reference genome (Wang et al., 2018).Finally, a total of 21,445,171 highquality SNPs were identified from all samples and employed for further analysis (Table S3).

| Population genetic structure and demographic history
To examine the population divergence of F. multistriata spanning from Hainan province in the south to Sichuan province in the north, we constructed a phylogenetic tree for 290 samples (10 out of 20 samples from location 3 were excluded due to insufficient SNP quality) utilising genome-wide SNPs (Figure 2a).Following the phylogenetic assessment, two primary geographical groups were discerned: one consisting of 80 samples from Sichuan and Guizhou provinces (SG group, location12-15; Table S1), and another comprising 210 samples from Hainan, Guangdong and Guangxi provinces (HGG group, location1-11, Table S1).To accommodate potential admixture among individuals, we further performed population structure analysis employing ADMIXTURE, determining that K = 2 represented the optimal number with the highest clustering probability (Figure 2b; Table S4).
At K = 2, a separation between the SG and HGG groups was also evident.In agreement with these findings, the principal component analysis (PCA) demonstrated that the two groups are genetically distinct.The PCA plot illustrated that the principal components 1 (PC1) and 2 (PC2) accounted for variances of 20.5% and 6.92%, respectively, with PC1 effectively segregating the SG and HGG groups (Figure 2c).Collectively, these results substantiate the genetic divergence of F. multistriata across distinct geographic regions.
To assess the demographic history of the SG and HGG groups, we further investigated the effective population size (Ne) of the two groups over time utilising the PSMC model.The demographic history can be traced back to approximately 125 Ka, with both the SG and HGG groups experiencing a significant and continuous decline from ~125 and 106 to 73 Ka, coinciding with the last interglacial period (Figure 2d).This decline indicates a severe population bottleneck in F. multistriata populations during this period that may be caused by higher temperature.Subsequently, the effective population sizes for both groups remained relatively stable from ~73 Ka to the present, as the decline trend decelerated, coincide with the last glacial period.By utilising stairwayplot2 based on the SNP frequency spectra to infer changes in the effective population size (Ne), we confirmed an overall population decline trend (Figure S3).
During the last interglacial period, Ne declined from approximately 120 × 10 3 to 4 × 10 3 , and during the last glacial period, it further declined from 4 × 10 3 to 1.4 × 10 3 .Notably, the decline of Ne during the last interglacial period was significantly greater than that observed during the last glacial period.These findings align with the results obtained from the PSMC model.

| Signatures of selection in F. multistriata
To uncover the evolutionary drivers shaping the genomic diversity of F. multistriata, we conducted a comprehensive search of the genome for regions exhibiting selective signals between the SG and HGG groups.By employing F ST and pairwise nucleotide diversity (π) analyses, we have successfully discerned 2740 genes in F. multistriata that exhibit the maximum 5% of F ST values, as well as the maximum and minimum 5% of π ratios (Figure 3a).The HGG group, encompassing 1941 genes, exhibits an approximate 2.43 times more genes under positive selection compared to the SG group, comprising 799 genes.

| Genome-climate association analysis
The objective of this study was to explore the genetic basis underlying climate adaptation by examining the associations between selected genomic regions and various climatic factors.To achieve this, we conducted an association analysis and PCA using climate data collected from 15 distinct locations (Tables S7 and S8).We identified 5 key climatic factors, namely BIO1 (annual mean temperature), BIO3 (isothermality), BIO7 (temperature annual range), BIO15 (precipitation seasonality), and BIO17 (precipitation of the driest quarter) (Figures S1 and S2).Subsequently, we identified 193 candidate genes in selected regions associated with these 5 climatic factors based on LMM analysis.Among these genes, 77 were linked to BIO1, 50 to BIO3, 98 to BIO7, 74 to BIO15 and 79 to BIO17 (Figure 4a).
We also performed LFMM analysis to find associations between genes and environment variables.Although the number of genes detected by LFMM was comparatively lower than that detected by LMM, it is noteworthy that five genes, namely BMP2K, NRAP, ZW10, MYH1 and PLB1, overlapped with those identified by LMM (Figure S4).Among these, BMP2K exhibited associations with BIO1, BIO3 and BIO7, while NRAP showed associations with BIO1 and BIO7.Furthermore, ZW10, MYH1 and PLB1 were associated with BIO7, BIO15 and BIO17, respectively.The BMP2K gene, located on Chromosome 3, encodes bone morphology protein 2-inducible kinase (Figure 4b).This gene participates in skeletal development and growth, and its variant has close correlation with high myopia (Liu et al., 2009;Zhao et al., 2017).

| DISCUSS ION
Natural selection exerts a significant effect on the evolution and divergence of amphibian populations (Buskirk, 2009;Wei et al., 2020).
As a result, amphibians develop genetic traits to adapt to diverse environments (Richter-Boix et al., 2013;Urban et al., 2014).In this study, we performed genomic resequencing of F. multistriata from 15 geographically distinct populations (300 individuals) in China.Although the previous study conducted at the aforementioned sampling sites did not detect morphological diagnostic characters that could be employed to distinguish the populations (Zhong et al., 2008), our analysis of population structure revealed intraspecific divergence between the SG and HGG groups.Moreover, this study revealed selected genes shaping genetic diversity and identified 5 candidate genes associated with climate change adaptation in F. multistriata through a genomeclimate association analysis.Our findings provide novel perspectives on the genetic basis of climate adaptation in F. multistriata.

| Intraspecific divergence among populations
Amphibians are known to be highly susceptible to climatic and geological changes, which have profoundly affected their distribution and diversity (Araújo et al., 2008;Capinha et al., 2017;Hof et al., 2011).F. multistriata, a rice frog widely distributed in Southern China, has been the subject of taxonomic debate over whether there is differentiation between populations inhabiting distinct geographic regions (Jiang, Lv, Liu, et al., 2020;Jiang, Lv, Jia, et al., 2020;Matsui et al., 2007;Yang et al., 2022).A study utilising mtDNA D-loop sequences has revealed two distinct biogeographical groups of F. multistriata within the Chinese mainland: the Yangtze and southern lineages (Zhong et al., 2008).The Yangtze lineage was further regarded as F. kawamurai due to close genetic affinity in allozymes and supported by mtDNA divergence analysis (Djong et al., 2011;Yang et al., 2022).Consistent with these studies, we have identified significant adaptive genetic diversity within 15 populations of F. multistriata, which can be divided into two distinct groups: HGG and SG.
Notably, HGG group sampled from localities mainly distributed by the southern lineages, while SG group sampled from upper localities close to those mainly distributed by the Yangtze lineage.Drawing from these observations, we have proposed a new point that the SG populations belongs to the Yangtze lineage and therefore exhibits a different evolutionary differentiation compared to the HGG group populations, which belongs to the southern lineage.

| Demographic history is closely related to climate change
Climate change during the last glacial-interglacial period significantly impacted the demography and distribution of amphibians (Araújo et al., 2008;Crégut-Bonnoure et al., 2014;de Oliveira et al., 2021).
Last interglacial period was the warmest interglacial of the past 1.5 million years, with maximum temperatures 2 to 5°C above present in the northern hemisphere (Bakker et al., 2014;Kukla et al., 2002;Yau et al., 2016).In response to warming climate, amphibians have to accelerate the metabolic rate to meet the elevated energy demands of survival and reproduction.This adaptation led to increased risks of hunger, reproductive failure and ultimately, a higher likelihood of extinction (Carter et al., 2023;Huey & Kingsolver, 2019;Riddell et al., 2018;Trisos et al., 2020).Moreover, increased environmental temperature results in fast growth, accelerated telomere shortening and ageing rate of ectotherms, and thus leads to intergenerational effects (Burraco et al., 2020).The hierarchical mechanisms of thermal limitation (HMTL) hypothesis posited that the limited performance of amphibians in elevated temperatures could be attributed to failures at the subcellular and organ-system levels (Gangloff & Telemeco, 2018).Several studies utilising species abundance prediction model have also indicated a strong correlation between the rising extinction risk of amphibians and global warming (Pigot et al., 2023;Wagner et al., 2023).Therefore, the observed significant decline trend of F. multistriata within Southwestern China during the last interglacial period in our study may be partially attributed to increased climatic temperature.During the last glacial period, the deceleration in the decline of F. multistriata populations may be attributed to the multifaceted effects of temperature reduction on the physiology and behaviour of amphibians.First, the low temperature environment diminishes metabolic rates, resulting in reduced energy consumption and prolonged activity duration (Buckley & Huey, 2016).Furthermore, females inhabiting colder environments tend to exhibit larger body sizes and produce either more or larger eggs, or both (Morrison & Hero, 2003).Additionally, the slower growth rates of pathogens in low temperature environments contribute to a decrease in infection risk (Muths et al., 2020;Ratkowsky et al., 1982;Woodhams et al., 2008).
Conversely, the adverse effects of low temperature environments include habitat loss and range contraction, thereby exacerbating the pressure on amphibians to adapt (Luo et al., 2021;Vieira et al., 2018).
It should be noted that other factors such as precipitation, food resources and glaciation may have also influenced population dynamics during this period.However, further research is needed to provide additional evidence supporting the specific role of these factors.

| Selected pathways related to the local adaptation
Climatic and geographical changes vary extensively across our sampling locations in southern and southwestern China, where mean temperature ranged from 14.72 to 25.73°C, annual precipitation from 836 to 1896 mm, and mean altitude from 40 to 666 m (Fick & Hijmans, 2017).The adaptation of F. multistriata to local environment under selective pressure is coincide with the selection of genetic traits, consequently contributing to the divergence of populations.
Our results have identified diverse pathways related to selective pressure.The extracellular matrix organisation pathway, including LAMA2, BMP7, TLL2, LTBP1 and A2M, and so forth, primarily affects tissue structure, cell proliferation, differentiation and migration (Frantz et al., 2010;Latimer & Jessen, 2010).The PI3K-Akt signalling pathway, which includes JAK3, TGFA, PIK3R5 and NGFR, and so forth, is responsible for extra-and intracellular signal transduction that regulates cell proliferation, apoptosis, metabolism and growth (Hemmings & Restuccia, 2012;Peltier et al., 2007).It is noteworthy that climatic factors, such as temperature and ultraviolet radiation, can influence the processes mentioned above by modulating the activity of the PI3K-Akt signalling pathway (Cheng et al., 2023;Kim et al., 2006;Wu et al., 2019;Zhu et al., 2022).The metabolism of lipids pathway, including HEXB, SAMD8, ACER3 and CEPT1, and so is involved in regulation of lipid digestion, absorption, storage, cholesterol and steroid of the organism (Balk et al., 2007;Fliesler & Anderson, 1983;Large et al., 2004).Recent studies have demonstrated that elevated temperatures can induce changes in the lipid composition of fish, while cold temperatures have been shown to have a similar effect on pigs (Xu et al., 2021;Yuan et al., 2022).Another study emphasises that climate change with urbanisation increased animal lipid demand (Becker & McCluney, 2021).Unexpectedly, some genes were enriched in human papillomavirus infection pathway (e.g.PXN, ITGA1, PSMC1, VEGFA), suggesting that immunologic process is influenced by nature selection.As far as we know, however, only two papillomavirus subfamilies have been identified in amphibians, and more relevant data are still needed (Xu et al., 2021).The MAPK signalling pathway, including RASA1, MAP2K6, IKK2, SOS1 and ARAF, and so forth, plays various roles in regulation of cell cycle, transcription and translation (Chen & Thorner, 2007;Morrison, 2012;Zhang & Liu, 2002).These pathways and genes indicate that nature selection impacts a broad range of physiological processes in vivo, thereby driving population diversity in F. multistriata.

| Genes associated with climate changes
In the selected genome regions of F. multistriata, we have identified 5 genes (BMP2K, NRAP, ZW10, MYH1 and PLB1) that are potentially involved in the process of adaptation to climate change.Notably, a variant of BMP2K has been confirmed to have a strong correlation with high myopia, consequently emerging as a risk factor in the progression of myopic pathogenesis (Liu et al., 2009).For anurans, the tadpole's eye power undergoes a progressive myopic shift due to changes in the shape of the lens (Schaeffel et al., 1994;Sivak & Warburg, 1983).It was hypothesised that climate change might influence the eye development of F. multistriata tadpoles before metamorphosis through selection at BMP2K loci, thus affecting the frog's vision.Additionally, BMP2K is also involved in BMP2-induced osteoblastic differentiation (Kearns et al., 2001;Lisse et al., 2013).The protein encoded by NRAP has high expression in cardiac and skeletal muscle (Lu et al., 2003).Knockdown of NRAP can result in myofibril disassembly in developmental cardiomyocytes, and its variants has been related to cardiomyopathy and myofibrillar myopathy (Dhume et al., 2006;Koskenvuo et al., 2021;Lu et al., 2003).Similarly, MYH1 deficiency has been associated with muscle weakness accompanied by a reduction in skeletal muscle mass (Acakpo-Satchivi et al., 1997).
In mice, the loss of MYH1 and MYH4 has been found to cause skeletal muscle hypoplasia, ultimately leading to death within 30 days after birth (Hitachi et al., 2023).These findings underscore the critical role of MYH1 in skeletal muscle development and function.The ZW10 gene encodes a highly conserved protein found across eukaryotes.It plays an essential role in both mitosis and meiosis processes by participating in chromosome segregation, which is achieved through the formation of a complex with Rod and Zwilch proteins.(Chan et al., 2009;Karess, 2005).The conservation and pivotal role of ZW10 in cellular division highlight its significance in maintaining genomic stability.The protein PLB1 is a membrane-associated phospholipase exhibiting both lysophospholipase and phospholipase A2 activities (Wright et al., 2007).In C. albicans, the expression of PLB1 has been observed to be modulated by changes in temperature and other environmental conditions (Mukherjee et al., 2003).This suggests that climate change may selectively impact PLB1 expression, potentially influencing metabolic regulation in amphibians.These 5 climate-associated genes likely contribute to various physiological processes in the sensory, locomotor, cellular division and metabolism of F. multistriata, enabling its adaptability to climate change.

| CON CLUS IONS
On the basis of genome resequencing, we have revealed adaptive intraspecific divergence in F. multistriata and further identified selected regions and genes driving population differentiation.The demographic history analysis showed that the effective population sizes of F. multistriata remained relatively stable over the past 73,000 years, following a significant decline in last interglacial period.Through the examination of genome-climate associations, we have identified five candidate genes (BMP2K, NRAP, ZW10, MYH1 and PLB1) that may play a role in mitigating the impact of various climate changes.In summary, our research contributes valuable insights into understanding the genetic basis underlying climate adaptation in F. multistriata.These findings are helpful to predict the potential fate of populations in the face of climatic challenges, thus aiding in effective conservation planning.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest.

PE E R R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/ api/ gatew ay/ wos/ peer-review/ 10. 1111/ ddi.13796 .

DATA AVA I L A B I L I T Y S TAT E M E N T
DNA sequences are available from the China National Centre for Bioinformation (CNCB, PRJCA018392).

F
I G U R E 2 Population structure and demographic history analysis of Fejervarya multistriata based on genome-wide SNPs.(a) Phylogenetic tree of 15 F. multistriata populations inferred from genome-wide SNPs, with Nanorana parkeri as two outgroups.Two distinct groups, SG and HGG, were identified.SG group included samples from Sichuan and Guizhou provinces, and HGG group included samples from Hainan, Guangdong and Guangxi provinces.The sampled locations are labelled with the corresponding number.(b) ADMIXTURE plot for 290 individuals with the number of ancestral clusters (K = 2, 3, 4 and 5).Each vertical bar corresponds to an individual, and the height of each colour corresponds to the likelihood of being assigned to that cluster.The solid blue line at the bottom represents the sampling locations of the HGG group, while the solid orange line represents the sampling locations of the SG group.(c) Principal components analysis of 15 populations.The blue and orange dots represent samples from HGG and SG group, respectively.(d) The demographic history of F. multistriata HGG and SG groups depicted by the solid blue and solid orange lines, respectively [Generation time (g) = 1, mutation rate (μ) = 1 × 10 −9 ].
differentiation of BMP2K, NRAP, ZW10 and PLB1 was observed in the HGG group (F ST ≥ .441,π ratio ≥ .651),whereas MYH1 exhibited greater differentiation in the SG group (F ST ≥ .441,π ratio ≤ .393),implying their potential roles in local climate adaptation.F I G U R E 3 Genomic regions under selection in Fejervarya multistriata.(a) Distribution of F ST and π ratio values (SG group/HGG group), calculated in a 50-kb window sliding in steps of 10 kb.The regions situated to the right (π ratio ≥ .651)and left (π ratio ≤ .393)side of the grey vertical dashed lines respectively represent the maximum and minimum 5% of the π ratio distribution, and the region above the grey horizontal dashed line represents the maximum 5% of the F ST distribution (F ST ≥ .441).These two genomic regions have been identified as the selected regions for the SG group (red dots, F ST ≥ .441,π ratio ≤ .393)and the HGG group (blue dots, F ST ≥ .441,π ratio ≥ .651).(b) GO and KEGG pathway enrichment of selected genes in F. multistriata.

F
Genome-climate association analysis in Fejervarya multistriata.(a) The Venn diagram indicated selected genes associated with climatic factors (BIO1, BIO3, BIO7, BIO15 and BIO17) using the linear mixed model.(b-e) Genes BMP2K, NRAP, ZW10 and PLB1 with strong selective sweep signals in F. multistriata (HGG group).Genomic regions above the blue horizontal dashed line, representative of the top 5% of F ST values (F ST = .441),and above the red upper dashed line, representative of the top 5% of maximum π ratios (π = .651),were identified as exhibiting strong sweep signals for F. multistriata.(f) Gene MYH1 with strong selective sweep signals in F. multistriata (SG group).Genomic regions above the blue horizontal dashed line, representative of the top 5% of F ST values (F ST = .441),and below the red lower dashed line, representative of the top 5% of minimum π ratios (π = .393),were identified as exhibiting strong selective sweep signals for F. multistriata.

Financial
support was provided by the National Natural Science Foundation of China (31970393) and the Key Project of Natural Science Foundation of Sichuan Province (22NSFSC0011).FU N D I N G I N FO R M ATI O N This research was funded by the National Natural Science Foundation of China (31970393) and the Key Project of Natural Science Foundation of Sichuan Province (22NSFSC0011).