Evolutionary Clues of North-South Divergence: Genomic Adaptation of the Intertidal Mudskipper Boleophthalmus pectinirostris

Evolutionary Clues of North-South Divergence: Genomic Adaptation of the Intertidal Mudskipper Boleophthalmus pectinirostris | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Evolutionary Clues of North-South Divergence: Genomic Adaptation of the Intertidal Mudskipper Boleophthalmus pectinirostris Na Song, Ruizhi Wang, Zengman Wu, Bailin Cong, Guangxun Du, Dongdong Xiao, and 1 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-6768090/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Boleophthalmus pectinirostris is an amphibious fish widely distributed in the coastal intertidal zone, which is a good object for evolutionary biology research. With the increasing application of whole-genome resequencing technology, genetic evolution analysis and gene screening can better clarify the genetic basis and molecular mechanisms of environmental adaptation. In the present study, B. pectinirostris samples were collected across five coastal regions of China, and their population genetic viarations, demographic history alteration, and environmental adaptation mechanisms were analyzed. After alignment with the reference genome, a total of 33,016,712 high-quality SNPs were obtained for 68 specimens. The results showed that population Qingdao (QD) exhibited strong genetic heterogeneity with other four populations, and its low genetic diversity may be related to its long-term population decline. The significant changes in the effective population size of different populations over history were revealed, with notable differences between northern and southern groups. In addition, Selective sweep analysis uncovered northern-adapted genes in B. pectinirostris, linked to energy metabolism and stress response under colder environments. The results of this study contribute to a better understanding of the adaptive evolution of intertidal fishes, particularly the mechanisms associated with temperature, seasonal fluctuations and other adaptations. These findings will provide an important scientific basis for the conservation and sustainable management of this species. Biological sciences/Evolution/Population genetics/Genetic variation Biological sciences/Genetics/Genome/Genetic variation Boleophthalmus pectinirostris population genetic divergence whole-genome resequencing adaptive evolution environmental adaptation Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Introduction Boleophthalmus pectinirostris , a member of the family Gobiidae, is a distinctive amphibious fish species highly adapted to intertidal mudflat ecosystems. It is widely distributed along the western Pacific coast, spanning regions such as China, Japan, Korea, Malaysia, and Indonesia (Murdy 2017). As a warm-temperate species, B. pectinirostris thrives in coastal intertidal mudflats and demonstrates exceptional tolerance to a wide range of salinities (Takegaki 2008). B. pectinirostris has developed a unique feeding strategy, primarily consuming benthic diatoms. This dietary specialization further influences its distribution in aquatic environments (Yang et al. 2003; Polgar et al. 2009). Due to its relatively low trophic position and short food chain, B. pectinirostris can rapidly respond to environmental changes, making it a valuable bioindicator for coastal ecosystem monitoring (Pan et al. 2022). Moreover, it has been recognized as a key species for assessing the ecological health of mudflat and mangrove ecosystems (Polgar 2009). Despite growing interest in the genetic characteristics of B. pectinirostris , previous studies have predominantly relied on traditional molecular markers, such as mitochondrial DNA sequences and microsatellites, to explore its genetic diversity, phylogenetic relationships, and population differentiation (Kanemori et al. 2006; Gong et al. 2017). However, these methods often suffer from limited resolution and insufficient sample coverage, hindering a comprehensive understanding of population genetic structure and adaptive evolution. In contrast, whole-genome resequencing (WGS) offers a powerful approach for analyzing genomic variations, including single nucleotide polymorphisms (SNPs) and structural variants, by aligning sequencing reads to a reference genome (Wang et al. 2024). With the rapid advancement of sequencing technologies, an increasing number of marine fish species have been fully sequenced, and WGS has emerged as a critical tool for elucidating population genetic diversity and local adaptation (Xiong et al. 2023). Furthermore, downstream analyses enable the identification of candidate genes and pathways associated with environmental adaptation, providing insights into the molecular mechanisms underlying ecological resilience in fish species (Li et al. 2023; Sun et al. 2023; Huang et al. 2024). Intertidal organisms such as B. pectinirostris respond to extreme environmental challenges, including temperature changes caused by Marine heat waves (Wang et al. 2020; Stillman et al. 2025). These selective pressures are likely to drive the accumulation of adaptive genetic variations through polygenic selection (You et al. 2014). Previous whole-genome sequencing studies on the intertidal snail Littorina brevicula have revealed two distinct genetic lineages along a latitudinal thermal gradient, with north-south divergence driven by thermal environmental variation (Lin et al., 2024). Similarly, mudskippers, as emblematic intertidal organisms exhibiting amphibious adaptations, demonstrate heightened sensitivity to environmental fluctuations, making them an ideal model for studying adaptive evolution in dynamic habitats. Previous studies based on traditional genetic markers (e.g., mitochondrial DNA and microsatellites) have shown no significant population differentiation within the southern Chinese populations, but revealed two distinct lineages between southern Chinese populations and those from Japan and South Korea (Wang et al. 2012; Chen et al. 2015). However, the limited resolution of these markers, combined with the lack of research on northern Chinese populations, has hindered a comprehensive understanding of population divergence and the genetic basis of their adaptive potential. Given the geographical proximity of Japanese and South Korean populations to northern China, we hypothesized that significant genetic structure may exist between northern and southern Chinese populations of B. pectinirostris , potentially reflecting distinct adaptive evolutionary responses to local environmental conditions. To address these gaps, in this study, we employed whole-genome resequencing to analyze populations of B. pectinirostris from five coastal locations, including a northern population from Qingdao, China. Using population genomic data, this study aims to unravel the demographic history, genetic diversity patterns, and adaptive genetic mechanisms of B. pectinirostris . By analyzing genome-wide variations, we can not only more accurately characterize population structure and divergence history but also identify key genomic regions associated with local adaptation. These findings will provide new insights into the adaptive evolution of intertidal species and offer a scientific basis for the conservation and management of this ecologically significant species. Materials And Methods Sample collection, DNA extraction, and genome resequencing A total of 68 B. pectinirostris specimens were collected from five coastal regions in China: Qingdao (QD) in Shandong Province, Ganyu (GY) in Jiangsu Province, Wenzhou (WZ) in Zhejiang Province, Xiamen (XM) in Fujian Province, and Zhanjiang (ZJ) in Guangdong Province (Fig. 1). For each region, 10–15 specimens were collected and identified through morphological examination and COⅠ gene analysis before being used for subsequent whole-genome studies. Dorsal muscle tissue was dissected from each specimen, preserved in absolute ethanol, and stored at -80°C in an ultra-low-temperature freezer until genomic DNA extraction. Genomic DNA was extracted using the standard phenol-chloroform method. The extracted DNA was randomly fragmented using an ultrasonic disruptor to construct sequencing libraries. The library preparation process included end repair, poly(A) tailing, adaptor ligation, and PCR amplification. The quality and concentration of the genomic libraries were initially assessed using a Qubit 2.0 fluorometer, followed by verification of insert fragment sizes using an Agilent 2100 Bioanalyzer. Finally, paired-end sequencing was performed on the Illumina HiSeq PE150 platform, generating reads with an average length of 150 bp. All experimental procedures were conducted at Novogene Bioinformatics Technology Co., Ltd. Sequencing data quality control and sequence alignment The raw sequencing data contained adaptor sequences and low-quality bases, which could compromise subsequent analyses. To ensure data quality, raw reads were filtered using Fastp software (Chen et al. 2018) to remove adaptors and filter out low-quality sequences, resulting in high-quality clean reads. For this study, the first chromosome-level genome of B. pectinirostris (GenBank accession number: GCA_026225935.1) was used as the reference genome (Bian et al. 2023). Clean reads from each individual were aligned to the reference genome using BWA (Li and Durbin 2009) with the parameters “mem -t 4 -k 32 -M.” The alignment results were saved in BAM format. Subsequently, SAMtools (Li et al. 2009) was employed to sort the BAM files, and duplicate reads were removed using Picard v2.18.17. SNP detection and annotation Based on the alignment results, variant calling was performed using GATK (McKenna et al. 2010) to generate a VCF file containing raw SNPs. The raw SNPs were filtered using stringent quality control criteria (parameters: QD < 2.0 || MQ 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0). High-quality SNPs retained after filtering were annotated using ANNOVAR software (Wang et al. 2010). The SNPs were categorized based on their genomic locations (e.g., upstream, downstream, intronic, intergenic, and splicing regions) and functional effects (e.g., synonymous and non-synonymous variants) (Jones et al. 2014). Population genetic structure and genetic diversity analysis After filtering, a set of high-quality SNPs was obtained for subsequent population genetic analyses. Three methods were employed to investigate population genomic relationships: Neighbor-Joining (NJ) tree construction, Principal Component Analysis (PCA) and ADMIXTURE analysis. First, a genetic distance matrix was calculated using Phylip software, and a phylogenetic tree was constructed using the Neighbor-Joining (NJ) method (Felsenstein 1993). The resulting tree was visualized using iTOL (https://itol.embl.de/). PCA was conducted on whole-genome SNPs from all 68 individuals using PLINK v1.9 (Purcell et al. 2007). Additionally, population genetic structure and lineage information were inferred using ADMIXTURE v1.3 software, with K values ranging from 2 to 6 (Alexander et al. 2009). To evaluate patterns of linkage disequilibrium (LD) decay across the genomes of the five populations, allele frequency correlation (r²) was calculated using PopLDdecay v3.42 (Zhang et al. 2019). To further explore genetic variation between northern and southern populations of B. pectinirostris , we assessed the genome-wide distribution of nucleotide diversity ( π ) and the genetic differentiation index ( F st ) using VCFtools v0.1.16. These analyses were performed with a sliding window approach, using a window size of 100 kb and a step size of 10 kb (Danecek et al., 2011). Effective population size analysis of evolutionary history To reconstruct the demographic history of different populations, Paired Sequentially Markovian Coalescent (PSMC) analysis was employed to infer changes in the effective population size of B. pectinirostris over its evolutionary history (Li et al. 2011). Input files for PSMC modeling were generated using the “fq2psmcfa” and “splitfa” tools within the PSMC software. The PSMC analysis was conducted with the following parameters: “-N25” for the number of algorithm iterations, “-t15” as the upper limit for the time to the most recent common ancestor (TMRCA), “-r5” for the initial θ/ρ ratio, and “-p 4 + 25 * 2 + 4 + 6” for atomic time intervals. The reconstructed population history was visualized using the “psmc_plot.pl” script. A mutation rate of 3.51×10⁻⁹ per generation and a generation time of 1 year were assumed for the analysis (You et al. 2014). Additionally, TreeMix v1.13 was used to infer the maximum-likelihood (ML) tree for the five populations. A window size of 1000 was applied to account for linkage disequilibrium (–k), and the “–global” option was used to generate the ML tree. Migration events were sequentially added to the tree by testing models with 0 to 4 migration events (–m 0 1 2 3 4) (Pickrell et al. 2012). Selection signature detection and functional enrichment analysis The π ratio between the two populations was calculated based on the previously estimated nucleotide diversity ( π ). The results of the population differentiation coefficient ( F st ) and π ratio were then integrated and summarized (Danecek et al. 2011). The northern population (QD) was compared collectively against the four southern populations (GY, WZ, XM, and ZJ) in the analysis. Genomic windows with values in the top 5% for both F st and π ratio were identified as candidate regions under strong selective pressure. Based on these candidate regions, corresponding SNPs and genes were annotated. Additionally, we performed an XP-EHH (cross population extended haplotype homozogysity) analysis using SELSCAN v1.3.0 with default settings to detect ongoing or nearly fixed selective sweeps by comparing haplotype patterns across the five populations (Szpiech et al. 2014). For the XP-EHH selection scan, genes located within 10 kb upstream and downstream of the identified loci were extracted for further analysis. For genes exhibiting selection signals in each population, functional annotation was performed using eggNOG-mapper software (Cantalapiedra et al. 2021) to assign Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. GO term classification and KEGG pathway enrichment analysis were further conducted using OmicShare tools (http://www.omicshare.com/tools) (Ashburner et al. 2000; Jones et al. 2014; Kanehisa et al. 2023). Results Sequencing and SNP calling The study conducted library construction and sequencing on 68 samples from five geographic populations of B. pectinirostris . A total of approximately 821.72 Gb of raw sequencing data was generated, with an average of 12.08 Gb per sample. The average sequencing depth was 12.54×, with an average GC content of 40.50%. The average Q20 and Q30 scores were 96.49% and 91.65%, respectively (Fig. S1). After quality filtering, approximately 805.11 Gb of high-quality clean data was retained for downstream analyses. By aligning the sequencing data with the reference genome, a total of 33,016,712 SNP variants were identified. The distribution of these SNPs provided valuable insights into the genomic variation across the population. The majority of the variants (63.6%) were located in intergenic regions, reflecting the extensive noncoding regions of the genome. Intronic variants accounted for 18.7%, while coding variants, including 200,311 synonymous and 490,528 non-synonymous mutations, were less frequent. Regulatory regions, encompassing both upstream and downstream areas, represented 15.5% of the SNPs, suggesting their potential role in gene expression regulation. Variants in splicing-related regions and untranslated regions (UTRs) were rare. To further validate the reliability of the data, transition/transversion (Ts/Tv) analysis will be conducted (Table 1). These findings provide a comprehensive overview of the SNP landscape within this population. Additionally, we analyzed the distribution of core SNP markers across the chromosomes. The SNP density distribution map revealed that SNPs are evenly distributed across all chromosomes (Fig. 2). Table 1 The distribution of SNP variants in the genome region. Genetic element SNP numbers Upstream 3,068,552 Downstream 3,061,303 CDS 753,445 Intron 7,335,438 Intergenic 25,129,719 Splicing 53,847 5' UTR 233 Synonymous 200,311 Non-synonymous 490,528 Transitions (Ts) 17,040,364 Transversions (Tv) 11,466,334 Ts/Tv Ratio 1.49 Population genetic structure analysis and gene flow To investigate the genetic structure of B. pectinirostris across different populations, phylogenetic analysis was conducted using high-quality SNP loci filtered from the sequencing data of samples collected from five geographic populations. The resulting phylogenetic tree revealed that populations generally clustered according to their geographic origins, with the exception of two individuals from the ZJ population, which grouped with the WZ population (Fig. 3A). A three-dimensional PCA plot was constructed based on the first (PCA1), second (PCA2), and third (PCA3) principal components. The PCA results indicated that the XM and ZJ populations clustered together, suggesting a closer genetic relationship and higher genetic similarity between these two populations. In contrast, the QD and GY populations formed distinct clusters, reflecting greater genetic divergence from the other populations. The WZ population occupied an intermediate position between the other groups, consistent with its geographic location (Fig. 3B). The Admixture analysis further supported the results of the phylogenetic tree and PCA. When K=2, the cross-validation (CV) error was minimized, and the populations could be broadly divided into two major clusters: Qingdao formed one cluster, while XM and ZJ clustered together. The GY and WZ populations were positioned between these two groups. As K increased, finer genetic distinctions between the clusters were observed. For K values between 3 and 6, QD remained a distinct cluster, while the WZ population showed a genetic composition more similar to the XM and ZJ populations, indicating closer genetic relatedness between them (Fig. 3C). We also calculated the F st between populations based on all SNPs and then averaged them. Both the other four populations showed significant genetic divergence from the QD population (Fig. S2). This comprehensive genetic analysis provides valuable insights into the population structure and evolutionary relationships of B. pectinirostris . Based on the results of genetic differentiation and cluster analysis, the gene flow was performed. Gene flow analysis showed that when three migration events occurred, the QD population gene flowed to GY, and GY population gene flowed to WZ, which was similar to the results of cluster analysis. This implies that the QD and GY populations perhaps originated from a single ancestor. It is worth noting that gene flow from GY to ZJ was detected despite their geographic separation, possibly due to anthropogenic influences (Fig. 4). Linkage disequilibrium analysis and history dynamics Based on the results of genetic differentiation and cluster analysis, the linkage disequilibrium analysis was performed. LD analysis was conducted for five populations. The attenuation distance of the four southern populations was smaller than that of the Qingdao population, which may be due to the high genetic diversity of these populations. This was consistent with the results of genetic diversity analysis (Fig. 5). Based on the PSMC model, we reconstructed the changes of the historical effective population size ( N e ) of four different geographic groups of B. pectinirostris to infer the population evolution history. The results showed that the N e changes of different populations were basically consistent about 150,000 years ago, and they all experienced a population contraction about 700,000 to 800,000 years ago. After that, the N e of QD and GY populations in the north was slightly larger than that of other populations in the south. In the period after 150,000 years ago, the trajectory of N e change of QD populations and GY populations deviated from that of the other three populations, and the historical effective population size of the other three groups showed a downward trend, while that of the other three populations showed an upward trend. In general, a clear latitudinal gradient in effective population size ( N e ) was evident by 10,000 years ago, with southern populations maintaining higher values than their northern counterparts, and the QD population exhibiting the lowest N e . (Fig. 6). Genome-wide selective signals and candidate gene analysis In this study, we identified genomic regions under selection using F st & π to uncover genes associated with environmental adaptation. Considering that there were no significant genetic differences between the southern populations, we selected the southern populations as control groups to reveal environmental selection genes in the northern population. With four southern populations as the control groups, we screened 1685 (GY, z ( F st ) ≥ 1.84 and log 2 ( π QD/GY) ≥ 1.36), 1964 (WZ, z ( F st ) ≥ 1.86 and log 2 ( π QD/GY) ≥ 1.50), 1859 (z ( F st ) ≥ 1.81 and log 2 ( π QD/GY) ≥ 1.47) and 1892 (ZJ, There were z ( F st ) ≥ 1.84 and log 2 ( π QD/GY) ≥ 1.54) candidate genes (Fig. 7). Therefore, we identified 977 overlapping genes from the QD population screened with the southern population as the control group as potential genes related to adaptation to the northern environment (Fig. S3A). Similarly, with QD population as the control group, we obtained 381, 98, 273 and 109 candidate genes in GY, WZ, XM and ZJ populations (Fig. 7). A total of 14 overlapping genes were identified as potential genes associated with southern environmental adaptation (Fig. S3B). In order to fully understand the molecular functions of these candidate genes related to environmental adaptation, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for these candidate genes. Since the southern population exhibited fewer overlapping genes, all subsequent analyses primarily focused on the candidate genes from the northern population. The GO analysis results indicated that the major GO terms associated with adaptation to the northern environments included Carboxylic acid biosynthetic process (GO:0046394) and Organic acid biosynthetic process (GO:0016053) in Biological Processes, as well as ATP binding (GO:0005524), Adenyl nucleotide binding (GO:0030554), and Adenyl ribonucleotide binding (GO:0032559) in Molecular Function (Fig. 8A). The KEGG enrichment results revealed that the candidate genes were significantly enriched in several metabolic pathways, including the PI3K-Akt signaling pathway (ko04151), Tight junction (ko04530), Rap1 signaling pathway (ko04015), Platelet activation (ko04611), Regulation of actin cytoskeleton (ko04810), and Focal adhesion (ko04510) (Fig. 8B). Genome-wide selection signal and candidate gene analysis based on XP-EHH Based on the XP-EHH values at each locus, we selected the top 1% regions as candidate regions. After filtering, we detected strong selection signals between the QD population and the other four populations. A total of 105 genes were identified between the QD and GY populations (Fig. S4A). Meanwhile, 314, 312, and 321 genes were identified between the QD and WZ, XM, and ZJ populations, respectively (Fig. S4B; Fig. S4C; Fig. S4D). KEGG enrichment analysis was performed on the genes screened from these four comparisons. The results showed that the genes between the QD and WZ populations, as well as the QD and ZJ populations, were significantly enriched in pathways such as Hepatitis C (ko05160), Two-component system (ko02020), Phototransduction-fly (ko04745), Huntington disease (ko05016), and Tight junction (ko04530) (Fig. 9). Among these, the Hepatitis C and Tight junction pathways were consistent with the F st & π analysis. Additionally, the genes between the QD and GY populations, as well as the QD and XM populations, were significantly enriched in pathways such as Two-component system (ko02020), Phototransduction-fly (ko04745), Amoebiasis (k0), and Nitrogen metabolism (ko00910) (Fig. S5). Discussion Genome-wide SNPs decipher phylogeographic structure Single nucleotide polymorphisms (SNPs) have been well established as a major source of genetic variation in model organisms (Langley et al. 2012). In the present study, we identified 33,016,712 high-quality SNPs, confirming their prevalence as important molecular markers. These genomic resources enable comprehensive investigations of genetic diversity, population structure, and environmental adaptation in natural populations. Such analyses provide critical insights for multiple applications, including germplasm resource conservation, phylogenetic reconstruction, and evolutionary studies (Zhou et al. 2023). Using three complementary methods—PCA, ADMIXTURE, and NJ tree analysis—we identified clear geographic clustering patterns among the five populations. Specifically, the QD population formed a distinct cluster, showing limited gene flow with the geographically proximate GY population. In contrast, the ZJ and XM populations exhibited closer genetic relationships, while the WZ population occupied an intermediate genetic position between the northern and southern groups. Whole-genome SNP analysis revealed significant population genetic structure and differentiation between the QD population in Shandong Province and the other four populations. Additionally, higher F st values were observed within the southern group, consistent with previous findings based on cytochrome b analyses (Chen et al. 2015). According to principles of genetics and evolutionary biology, genetic diversity is the foundation for species to adapt to their environments over time and to evolve in response to environmental changes. It is also a fundamental determinant of a species' evolutionary potential (Ma et al. 2022). Our results showed that the ZJ population exhibited the highest nucleotide diversity, while the QD population had the lowest. Previous studies have also indicated that B. pectinirostris generally exhibits low genetic diversity, with significant variation among populations and a gradual decline in diversity from south to north, which aligns with the findings of this study (Chen et al. 2014). In this study, the NJ tree revealed that all individuals clustered into three major branches: the southernmost WZ, XM, and ZJ populations formed one branch, the QD population formed another, and the GY population occupied an intermediate branch between the two. This clustering pattern aligns with the geographical distribution of the populations. Previous analyses based on the COI gene have also indicated frequent gene flow among B. pectinirostris populations across several southern coastal provinces, including Guangdong, Fujian, and Jiangsu, which is consistent with the findings of this study (Yang et al. 2012). Population history and evolutionary dynamics The population genetic structure and genetic diversity of species are closely tied to their historical demographic dynamics (McKeown et al. 2019; Wang et al. 2020). In this study, we reconstructed the demographic history of different geographic populations of B. pectinirostris using the PSMC model. The results revealed that, prior to the early Pleistocene of the Quaternary, the effective population sizes of the various B. pectinirostris populations were relatively stable. However, during the Glacial-Interglacial Cycle, B. pectinirostris experienced both population bottlenecks and expansion events. Notably, between 120,000 and 150,000 years ago, the demographic trajectories of the populations began to diverge. During this period, the QD and GY populations started to separate from the other three populations, with the QD population undergoing a more rapid decline in effective population size compared to the GY population. This suggests that genetic differentiation between the QD population and the other four populations likely began around this time. Additionally, we observed that during the Last Glacial Period, extreme environmental conditions led to a significant reduction in the effective population sizes of all populations, a pattern consistent with findings in many other marine species (Wu et al. 2022). Furthermore, over the past 10,000 years, the QD population exhibited a notably smaller effective population size compared to the other populations, which is partially supported by the results of genetic diversity analysis. This may be attributed to the QD population being in a state of continuous decline without experiencing population expansion, resulting in lower genetic diversity (Bernos et al. 2022; Gao et al. 2024). Strong selective signals based on F st & π underlie north-south adaptation divergence Recent studies have demonstrated that whole-genome resequencing technology provides a powerful tool for investigating the environmental adaptation of fish populations across different geographic regions (Kon et al. 2021; Li et al. 2023; Huang et al. 2024). Previous studies based on mitochondrial gene fragments suggested a lack of distinct genetic structure among B. pectinirostris populations along the Chinese coast (Chen et al. 2015). In contrast, this study revealed significant genetic differentiation among B. pectinirostris populations in coastal China using whole-genome SNP markers. We hypothesize that this differentiation may result from the species' wide distribution range and limited dispersal ability, making it highly susceptible to local environmental heterogeneity (Zhang et al. 2012). Our genetic structure analyses, F st -based population differentiation metrics, and historical demographic inference indicated that the QD population exhibited the most pronounced genetic divergence from the other four populations. Therefore, we selected the QD population as the reference group and compared it with the four southern populations (GY, WZ, XM, and ZJ) to identify genomic regions associated with environmental adaptation. We detected 997 candidate genes in the northern population and 14 in the southern population. The higher number of selected genes in the northern population suggests stronger natural selection pressures, potentially driven by colder temperatures or greater seasonal variability. Populations under intense natural selection often exhibit lower genetic diversity, which aligns with our finding that the QD population had the lowest genetic diversity (Goh et al. 2023). Compared to the warm, humid climate and mangrove-rich habitats in the south, the northern environment is harsher, requiring B. pectinirostris to employ more diverse strategies to cope with environmental challenges (Zhang et al. 2014). Given the weaker selection signals in the southern populations, we focused on analyzing the selection signals in the northern population. Functional annotation and clustering of the selected genes revealed their involvement in energy metabolism regulation, DNA damage repair, and signal transduction. Notably, the PI3K-Akt signaling pathway, known for its role in skin wound healing, has not been extensively studied in marine fish (Yang et al. 2024). This pathway is also critical for regulating cellular processes such as mitosis and melanin production (Martini et al. 2014). Activation of the PI3K/Akt signaling pathway has been shown to reduce melanin production by downregulating genes such as TYR and TYRP1 (Shin et al. 2014; Tang et al. 2023). The selection of genes related to skin pigmentation and DNA damage repair may reflect differences in photoperiod and light intensity between northern and southern environments. Additionally, the mangrove ecosystems in the south may provide protection against ultraviolet radiation, a factor absent in the north (Coelho et al. 2010; Qian et al. 2013). The selection of energy metabolism-related genes may be linked to the QD population's adaptation to higher-latitude environments, where they undergo a 3–4 month hibernation period in caves during winter. This requires precise regulation of energy metabolism to survive low temperatures and other harsh conditions (Nie et al. 2022). The enrichment of GO terms related to ATP binding further supported this hypothesis. In conclusion, our adaptive evolution analysis identified candidate genes associated with the adaptation of B. pectinirostris to northern and southern environments. These findings provided valuable insights for further research into the evolutionary history and environmental adaptation mechanisms of this species. Population selection signals based on XP-EHH uncovers vision and immunity-related adaptations In addition to selective sweep analysis, this study employed the XP-EHH method to detect strong selection signals between the northern and southern populations. A total of 105 to 321 candidate genes were identified between the QD population and the other four populations. These genes were primarily associated with immune response, signal transduction, and DNA damage repair. KEGG enrichment analysis revealed that genes identified in comparisons between the QD and WZ populations, as well as the QD and ZJ populations, were enriched in the Tight Junction pathway. Additionally, genes from all four comparisons were enriched in the Phototransduction-fly pathway, which has been implicated in visual adaptation in fish and may reflect the visual adaptation of Boleophthalmus pectinirostris during terrestrial activities (Zhou et al. 2025). Notably, the Guanine Nucleotide-Binding Protein (GNAQ) gene and the G protein-coupled receptor 18 (GPR18) gene were identified among the candidate genes. GNAQ, a regulator or transducer, plays a role in various transmembrane signaling systems (Haase et al. 2021). In fish, GNAQ has also been linked to olfactory development and immune function (Wang et al. 2022). Similarly, GPR18 has been shown to be closely associated with immune responses in fish, providing protection against certain pathogenic infections (Pridgeon et al. 2013). Furthermore, LPSenhR-1, a member of the same GPR family, has been demonstrated to play a role in immunity in Oncorhynchus mykiss (Kiryu et al. 2003). These findings suggest that these genes may contribute to the adaptation of northern populations to their environment. Conclusions This study employed whole-genome sequencing to elucidate the genetic divergence patterns among five geographic populations of B. pectinirostris along China's coastline. The results indicated that there is a pronounced genetic isolation and reduced diversity in the northern Qingdao (QD) population, contrasting with higher gene flow in southern populations. Selective sweep analysis revealed strong signatures of selection in the northern population, particularly in genes associated with energy metabolism and immune response, reflecting adaptive evolution to colder climates and intense UV radiation. In contrast, southern populations exhibited weaker selection signals, likely due to the buffering effect of mangrove habitats. These findings provide compelling evidence that both geographic isolation and environmental gradients drive population differentiation in this ecologically important species. From a conservation perspective, our results highlight the urgent need for prioritized protection of the genetically distinct QD population to prevent further loss of genetic diversity, while maintaining gene flow among southern populations to preserve their adaptive potential. This study offers novel genomic insights into the mechanisms of adaptive evolution in intertidal organisms and establishes a foundation for evidence-based conservation strategies. Declarations Acknowledgements This study was funded by the National Key Research and Development Program of China (2023YFD2401103). Thank all the members who participated in this research for their contributions to this paper Research Ethics Statement All experimental procedures involving Boleophthalmus pectinirostris were conducted in accordance with relevant guidelines and regulations. The animal study was reviewed and approved by the Ocean University of China Academic Committee. Data archiving The genome resequencing data for Boleophthalmus pectinirostris have been deposited at NCBI under the BioProject accession number PRJNA1258070. Conflict of interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References Murdy EO, Jaafar Z (2017) Taxonomy and systematics review. In: Fishes Out of Water, 1st edn. CRC Press, pp 36 Takegaki T (2008) Threatened fishes of the world: Boleophthalmus pectinirostris (Linnaeus 1758) (Gobiidae). Environ Biol Fish 81:373–374 Yang KY, Lee SY, Williams GA (2003) Selective feeding by the mudskipper ( Boleophthalmus pectinirostris ) on the microalgal assemblage of a tropical mudflat. Mar Biol 143:245–25 Polgar G, Crosa G (2009) Multivariate characterisation of the habitats of seven species of Malayan mudskippers (Gobiidae: Oxudercinae). Mar Biol 156:1475–1486 Pan Y, Tian L, Zhao Q, Tao Z, Yang J, Zhou Y et al. (2022) Evaluation of the acute toxic effects of crude oil on intertidal mudskipper ( Boleophthalmus pectinirostris ) based on antioxidant enzyme activity and the integrated biomarker response. Environ Pollut 292:118341 Polgar G (2009) Species-area relationship and potential role as a biomonitor of mangrove communities of Malayan mudskippers. Wetl Ecol Manag 17:157–164 Kanemori Y, Takegaki T, Natsukari Y (2006) Genetic population structure of the mudskipper Boleophthalmus pectinirostris inferred from mitochondrial DNA sequences. Jpn J Ichthyol 53:133–141 Gong W, Huang Y, Xie J (2017) Genome-wide identification of novel microRNAs from genome sequences using computational approach in the mudskipper ( Boleophthalmus pectinirostris ). Russ J Bioorg Chem 43:397–408 Wang QL, Zhang JY, Liu ZS, Duan YY, Li CY (2024) Integrative approaches based on genomic techniques in the functional studies on enhancers. Brief Bioinform 25:bbad442 Xiong X, Liu JK, Rao YL (2023) Whole genome resequencing helps study important traits in chickens. Genes 14:1198 Li L, He S, Lin MH, Zhang YP, Kuhl H, Liang XF (2023) Whole-genome resequencing and bisulfite sequencing provide new insights into the feeding habit domestication in mandarin fish ( Siniperca chuatsi ). Front Genet 13:1088081 Sun CF, Zhang XH, Dong JJ, You XX, Tian YY, Gao FY et al. (2023) Whole-genome resequencing reveals recent signatures of selection in five populations of largemouth bass ( Micropterus salmoides ). Zool Res 44:78–89 Huang WW, Lai J, Liang WQ, Ye SZ, Li JW, Zhou JW et al. (2024) Identification of sex-specific DNA markers in the army fish ( Spinibarbus hollandi ) by whole genome re-sequencing method. Aquaculture 583:740605 Wang J, Ma LX, Dong YW (2023) Coping with harsh heat environments: molecular adaptation of metabolic depression in the intertidal snail Echinolittorina radiata . Cell Stress Chaperon 28:477–491 Stillman JH, Amri AB, Holdreith JM, Hooper A, Leon RV, Pruett LR et al. (2025) Ecophysiological responses to heat waves in the marine intertidal zone. J Exp Biol 228:JEB246503 You X, Bian C, Zan Q (2014) Mudskipper genomes provide insights into the terrestrial adaptation of amphibious fishes. Nat Commun 5:5594 Lin XN, Ma CY, Hu LS, Liao ML, Ma LX, Teske PR et al. (2024) Genomics-informed range predictions under global warming reveal reduced adaptive diversity whilst buffering range shifts for a marine snail. Glob Change Biol 30:e17571 Wang CT, Liu ZZ, Ma LB (2012) Isolation and characterization of 16 polymorphic microsatellite markers for the mudskipper ( Boleophthalmus pectinirostris ). Conserv Genet Resour 4:1059–1061 Chen W, Hong W, Chen S (2015) Population genetic structure and demographic history of the mudskipper Boleophthalmus pectinirostris on the northwestern Pacific coast. Environ Biol Fish 98:845–856 Chen S, Zhou Y, Chen Y, Gu J (2018) fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34:i884–i890 Bian C, Huang Y, Li RH, Xu PW, You XX, Lv YY et al. (2024) Genomics comparisons of three chromosome-level mudskipper genome assemblies reveal molecular clues for water-to-land evolution and adaptation. J Adv Res 58:93–104 Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25:1754–1760 Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N et al. (2009) The Sequence Alignment/Map format and SAMtools. Bioinformatics 25:2078–2079 McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A et al. (2010) The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res 20:1297–1303 Wang K, Li M, Hakonarson H (2010) ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res 38:e164 Jones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C et al. (2014) InterProScan 5: genome-scale protein function classification. Bioinformatics 30:1236–1240 Felsenstein J (1993) PHYLIP (phylogeny inference package), version 3.5 c. University of Washington, Seattle. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D et al. (2007) PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet 81:559–575 Alexander DH, Novembre J, Lange K (2009) Fast model-based estimation of ancestry in unrelated individuals. Genome Res 19:1655–1664 Zhang C, Dong SS, Xu JY, He WM, Yang TL (2019) PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics 35:1786–1788 Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA et al. (2011) The variant call format and VCFtools. Bioinformatics 27:2156–2158 Li H, Durbin R (2011) Inference of human population history from individual whole-genome sequences. Nature 475:493–496 Pickrell JK, Pritchard JK (2012) Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet 8:e1002967 Szpiech ZA, Hernandez RD (2014) selscan: An Efficient Multithreaded Program to Perform EHH-Based Scans for Positive Selection. Mol Biol Evol 31:2824–2827 Cantalapiedra CP, Hernández-Plaza A, Letunic I, Bork P, Huerta-Cepas J (2021) eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale. Mol Biol Evol 38:5825–5829 Ashburner M, Ball C, Blake J (2000) Gene Ontology: tool for the unification of biology. Nat Genet 25:25–29 Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M (2023) KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res 51:D587–D592 Langley CH, Stevens K, Cardeno C, Lee YCG, Schrider DR, Pool JE et al. (2012) Genomic variation in natural populations of Drosophila melanogaster. Genetics 192:533–598 Zhou J, He W, Wang J (2023) The pan-plastome of tartary buckwheat ( Fagopyrum tataricum ): key insights into genetic diversity and the history of lineage divergence. BMC Plant Biol 23:212 Ma DN, Lai ZF, Ding QS, Zhang K, Chang KZ, Li SH et al. (2022) Identification, characterization and function of orphan genes among the current Cucurbitaceae genomes. Front Plant Sci 13:872137 Chen H, Polgar G, Yin W, Fu CZ (2014) Cryptic species and evolutionary history of Boleophthalmus pectinirostris complex along the northwestern Pacific coast. Acta Hydrobiol Sin 38:75–86 Yang F, He LJ, Lei GC, Zhang AB (2012) Genetic diversity and DNA barcoding of mudskipper common species along southeast coasts of China. Chin J Ecol 31:676–683 McKeown NJ, Arkhipkin AI, Shaw PW (2019) Genetic analysis reveals historical and contemporary population dynamics in the longfin squid Doryteuthis gahi: implications for cephalopod management and conservation. ICES J Mar Sci 76:1019-1027 Wang D, Gao L, Tian H, et al. (2020) Population genetics and sympatric divergence of the freshwater gudgeon, Gobiobotia filifer, in the Yangtze River inferred from mitochondrial DNA. Ecol Evol 10:50–58 Wu B, Chen X, Yu M, Ren J, Hu J, Shao C et al. (2022) Chromosome-level genome and population genomic analysis provide insights into the evolution and environmental adaptation of Jinjiang oyster Crassostrea ariakensis . Mol Ecol Resour 22:1529–1544 Bernos TA, Avlijaš S, Hill J, Morissette O, Ricciardi A, Mandrak NE et al. (2022) Genetic diversity and structure of a recent fish invasion: Tench ( Tinca tinca ) in eastern North America. Evol Appl 16:173–188 Gao K, He Z, Xiong J (2024) Population structure and adaptability analysis of Schizothorax o'connori based on whole-genome resequencing. BMC Genomics 25:145 Kon T, Pei L, Ichikawa R (2021) Whole-genome resequencing of large yellow croaker ( Larimichthys crocea ) reveals the population structure and signatures of environmental adaptation. Sci Rep 11:11235 Zhang G, Fang X, Guo X (2012) The oyster genome reveals stress adaptation and complexity of shell formation. Nature 490:49–54 Goh KM, González-Siso MI, Sani RK (2023) Genomics of extreme environments: unveiling the secrets of survival. Sci Rep 13:21441 Zhang X, Yan X (2014) Spatiotemporal change in geographical distribution of global climate types in the context of climate warming. Clim Dyn 43:595–605 Yang D, Li W, Xiang P, Ge T, Li H, Zhang Y (2024) Rhein promotes skin wound healing by activating the PI3K/AKT signaling pathway. Open Med 19:20241116 Martini M, De Santis MC, Braccini L, Gulluni F, Hirsch E (2014) PI3K/AKT signaling pathway and cancer: an updated review. Ann Med 46:372–383 Shin JM, Kim MY, Sohn KC, Jung SY, Lee HE, Lim JW et al. (2014) Nrf2 negatively regulates melanogenesis by modulating PI3K/Akt signaling. PLoS One 9:e96035 Tang H, Liu J, Wang Z, Zhang L, Yang M, Huang J et al. (2023) Genome-wide association study (GWAS) analysis of black color trait in the leopard coral grouper ( Plectropomus leopardus ) using whole genome resequencing. Comp Biochem Physiol Part D Genomics Proteomics 48:101138 Coelho SG, Hearing VJ (2010) UVA tanning is involved in the increased incidence of skin cancers in fair-skinned young women. Pigment Cell Melanoma Res 23:57–63 Qian W, Deng L, Lu D, Xu S (2013) Genome-wide landscapes of human local adaptation in Asia. PLoS One 8:e54224 Nie MM, Ni WL, Wang LH, Gao Q, Liu D, Tian F et al. (2022) Insights into miRNA-mRNA regulatory mechanisms of cold adaptation in Gymnocypris eckloni : ubiquitin-mediated proteolysis is pivotal for adaptive energy metabolism. Front Genet 13:903995 Zhou X, Cao J, Xie J, Tong W, Jia B, Fu J (2025) Effect of Dync1h1 on phototransduction protein transport and the development and maintenance of photoreceptor cells in zebrafish. Invest Ophthalmol Vis Sci 66:38 Haase J, Jones AKC, McVeigh CJ, Brown E, Clarke G, Ahnert-Hilger G (2021) Sex and brain region-specific regulation of serotonin transporter activity in synaptosomes in guanine nucleotide-binding protein G(q) alpha knockout mice. J Neurochem 159:156–171 Wang C, Yang L, Xiao TY, Li JH, Liu QL, Xiong ST (2022) Identification and expression analysis of zebrafish gnaq in the hypothalamic-pituitary-gonadal axis. Front Genet 13:1015796 Pridgeon JW, Klesius PH (2013) G-protein coupled receptor 18 (GPR18) in channel catfish: expression analysis and efficacy as immunostimulant against Aeromonas hydrophila infection. Fish Shellfish Immunol 35:1070–1078 Kiryu I, Köllner B, Kuroda A, Ototake M, Dijkstra JM (2003) A new putative G-protein coupled receptor gene associated with the immune system of rainbow trout ( Oncorhynchus mykiss ). Fish Shellfish Immunol 15:117–127 Additional Declarations There is a duality of interest Supplementary Files SupplementaryMaterial.docx Supplementary Material Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-6768090","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":470006575,"identity":"54ba98c5-765c-462b-8071-bfe145691fda","order_by":0,"name":"Na Song","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAzklEQVRIiWNgGAWjYPACGznG9sbGhx9I0JJmzNxzuNlYggQthxLbZ6S3CfAQo1a+vffw64KaA4y9Mx+2MUgw2MnpNhDQYnDmXJr1jGN3mCVnJ7Y9KGBINjY7QEiLRI6ZMQ/bMzbD2YntBhIMBxK3EdIiPwOk5d9hHvubB9skeIjRwnAjx/gxb9thCcYZjERqMThzxoyZty/NgLEnERjIBkT4Rb69x/gzzzeb+sb24w8ffqiwkyOoBQjYkCLQgLByEGAmJZmMglEwCkbBSAQANQNEMNPMVSIAAAAASUVORK5CYII=","orcid":"","institution":"Ocean University of China","correspondingAuthor":true,"prefix":"","firstName":"Na","middleName":"","lastName":"Song","suffix":""},{"id":470006576,"identity":"48ee382b-529f-4299-ab57-0e1fb98ff5f9","order_by":1,"name":"Ruizhi Wang","email":"","orcid":"","institution":"Ocean University of China","correspondingAuthor":false,"prefix":"","firstName":"Ruizhi","middleName":"","lastName":"Wang","suffix":""},{"id":470006577,"identity":"a4e10267-4d2f-49c7-9779-8ba085fd3118","order_by":2,"name":"Zengman Wu","email":"","orcid":"","institution":"Ocean University of China","correspondingAuthor":false,"prefix":"","firstName":"Zengman","middleName":"","lastName":"Wu","suffix":""},{"id":470006578,"identity":"0b520d49-ae62-4c6d-aa0d-70398391e063","order_by":3,"name":"Bailin Cong","email":"","orcid":"","institution":"First Institute of Oceanography, Ministry of Natural Resources","correspondingAuthor":false,"prefix":"","firstName":"Bailin","middleName":"","lastName":"Cong","suffix":""},{"id":470006579,"identity":"9ae4e862-ff4d-4a68-bdcb-a061b66d7009","order_by":4,"name":"Guangxun Du","email":"","orcid":"","institution":"Ministry of Natural Resources","correspondingAuthor":false,"prefix":"","firstName":"Guangxun","middleName":"","lastName":"Du","suffix":""},{"id":470006580,"identity":"9d16c937-a3d0-472c-b21b-5d86ac7e05bc","order_by":5,"name":"Dongdong Xiao","email":"","orcid":"","institution":"Ocean University of China","correspondingAuthor":false,"prefix":"","firstName":"Dongdong","middleName":"","lastName":"Xiao","suffix":""},{"id":470006581,"identity":"cab03bee-71d3-4ac6-a208-749a1dba41e5","order_by":6,"name":"Linlin Zhao","email":"","orcid":"","institution":"First Institute of Oceanography, Ministry of Natural Resources","correspondingAuthor":false,"prefix":"","firstName":"Linlin","middleName":"","lastName":"Zhao","suffix":""}],"badges":[],"createdAt":"2025-05-28 12:31:20","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-6768090/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-6768090/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":84692199,"identity":"edee8f28-331c-451f-95ae-87c5148282e8","added_by":"auto","created_at":"2025-06-16 10:00:35","extension":"jpeg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":71246,"visible":true,"origin":"","legend":"\u003cp\u003eSchematic of five sampling locations for \u003cem\u003eB. pectinirostris\u003c/em\u003e\u003c/p\u003e","description":"","filename":"image1.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/b5a40494428080c991d2745b.jpeg"},{"id":84692200,"identity":"4b5fe754-0ba8-4616-b7c5-40417e6dbff6","added_by":"auto","created_at":"2025-06-16 10:00:35","extension":"jpeg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":107932,"visible":true,"origin":"","legend":"\u003cp\u003eThe distribution of SNPs across \u003cem\u003eB. pectinirostris\u003c/em\u003e. Chromosome length is displayed on the x-axis, with darker areas indicate regions where SNPs are concentrated\u003c/p\u003e","description":"","filename":"image2.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/be699ba8f9ca1a7fdd26d3fb.jpeg"},{"id":84692198,"identity":"89d276f0-d59a-4d11-ba4d-bd96df7eb13c","added_by":"auto","created_at":"2025-06-16 10:00:35","extension":"jpeg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":78772,"visible":true,"origin":"","legend":"\u003cp\u003e(\u003cstrong\u003eA\u003c/strong\u003e) Phylogenetic tree analysis of the 68 samples. Each colour represents a sampling area; (\u003cstrong\u003eB\u003c/strong\u003e) PCA three-dimensional cluster diagram of the 68 samples showing PC1, PC2, and PC3. Each point represents a sample, and points of different color to different sample groups; (\u003cstrong\u003eC\u003c/strong\u003e) Cluster graphs of K=2-6 with 68 samples. The horizontal axis represents the sample groups. Each color in the figure represents a cluster, with each row corresponding to a cluster value. The CV error was minimized when K=2\u003c/p\u003e","description":"","filename":"image3.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/40a3cfd8d1dc28c29a7ab48e.jpeg"},{"id":84693005,"identity":"a52382c7-3404-4b26-8744-c28e0511e372","added_by":"auto","created_at":"2025-06-16 10:08:35","extension":"jpeg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":23904,"visible":true,"origin":"","legend":"\u003cp\u003eGene flow analysis of five \u003cem\u003eB. pectinirostris\u003c/em\u003e populations (three gene migrations)\u003c/p\u003e","description":"","filename":"image4.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/8de51939ae8f106df8929221.jpeg"},{"id":84692201,"identity":"743f362b-d9a7-4362-956e-35ed446aadfb","added_by":"auto","created_at":"2025-06-16 10:00:35","extension":"jpeg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":30289,"visible":true,"origin":"","legend":"\u003cp\u003eInkage disequilibrium analysis five \u003cem\u003eB. pectinirostris\u003c/em\u003e populations\u003c/p\u003e","description":"","filename":"image5.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/3b58b4fa0c0450fca2a05c2c.jpeg"},{"id":84693006,"identity":"a3f880d7-229e-4ee5-a188-f651433a0921","added_by":"auto","created_at":"2025-06-16 10:08:35","extension":"jpeg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":58619,"visible":true,"origin":"","legend":"\u003cp\u003eDemographic history of five \u003cem\u003eB. pectinirostris\u003c/em\u003e populations\u003c/p\u003e","description":"","filename":"image6.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/2f210567ef5753d762f379df.jpeg"},{"id":84692202,"identity":"81c1b47c-ee76-4f09-9f69-885067d24bfd","added_by":"auto","created_at":"2025-06-16 10:00:35","extension":"jpeg","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":136781,"visible":true,"origin":"","legend":"\u003cp\u003eGenomic regions with strong selective signals in populations of\u0026nbsp;\u003cem\u003eB. pectinirostris\u003c/em\u003e.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003eA\u003c/strong\u003e) Distribution of log\u003csub\u003e2 \u003c/sub\u003e(\u003cem\u003eπ\u003c/em\u003e ratios) and\u0026nbsp;\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e values calculated in 100\u0026nbsp;kb sliding windows with 10\u0026nbsp;kb increments between QD/GY populations; (\u003cstrong\u003eB\u003c/strong\u003e) Distribution of log\u003csub\u003e2 \u003c/sub\u003e(\u003cem\u003eπ\u003c/em\u003e ratios) and\u0026nbsp;\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e values between QD/WZ populations; (\u003cstrong\u003eC\u003c/strong\u003e) Distribution of log\u003csub\u003e2 \u003c/sub\u003e(\u003cem\u003eπ\u003c/em\u003e ratios) and\u0026nbsp;\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e\u0026nbsp;values between QD/XM populations; (\u003cstrong\u003eD\u003c/strong\u003e) Distribution of log\u003csub\u003e2 \u003c/sub\u003e(\u003cem\u003eπ\u003c/em\u003e ratios) and\u0026nbsp;\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e\u0026nbsp;values between QD/ZJ populations\u003c/p\u003e","description":"","filename":"image7.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/1c707251876673dd76e290a2.jpeg"},{"id":84693578,"identity":"1ca43d62-0e79-4a23-8d57-6d1f7ab75216","added_by":"auto","created_at":"2025-06-16 10:16:35","extension":"jpeg","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":63168,"visible":true,"origin":"","legend":"\u003cp\u003e(\u003cstrong\u003eA\u003c/strong\u003e)\u003cstrong\u003e \u003c/strong\u003eGO enrichment for candidate genes; (\u003cstrong\u003eB\u003c/strong\u003e) KEGG enrichment for candidate genes\u003c/p\u003e","description":"","filename":"image8.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/5a77edb6b466350d401bba2f.jpeg"},{"id":84692207,"identity":"fc610d68-1bc1-43cb-ac7e-ceadf1cdf2bc","added_by":"auto","created_at":"2025-06-16 10:00:35","extension":"jpeg","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":53886,"visible":true,"origin":"","legend":"\u003cp\u003e(\u003cstrong\u003eA\u003c/strong\u003e) KEGG enrichment for the selected genes based on XP-EHH values between the QD and WZ populations; (\u003cstrong\u003eB\u003c/strong\u003e) KEGG enrichment for the selected genes based on XP-EHH values between the QD and ZJ populations\u003c/p\u003e","description":"","filename":"image9.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/8a7edfe9d90d95c8583cd635.jpeg"},{"id":91983242,"identity":"c683fdde-5c8d-4230-bf02-84fce1d55081","added_by":"auto","created_at":"2025-09-23 11:28:45","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1638609,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/5abc2e62-3300-467b-ad22-972b0e2a55eb.pdf"},{"id":84692204,"identity":"ed9a8711-6d0f-481d-92cc-1f8fc0d63a56","added_by":"auto","created_at":"2025-06-16 10:00:35","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":322966,"visible":true,"origin":"","legend":"Supplementary Material","description":"","filename":"SupplementaryMaterial.docx","url":"https://assets-eu.researchsquare.com/files/rs-6768090/v1/e0cf40d13307333c4d667cf2.docx"}],"financialInterests":"There is a duality of interest","formattedTitle":"Evolutionary Clues of North-South Divergence: Genomic Adaptation of the Intertidal Mudskipper Boleophthalmus pectinirostris","fulltext":[{"header":"Introduction","content":"\u003cp\u003e\u003cem\u003eBoleophthalmus pectinirostris\u003c/em\u003e, a member of the family Gobiidae, is a distinctive amphibious fish species highly adapted to intertidal mudflat ecosystems. It is widely distributed along the western Pacific coast, spanning regions such as China, Japan, Korea, Malaysia, and Indonesia\u0026nbsp;(Murdy 2017). As a warm-temperate species, \u003cem\u003eB. pectinirostris\u003c/em\u003e thrives in coastal intertidal mudflats and demonstrates exceptional tolerance to a wide range of salinities (Takegaki 2008). \u003cem\u003eB. pectinirostris\u003c/em\u003e has developed a unique feeding strategy, primarily consuming benthic diatoms. This dietary specialization further influences its distribution in aquatic environments (Yang et al. 2003; Polgar et al. 2009). Due to its relatively low trophic position and short food chain, \u003cem\u003eB. pectinirostris\u003c/em\u003e can rapidly respond to environmental changes, making it a valuable bioindicator for coastal ecosystem monitoring (Pan et al. 2022). Moreover, it has been recognized as a key species for assessing the ecological health of mudflat and mangrove ecosystems (Polgar 2009).\u003c/p\u003e\n\u003cp\u003eDespite growing interest in the genetic characteristics of \u003cem\u003eB. pectinirostris\u003c/em\u003e, previous studies have predominantly relied on traditional molecular markers, such as mitochondrial DNA sequences and microsatellites, to explore its genetic diversity, phylogenetic relationships, and population differentiation\u0026nbsp;(Kanemori et al. 2006; Gong et al. 2017). However, these methods often suffer from limited resolution and insufficient sample coverage, hindering a comprehensive understanding of population genetic structure and adaptive evolution. In contrast, whole-genome resequencing (WGS) offers a powerful approach for analyzing genomic variations, including single nucleotide polymorphisms (SNPs) and structural variants, by aligning sequencing reads to a reference genome (Wang et al. 2024). With the rapid advancement of sequencing technologies, an increasing number of marine fish species have been fully sequenced, and WGS has emerged as a critical tool for elucidating population genetic diversity and local adaptation (Xiong et al. 2023). Furthermore, downstream analyses enable the identification of candidate genes and pathways associated with environmental adaptation, providing insights into the molecular mechanisms underlying ecological resilience in fish species (Li et al. 2023; Sun et al. 2023; Huang et al. 2024).\u003c/p\u003e\n\u003cp\u003eIntertidal organisms such as \u003cem\u003eB. pectinirostris\u0026nbsp;\u003c/em\u003erespond to extreme environmental challenges, including temperature changes caused by Marine heat waves\u0026nbsp;(Wang et al. 2020;\u0026nbsp;Stillman\u0026nbsp;et al. 2025). These selective pressures are likely to drive the accumulation of adaptive genetic variations through polygenic selection (You et al. 2014). Previous whole-genome sequencing studies on the intertidal snail Littorina brevicula have revealed two distinct genetic lineages along a latitudinal thermal gradient, with north-south divergence driven by thermal environmental variation\u0026nbsp;(Lin et al., 2024).\u0026nbsp;Similarly, mudskippers, as emblematic intertidal organisms exhibiting amphibious adaptations, demonstrate heightened sensitivity to environmental fluctuations, making them an ideal model for studying adaptive evolution in dynamic habitats.\u003c/p\u003e\n\u003cp\u003ePrevious studies based on traditional genetic markers (e.g., mitochondrial DNA and microsatellites) have shown no significant population differentiation within the southern Chinese populations, but revealed two distinct lineages between southern Chinese populations and those from Japan and South Korea\u0026nbsp;(Wang et al. 2012; Chen et al. 2015). However, the limited resolution of these markers, combined with the lack of research on northern Chinese populations, has hindered a comprehensive understanding of population divergence and the genetic basis of their adaptive potential. Given the geographical proximity of Japanese and South Korean populations to northern China, we hypothesized that significant genetic structure may exist between northern and southern Chinese populations of \u003cem\u003eB. pectinirostris\u003c/em\u003e, potentially reflecting distinct adaptive evolutionary responses to local environmental conditions. To address these gaps, in this study, we employed whole-genome resequencing to analyze populations of \u003cem\u003eB. pectinirostris\u003c/em\u003e from five coastal locations, including a northern population from Qingdao, China. Using population genomic data, this study aims to unravel the demographic history, genetic diversity patterns, and adaptive genetic mechanisms of \u003cem\u003eB. pectinirostris\u003c/em\u003e. By analyzing genome-wide variations, we can not only more accurately characterize population structure and divergence history but also identify key genomic regions associated with local adaptation. These findings will provide new insights into the adaptive evolution of intertidal species and offer a scientific basis for the conservation and management of this ecologically significant species.\u003c/p\u003e"},{"header":"Materials And Methods","content":"\u003cp\u003e\u003cstrong\u003eSample collection, DNA extraction, and genome resequencing\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eA total of 68 \u003cem\u003eB. pectinirostris\u003c/em\u003e specimens were collected from five coastal regions in China: Qingdao (QD) in Shandong Province, Ganyu (GY) in Jiangsu Province, Wenzhou (WZ) in Zhejiang Province, Xiamen (XM) in Fujian Province, and Zhanjiang (ZJ) in Guangdong Province (Fig. 1).\u0026nbsp;For each region, 10\u0026ndash;15 specimens were collected and identified through morphological examination and COⅠ gene analysis before being used for subsequent whole-genome studies.\u0026nbsp;Dorsal muscle tissue was dissected from each specimen, preserved in absolute ethanol, and stored at -80\u0026deg;C in an ultra-low-temperature freezer until genomic DNA extraction.\u003c/p\u003e\n\u003cp\u003eGenomic DNA was extracted using the standard phenol-chloroform method. The extracted DNA was randomly fragmented using an ultrasonic disruptor to construct sequencing libraries. The library preparation process included end repair, poly(A) tailing, adaptor ligation, and PCR amplification. The quality and concentration of the genomic libraries were initially assessed using a Qubit 2.0 fluorometer, followed by verification of insert fragment sizes using an Agilent 2100 Bioanalyzer. Finally, paired-end sequencing was performed on the Illumina HiSeq PE150 platform, generating reads with an average length of 150 bp. All experimental procedures were conducted at Novogene Bioinformatics Technology Co., Ltd.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eSequencing data quality control and sequence alignment\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe raw sequencing data contained adaptor sequences and low-quality bases, which could compromise subsequent analyses. To ensure data quality, raw reads were filtered using Fastp software\u0026nbsp;(Chen et al. 2018) to remove adaptors and filter out low-quality sequences, resulting in high-quality clean reads. For this study, the first chromosome-level genome of \u003cem\u003eB. pectinirostris\u003c/em\u003e (GenBank accession number: GCA_026225935.1) was used as the reference genome (Bian et al. 2023). Clean reads from each individual were aligned to the reference genome using BWA (Li and Durbin 2009) with the parameters \u0026ldquo;mem -t 4 -k 32 -M.\u0026rdquo; The alignment results were saved in BAM format. Subsequently, SAMtools (Li et al. 2009) was employed to sort the BAM files, and duplicate reads were removed using Picard v2.18.17.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eSNP detection and annotation\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eBased on the alignment results, variant calling was performed using GATK\u0026nbsp;(McKenna et al. 2010) to generate a VCF file containing raw SNPs. The raw SNPs were filtered using stringent quality control criteria (parameters: QD \u0026lt; 2.0 || MQ \u0026lt; 40.0 || FS \u0026gt; 60.0 || SOR \u0026gt; 3.0 || MQRankSum \u0026lt; -12.5 || ReadPosRankSum \u0026lt; -8.0). High-quality SNPs retained after filtering were annotated using ANNOVAR software (Wang et al. 2010). The SNPs were categorized based on their genomic locations (e.g., upstream, downstream, intronic, intergenic, and splicing regions) and functional effects (e.g., synonymous and non-synonymous variants) (Jones et al. 2014).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePopulation genetic structure and\u0026nbsp;\u003c/strong\u003e\u003cstrong\u003egenetic diversity analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAfter filtering, a set of high-quality SNPs was obtained for subsequent population genetic analyses. Three methods were employed to investigate population genomic relationships: Neighbor-Joining (NJ) tree construction, Principal Component Analysis (PCA) and ADMIXTURE analysis. First, a genetic distance matrix was calculated using Phylip software, and a phylogenetic tree was constructed using the Neighbor-Joining (NJ) method\u0026nbsp;(Felsenstein 1993). The resulting tree was visualized using iTOL (https://itol.embl.de/). PCA was conducted on whole-genome SNPs from all 68 individuals using PLINK v1.9\u0026nbsp;(Purcell et al. 2007). Additionally, population genetic structure and lineage information were inferred using ADMIXTURE v1.3 software, with K values ranging from 2 to 6 (Alexander et al. 2009).\u003c/p\u003e\n\u003cp\u003eTo evaluate patterns of linkage disequilibrium (LD) decay across the genomes of the five populations, allele frequency correlation (r\u0026sup2;) was calculated using PopLDdecay v3.42\u0026nbsp;(Zhang et al. 2019). To further explore genetic variation between northern and southern populations of \u003cem\u003eB. pectinirostris\u003c/em\u003e, we assessed the genome-wide distribution of nucleotide diversity (\u003cem\u003e\u0026pi;\u003c/em\u003e) and the genetic differentiation index (\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e) using VCFtools v0.1.16. These analyses were performed with a sliding window approach, using a window size of 100 kb and a step size of 10 kb\u0026nbsp;(Danecek et al., 2011).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEffective population size analysis of evolutionary history\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo reconstruct the demographic history of different populations, Paired Sequentially Markovian Coalescent (PSMC) analysis was employed to infer changes in the effective population size of \u003cem\u003eB. pectinirostris\u003c/em\u003e over its evolutionary history\u0026nbsp;(Li et al. 2011). Input files for PSMC modeling were generated using the \u0026ldquo;fq2psmcfa\u0026rdquo; and \u0026ldquo;splitfa\u0026rdquo; tools within the PSMC software. The PSMC analysis was conducted with the following parameters: \u0026ldquo;-N25\u0026rdquo; for the number of algorithm iterations, \u0026ldquo;-t15\u0026rdquo; as the upper limit for the time to the most recent common ancestor (TMRCA), \u0026ldquo;-r5\u0026rdquo; for the initial \u0026theta;/\u0026rho; ratio, and \u0026ldquo;-p 4\u0026thinsp;+\u0026thinsp;25 * 2\u0026thinsp;+\u0026thinsp;4\u0026thinsp;+\u0026thinsp;6\u0026rdquo; for atomic time intervals. The reconstructed population history was visualized using the \u0026ldquo;psmc_plot.pl\u0026rdquo; script. A mutation rate of 3.51\u0026times;10⁻⁹ per generation and a generation time of 1 year were assumed for the analysis (You et al. 2014). Additionally, TreeMix v1.13 was used to infer the maximum-likelihood (ML) tree for the five populations. A window size of 1000 was applied to account for linkage disequilibrium (\u0026ndash;k), and the \u0026ldquo;\u0026ndash;global\u0026rdquo; option was used to generate the ML tree. Migration events were sequentially added to the tree by testing models with 0 to 4 migration events (\u0026ndash;m 0 1 2 3 4) (Pickrell et al. 2012).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eSelection signature detection and functional enrichment analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe \u003cem\u003e\u0026pi;\u003c/em\u003e ratio between the two populations was calculated based on the previously estimated nucleotide diversity (\u003cem\u003e\u0026pi;\u003c/em\u003e). The results of the population differentiation coefficient (\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e) and \u003cem\u003e\u0026pi;\u003c/em\u003e ratio were then integrated and summarized\u0026nbsp;(Danecek et al. 2011). The northern population (QD) was compared collectively against the four southern populations (GY, WZ, XM, and ZJ) in the analysis. Genomic windows with values in the top 5% for both \u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e and \u003cem\u003e\u0026pi;\u003c/em\u003e ratio were identified as candidate regions under strong selective pressure. Based on these candidate regions, corresponding SNPs and genes were annotated. Additionally, we performed an XP-EHH (cross population extended haplotype homozogysity)\u0026nbsp;analysis using SELSCAN v1.3.0 with default settings to detect ongoing or nearly fixed selective sweeps by comparing haplotype patterns across the five populations\u0026nbsp;(Szpiech et al. 2014). For the XP-EHH selection scan, genes located within 10 kb upstream and downstream of the identified loci were extracted for further analysis.\u003c/p\u003e\n\u003cp\u003eFor genes exhibiting selection signals in each population, functional annotation was performed using eggNOG-mapper software (Cantalapiedra et al. 2021) to assign Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. GO term classification and KEGG pathway enrichment analysis were further conducted using OmicShare tools (http://www.omicshare.com/tools) (Ashburner et al. 2000; Jones et al. 2014; Kanehisa et al. 2023).\u003c/p\u003e"},{"header":"Results","content":"\u003cp\u003e\u003cstrong\u003eSequencing and SNP calling\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe study conducted library construction and sequencing on 68 samples from five geographic populations of \u003cem\u003eB. pectinirostris\u003c/em\u003e. A total of approximately 821.72 Gb of raw sequencing data was generated, with an average of 12.08 Gb per sample. The average sequencing depth was 12.54\u0026times;,\u0026nbsp;with an average GC content of\u0026nbsp;40.50%. The average Q20 and Q30 scores were 96.49% and 91.65%, respectively (Fig. S1). After quality filtering, approximately 805.11 Gb of high-quality clean data was retained for downstream analyses.\u003c/p\u003e\n\u003cp\u003eBy aligning the sequencing data with the reference genome, a total of 33,016,712 SNP variants were identified. The distribution of these SNPs provided valuable insights into the genomic variation across the population. The majority of the variants (63.6%) were located in intergenic regions, reflecting the extensive noncoding regions of the genome. Intronic variants accounted for 18.7%, while coding variants, including 200,311 synonymous and 490,528 non-synonymous mutations, were less frequent. Regulatory regions, encompassing both upstream and downstream areas, represented 15.5% of the SNPs, suggesting their potential role in gene expression regulation. Variants in splicing-related regions and untranslated regions (UTRs) were rare. To further validate the reliability of the data, transition/transversion (Ts/Tv) analysis will be conducted (Table 1). These findings provide a comprehensive overview of the SNP landscape within this population. Additionally, we analyzed the distribution of core SNP markers across the chromosomes. The SNP density distribution map revealed that SNPs are evenly distributed across all chromosomes (Fig. 2).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eTable\u0026nbsp;\u003c/strong\u003e\u003cstrong\u003e1\u003c/strong\u003e The distribution of SNP variants in the genome region.\u003c/p\u003e\n\u003ctable border=\"1\" cellspacing=\"0\" cellpadding=\"0\"\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eGenetic element\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eSNP numbers\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eUpstream\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e3,068,552\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eDownstream\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e3,061,303\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eCDS\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e753,445\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eIntron\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e7,335,438\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eIntergenic\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e25,129,719\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eSplicing\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e53,847\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e5\u0026apos; UTR\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e233\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eSynonymous\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e200,311\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eNon-synonymous\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e490,528\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eTransitions (Ts)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e17,040,364\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eTransversions (Tv)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e11,466,334\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003eTs/Tv Ratio\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 284px;\"\u003e\n \u003cp\u003e1.49\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n\u003c/table\u003e\n\u003cp\u003e\u003cbr\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePopulation genetic structure analysis and gene flow\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo investigate the genetic structure of \u003cem\u003eB. pectinirostris\u003c/em\u003e across different populations, phylogenetic analysis was conducted using high-quality SNP loci filtered from the sequencing data of samples collected from five geographic populations. The resulting phylogenetic tree revealed that populations generally clustered according to their geographic origins, with the exception of two individuals from the ZJ population, which grouped with the WZ population (Fig. 3A). A three-dimensional PCA plot was constructed based on the first (PCA1), second (PCA2), and third (PCA3) principal components. The PCA results indicated that the XM and ZJ populations clustered together, suggesting a closer genetic relationship and higher genetic similarity between these two populations. In contrast, the QD and GY populations formed distinct clusters, reflecting greater genetic divergence from the other populations. The WZ population occupied an intermediate position between the other groups, consistent with its geographic location (Fig. 3B). The Admixture analysis further supported the results of the phylogenetic tree and PCA. When K=2, the cross-validation (CV) error was minimized, and the populations could be broadly divided into two major clusters: Qingdao formed one cluster, while XM and ZJ clustered together. The GY and WZ populations were positioned between these two groups. As K increased, finer genetic distinctions between the clusters were observed. For K values between 3 and 6, QD remained a distinct cluster, while the WZ population showed a genetic composition more similar to the XM and ZJ populations, indicating closer genetic relatedness between them (Fig. 3C). We also calculated the \u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e between populations based on all SNPs and then averaged them. Both the other four populations showed significant genetic divergence from the QD population (Fig. S2). This comprehensive genetic analysis provides valuable insights into the population structure and evolutionary relationships of \u003cem\u003eB. pectinirostris\u003c/em\u003e.\u003c/p\u003e\n\u003cp\u003eBased on the results of genetic differentiation and cluster analysis, the gene flow was performed. Gene flow analysis showed that when three migration events occurred, the QD population gene flowed to GY, and GY population gene flowed to WZ, which was similar to the results of cluster analysis. This implies that the QD and GY populations perhaps originated from a single ancestor. It is worth noting that gene flow from GY to ZJ was detected despite their geographic separation, possibly due to anthropogenic influences (Fig. 4).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eLinkage disequilibrium analysis and history dynamics\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eBased on the results of genetic differentiation and cluster analysis, the linkage disequilibrium analysis was performed. LD analysis was conducted for five populations. The attenuation distance of the four southern populations was smaller than that of the Qingdao population, which may be due to the high genetic diversity of these populations. This was consistent with the results of genetic diversity analysis (Fig. 5).\u003c/p\u003e\n\u003cp\u003eBased on the PSMC model, we reconstructed the changes of the historical effective population size (\u003cem\u003eN\u003c/em\u003e\u003csub\u003ee\u003c/sub\u003e) of four different geographic groups of \u003cem\u003eB. pectinirostris\u003c/em\u003e to infer the population evolution history. The results showed that the \u003cem\u003eN\u003c/em\u003e\u003csub\u003ee\u003c/sub\u003e changes of different populations were basically consistent about 150,000 years ago, and they all experienced a population contraction about 700,000 to 800,000 years ago. After that, the \u003cem\u003eN\u003c/em\u003e\u003csub\u003ee\u003c/sub\u003e of QD and GY populations in the north was slightly larger than that of other populations in the south. In the period after 150,000 years ago, the trajectory of \u003cem\u003eN\u003c/em\u003e\u003csub\u003ee\u003c/sub\u003e change of QD populations and GY populations deviated from that of the other three populations, and the historical effective population size of the other three groups showed a downward trend, while that of the other three populations showed an upward trend. In general, a clear latitudinal gradient in effective population size (\u003cem\u003eN\u003c/em\u003e\u003csub\u003ee\u003c/sub\u003e) was evident by 10,000 years ago, with southern populations maintaining higher values than their northern counterparts, and the QD population exhibiting the lowest \u003cem\u003eN\u003c/em\u003e\u003csub\u003ee\u003c/sub\u003e. (Fig. 6).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eGenome-wide selective signals and candidate gene analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eIn this study, we identified genomic regions under selection using \u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e \u0026amp; \u003cem\u003e\u0026pi;\u003c/em\u003e to uncover genes associated with environmental adaptation. Considering that there were no significant genetic differences between the southern populations, we selected the southern populations as control groups to reveal environmental selection genes in the northern population. With four southern populations as the control groups, we screened 1685 (GY, z (\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e)\u0026nbsp;\u0026ge;\u0026nbsp;1.84 and log\u003csub\u003e2\u003c/sub\u003e (\u003cem\u003e\u0026pi;\u0026nbsp;\u003c/em\u003eQD/GY)\u0026nbsp;\u0026ge;\u0026nbsp;1.36), 1964 (WZ, z (\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e)\u0026nbsp;\u0026ge;\u0026nbsp;1.86 and log\u003csub\u003e2\u003c/sub\u003e (\u003cem\u003e\u0026pi;\u0026nbsp;\u003c/em\u003eQD/GY)\u0026nbsp;\u0026ge;\u0026nbsp;1.50), 1859 (z (\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e)\u0026nbsp;\u0026ge;\u0026nbsp;1.81 and log\u003csub\u003e2\u003c/sub\u003e (\u003cem\u003e\u0026pi;\u0026nbsp;\u003c/em\u003eQD/GY)\u0026nbsp;\u0026ge;\u0026nbsp;1.47) and 1892 (ZJ, There were z (\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e)\u0026nbsp;\u0026ge;\u0026nbsp;1.84 and log\u003csub\u003e2\u003c/sub\u003e (\u003cem\u003e\u0026pi;\u0026nbsp;\u003c/em\u003eQD/GY)\u0026nbsp;\u0026ge;\u0026nbsp;1.54) candidate genes (Fig. 7). Therefore, we identified 977 overlapping genes from the QD population screened with the southern population as the control group as potential genes related to adaptation to the northern environment (Fig. S3A). Similarly, with QD population as the control group, we obtained 381, 98, 273 and 109 candidate genes in GY, WZ, XM and ZJ populations (Fig. 7). A total of 14 overlapping genes were identified as potential genes associated with southern environmental adaptation (Fig. S3B).\u003c/p\u003e\n\u003cp\u003eIn order to fully understand the molecular functions of these candidate genes related to environmental adaptation, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for these candidate genes. Since the southern population exhibited fewer overlapping genes, all subsequent analyses primarily focused on the candidate genes from the northern population. The GO analysis results indicated that the major GO terms associated with adaptation to the northern environments included \u003cstrong\u003eCarboxylic acid biosynthetic process\u003c/strong\u003e (GO:0046394) and \u003cstrong\u003eOrganic acid biosynthetic process\u003c/strong\u003e (GO:0016053) in Biological Processes, as well as \u003cstrong\u003eATP binding\u003c/strong\u003e (GO:0005524), \u003cstrong\u003eAdenyl nucleotide binding\u003c/strong\u003e (GO:0030554), and \u003cstrong\u003eAdenyl ribonucleotide binding\u003c/strong\u003e (GO:0032559) in Molecular Function (Fig. 8A). The KEGG enrichment results revealed that the candidate genes were significantly enriched in several metabolic pathways, including the \u003cstrong\u003ePI3K-Akt signaling pathway\u003c/strong\u003e (ko04151), \u003cstrong\u003eTight junction\u003c/strong\u003e (ko04530), \u003cstrong\u003eRap1 signaling pathway\u003c/strong\u003e (ko04015), \u003cstrong\u003ePlatelet activation\u003c/strong\u003e (ko04611), \u003cstrong\u003eRegulation of actin cytoskeleton\u003c/strong\u003e (ko04810), and \u003cstrong\u003eFocal adhesion\u003c/strong\u003e (ko04510) (Fig. 8B).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eGenome-wide selection signal and candidate gene analysis based on XP-EHH\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eBased on the XP-EHH values at each locus, we selected the top 1% regions as candidate regions. After filtering, we detected strong selection signals between the QD population and the other four populations. A total of 105 genes were identified between the QD and GY populations (Fig. S4A). Meanwhile, 314, 312, and 321 genes were identified between the QD and WZ, XM, and ZJ populations, respectively (Fig. S4B; Fig. S4C; Fig. S4D). KEGG enrichment analysis was performed on the genes screened from these four comparisons. The results showed that the genes between the QD and WZ populations, as well as the QD and ZJ populations, were significantly enriched in pathways such as Hepatitis C (ko05160), Two-component system (ko02020), Phototransduction-fly (ko04745), Huntington disease (ko05016), and Tight junction (ko04530) (Fig. 9). Among these, the Hepatitis C and Tight junction pathways were consistent with the \u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e \u0026amp; \u003cem\u003e\u0026pi;\u003c/em\u003e analysis. Additionally, the genes between the QD and GY populations, as well as the QD and XM populations, were significantly enriched in pathways such as Two-component system (ko02020), Phototransduction-fly (ko04745), Amoebiasis (k0), and Nitrogen metabolism (ko00910) (Fig. S5).\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003e\u003cstrong\u003eGenome-wide SNPs decipher phylogeographic structure\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eSingle nucleotide polymorphisms (SNPs) have been well established as a major source of genetic variation in model organisms\u0026nbsp;(Langley et al. 2012). In the present study, we identified 33,016,712 high-quality SNPs, confirming their prevalence as important molecular markers. These genomic resources enable comprehensive investigations of genetic diversity, population structure, and environmental adaptation in natural populations. Such analyses provide critical insights for multiple applications, including germplasm resource conservation, phylogenetic reconstruction, and evolutionary studies (Zhou et al. 2023).\u003c/p\u003e\n\u003cp\u003eUsing three complementary methods\u0026mdash;PCA, ADMIXTURE, and NJ tree analysis\u0026mdash;we identified clear geographic clustering patterns among the five populations. Specifically, the QD population formed a distinct cluster, showing limited gene flow with the geographically proximate GY population. In contrast, the ZJ and XM populations exhibited closer genetic relationships, while the WZ population occupied an intermediate genetic position between the northern and southern groups. Whole-genome SNP analysis revealed significant population genetic structure and differentiation between the QD population in Shandong Province and the other four populations. Additionally, higher\u0026nbsp;\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e values were observed within the southern group, consistent with previous findings based on \u003cem\u003ecytochrome b\u003c/em\u003e analyses\u0026nbsp;(Chen et al. 2015). According to principles of genetics and evolutionary biology, genetic diversity is the foundation for species to adapt to their environments over time and to evolve in response to environmental changes. It is also a fundamental determinant of a species\u0026apos; evolutionary potential (Ma et al. 2022). Our results showed that the ZJ population exhibited the highest nucleotide diversity, while the QD population had the lowest. Previous studies have also indicated that \u003cem\u003eB. pectinirostris\u003c/em\u003e generally exhibits low genetic diversity, with significant variation among populations and a gradual decline in diversity from south to north, which aligns with the findings of this study (Chen et al. 2014).\u003c/p\u003e\n\u003cp\u003eIn this study, the NJ tree revealed that all individuals clustered into three major branches: the southernmost WZ, XM, and ZJ populations formed one branch, the QD population formed another, and the GY population occupied an intermediate branch between the two. This clustering pattern aligns with the geographical distribution of the populations. Previous analyses based on the COI gene have also indicated frequent gene flow among \u003cem\u003eB. pectinirostris\u003c/em\u003e populations across several southern coastal provinces, including Guangdong, Fujian, and Jiangsu, which is consistent with the findings of this study\u0026nbsp;(Yang et al. 2012).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePopulation history and evolutionary dynamics\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe population genetic structure and genetic diversity of species are closely tied to their historical demographic dynamics (McKeown et al. 2019; Wang et al. 2020). In this study, we reconstructed the demographic history of different geographic populations of \u003cem\u003eB. pectinirostris\u003c/em\u003e using the PSMC model. The results revealed that, prior to the early Pleistocene of the Quaternary, the effective population sizes of the various \u003cem\u003eB. pectinirostris\u003c/em\u003e populations were relatively stable. However, during the Glacial-Interglacial Cycle, \u003cem\u003eB. pectinirostris\u003c/em\u003e experienced both population bottlenecks and expansion events. Notably, between 120,000 and 150,000 years ago, the demographic trajectories of the populations began to diverge. During this period, the QD and GY populations started to separate from the other three populations, with the QD population undergoing a more rapid decline in effective population size compared to the GY population. This suggests that genetic differentiation between the QD population and the other four populations likely began around this time. Additionally, we observed that during the Last Glacial Period, extreme environmental conditions led to a significant reduction in the effective population sizes of all populations, a pattern consistent with findings in many other marine species (Wu et al. 2022). Furthermore, over the past 10,000 years, the QD population exhibited a notably smaller effective population size compared to the other populations, which is partially supported by the results of genetic diversity analysis. This may be attributed to the QD population being in a state of continuous decline without experiencing population expansion, resulting in lower genetic diversity (Bernos et al. 2022; Gao et al. 2024).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eStrong selective signals\u0026nbsp;\u003c/strong\u003e\u003cstrong\u003ebased on \u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e\u0026amp;\u003cem\u003e\u0026pi;\u0026nbsp;\u003c/em\u003e\u003c/strong\u003e\u003cstrong\u003eunderlie north-south adaptation divergence\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eRecent studies have demonstrated that whole-genome resequencing technology provides a powerful tool for investigating the environmental adaptation of fish populations across different geographic regions\u0026nbsp;(Kon et al. 2021; Li et al. 2023; Huang et al. 2024). Previous studies based on mitochondrial gene fragments suggested a lack of distinct genetic structure among \u003cem\u003eB. pectinirostris\u003c/em\u003e populations along the Chinese coast (Chen et al. 2015). In contrast, this study revealed significant genetic differentiation among \u003cem\u003eB. pectinirostris\u003c/em\u003e populations in coastal China using whole-genome SNP markers. We hypothesize that this differentiation may result from the species\u0026apos; wide distribution range and limited dispersal ability, making it highly susceptible to local environmental heterogeneity (Zhang et al. 2012).\u003c/p\u003e\n\u003cp\u003eOur genetic structure analyses,\u0026nbsp;\u003cem\u003eF\u003c/em\u003e\u003csub\u003est\u003c/sub\u003e-based population differentiation metrics, and historical demographic inference indicated that the QD population exhibited the most pronounced genetic divergence from the other four populations. Therefore, we selected the QD population as the reference group and compared it with the four southern populations (GY, WZ, XM, and ZJ) to identify genomic regions associated with environmental adaptation. We detected 997 candidate genes in the northern population and 14 in the southern population. The higher number of selected genes in the northern population suggests stronger natural selection pressures, potentially driven by colder temperatures or greater seasonal variability. Populations under intense natural selection often exhibit lower genetic diversity, which aligns with our finding that the QD population had the lowest genetic diversity\u0026nbsp;(Goh et al. 2023). Compared to the warm, humid climate and mangrove-rich habitats in the south, the northern environment is harsher, requiring \u003cem\u003eB. pectinirostris\u003c/em\u003e to employ more diverse strategies to cope with environmental challenges (Zhang et al. 2014).\u003c/p\u003e\n\u003cp\u003eGiven the weaker selection signals in the southern populations, we focused on analyzing the selection signals in the northern population. Functional annotation and clustering of the selected genes revealed their involvement in energy metabolism regulation, DNA damage repair, and signal transduction. Notably, the PI3K-Akt signaling pathway, known for its role in skin wound healing, has not been extensively studied in marine fish\u0026nbsp;(Yang et al. 2024). This pathway is also critical for regulating cellular processes such as mitosis and melanin production (Martini et al. 2014). Activation of the PI3K/Akt signaling pathway has been shown to reduce melanin production by downregulating genes such as TYR and TYRP1 (Shin et al. 2014; Tang et al. 2023). The selection of genes related to skin pigmentation and DNA damage repair may reflect differences in photoperiod and light intensity between northern and southern environments. Additionally, the mangrove ecosystems in the south may provide protection against ultraviolet radiation, a factor absent in the north (Coelho et al. 2010; Qian et al. 2013). The selection of energy metabolism-related genes may be linked to the QD population\u0026apos;s adaptation to higher-latitude environments, where they undergo a 3\u0026ndash;4 month hibernation period in caves during winter. This requires precise regulation of energy metabolism to survive low temperatures and other harsh conditions (Nie et al. 2022). The enrichment of GO terms related to ATP binding further supported this hypothesis. In conclusion, our adaptive evolution analysis identified candidate genes associated with the adaptation of \u003cem\u003eB. pectinirostris\u003c/em\u003e to northern and southern environments. These findings provided valuable insights for further research into the evolutionary history and environmental adaptation mechanisms of this species.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePopulation selection signals based on XP-EHH\u0026nbsp;\u003c/strong\u003e\u003cstrong\u003euncovers vision and immunity-related adaptations\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eIn addition to selective sweep analysis, this study employed the XP-EHH method to detect strong selection signals between the northern and southern populations. A total of 105 to 321 candidate genes were identified between the QD population and the other four populations. These genes were primarily associated with immune response, signal transduction, and DNA damage repair. KEGG enrichment analysis revealed that genes identified in comparisons between the QD and WZ populations, as well as the QD and ZJ populations, were enriched in the Tight Junction pathway. Additionally, genes from all four comparisons were enriched in the Phototransduction-fly pathway, which has been implicated in visual adaptation in fish and may reflect the visual adaptation of Boleophthalmus pectinirostris during terrestrial activities (Zhou et al. 2025). Notably, the Guanine Nucleotide-Binding Protein (GNAQ) gene and the G protein-coupled receptor 18 (GPR18) gene were identified among the candidate genes. GNAQ, a regulator or transducer, plays a role in various transmembrane signaling systems (Haase et al. 2021). In fish, GNAQ has also been linked to olfactory development and immune function (Wang et al. 2022). Similarly, GPR18 has been shown to be closely associated with immune responses in fish, providing protection against certain pathogenic infections (Pridgeon et al. 2013). Furthermore, LPSenhR-1, a member of the same GPR family, has been demonstrated to play a role in immunity in Oncorhynchus mykiss (Kiryu et al. 2003). These findings suggest that these genes may contribute to the adaptation of northern populations to their environment.\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003eThis study employed whole-genome sequencing to elucidate the genetic divergence patterns among five geographic populations of \u003cem\u003eB. pectinirostris\u003c/em\u003e along China\u0026apos;s coastline. The results indicated that there is a pronounced genetic isolation and reduced diversity in the northern Qingdao (QD) population, contrasting with higher gene flow in southern populations. Selective sweep analysis revealed strong signatures of selection in the northern population, particularly in genes associated with energy metabolism and immune response, reflecting adaptive evolution to colder climates and intense UV radiation. In contrast, southern populations exhibited weaker selection signals, likely due to the buffering effect of mangrove habitats. These findings provide compelling evidence that both geographic isolation and environmental gradients drive population differentiation in this ecologically important species. From a conservation perspective, our results highlight the urgent need for prioritized protection of the genetically distinct QD population to prevent further loss of genetic diversity, while maintaining gene flow among southern populations to preserve their adaptive potential. This study offers novel genomic insights into the mechanisms of adaptive evolution in intertidal organisms and establishes a foundation for evidence-based conservation strategies.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgements\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis study was funded by\u0026nbsp;the National Key Research and Development Program of China (2023YFD2401103). Thank all the members who participated in this research for their contributions to this paper\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eResearch\u0026nbsp;\u003c/strong\u003e\u003cstrong\u003eEthics\u0026nbsp;\u003c/strong\u003e\u003cstrong\u003eStatement\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll experimental procedures involving \u003cem\u003eBoleophthalmus pectinirostris\u003c/em\u003e were conducted in accordance with relevant guidelines and regulations.\u0026nbsp;The animal study was reviewed and approved by the Ocean University of China Academic Committee.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData archiving\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe genome resequencing data for\u0026nbsp;\u003cem\u003eBoleophthalmus pectinirostris\u003c/em\u003e have been deposited at NCBI under the BioProject accession number PRJNA1258070.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConflict of interest\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eMurdy EO, Jaafar Z (2017) Taxonomy and systematics review. In: Fishes Out of Water, 1st edn. CRC Press, pp 36\u003c/li\u003e\n\u003cli\u003eTakegaki T (2008) Threatened fishes of the world: Boleophthalmus pectinirostris (Linnaeus 1758) (Gobiidae). Environ Biol Fish 81:373\u0026ndash;374\u003c/li\u003e\n\u003cli\u003eYang KY, Lee SY, Williams GA (2003) Selective feeding by the mudskipper (\u003cem\u003eBoleophthalmus pectinirostris\u003c/em\u003e) on the microalgal assemblage of a tropical mudflat. Mar Biol 143:245\u0026ndash;25\u003c/li\u003e\n\u003cli\u003ePolgar G, Crosa G (2009) Multivariate characterisation of the habitats of seven species of Malayan mudskippers (Gobiidae: Oxudercinae). Mar Biol 156:1475\u0026ndash;1486\u003c/li\u003e\n\u003cli\u003ePan Y, Tian L, Zhao Q, Tao Z, Yang J, Zhou Y et al. (2022) Evaluation of the acute toxic effects of crude oil on intertidal mudskipper (\u003cem\u003eBoleophthalmus pectinirostris\u003c/em\u003e) based on antioxidant enzyme activity and the integrated biomarker response. Environ Pollut 292:118341\u003c/li\u003e\n\u003cli\u003ePolgar G (2009) Species-area relationship and potential role as a biomonitor of mangrove communities of Malayan mudskippers. Wetl Ecol Manag 17:157\u0026ndash;164\u003c/li\u003e\n\u003cli\u003eKanemori Y, Takegaki T, Natsukari Y (2006) Genetic population structure of the mudskipper \u003cem\u003eBoleophthalmus pectinirostris\u003c/em\u003e inferred from mitochondrial DNA sequences. Jpn J Ichthyol 53:133\u0026ndash;141\u003c/li\u003e\n\u003cli\u003eGong W, Huang Y, Xie J (2017) Genome-wide identification of novel microRNAs from genome sequences using computational approach in the mudskipper (\u003cem\u003eBoleophthalmus pectinirostris\u003c/em\u003e). Russ J Bioorg Chem 43:397\u0026ndash;408\u003c/li\u003e\n\u003cli\u003eWang QL, Zhang JY, Liu ZS, Duan YY, Li CY (2024) Integrative approaches based on genomic techniques in the functional studies on enhancers. Brief Bioinform 25:bbad442\u003c/li\u003e\n\u003cli\u003eXiong X, Liu JK, Rao YL (2023) Whole genome resequencing helps study important traits in chickens. Genes 14:1198\u003c/li\u003e\n\u003cli\u003eLi L, He S, Lin MH, Zhang YP, Kuhl H, Liang XF (2023) Whole-genome resequencing and bisulfite sequencing provide new insights into the feeding habit domestication in mandarin fish (\u003cem\u003eSiniperca chuatsi\u003c/em\u003e). Front Genet 13:1088081\u003c/li\u003e\n\u003cli\u003eSun CF, Zhang XH, Dong JJ, You XX, Tian YY, Gao FY et al. (2023) Whole-genome resequencing reveals recent signatures of selection in five populations of largemouth bass (\u003cem\u003eMicropterus salmoides\u003c/em\u003e). Zool Res 44:78\u0026ndash;89\u003c/li\u003e\n\u003cli\u003eHuang WW, Lai J, Liang WQ, Ye SZ, Li JW, Zhou JW et al. (2024) Identification of sex-specific DNA markers in the army fish (\u003cem\u003eSpinibarbus hollandi\u003c/em\u003e) by whole genome re-sequencing method. Aquaculture 583:740605\u003c/li\u003e\n\u003cli\u003eWang J, Ma LX, Dong YW (2023) Coping with harsh heat environments: molecular adaptation of metabolic depression in the intertidal snail \u003cem\u003eEchinolittorina radiata\u003c/em\u003e. Cell Stress Chaperon 28:477\u0026ndash;491\u003c/li\u003e\n\u003cli\u003eStillman JH, Amri AB, Holdreith JM, Hooper A, Leon RV, Pruett LR et al. (2025) Ecophysiological responses to heat waves in the marine intertidal zone. J Exp Biol 228:JEB246503\u003c/li\u003e\n\u003cli\u003eYou X, Bian C, Zan Q (2014) Mudskipper genomes provide insights into the terrestrial adaptation of amphibious fishes. Nat Commun 5:5594\u003c/li\u003e\n\u003cli\u003eLin XN, Ma CY, Hu LS, Liao ML, Ma LX, Teske PR et al. (2024) Genomics-informed range predictions under global warming reveal reduced adaptive diversity whilst buffering range shifts for a marine snail. Glob Change Biol 30:e17571\u003c/li\u003e\n\u003cli\u003eWang CT, Liu ZZ, Ma LB (2012) Isolation and characterization of 16 polymorphic microsatellite markers for the mudskipper (\u003cem\u003eBoleophthalmus pectinirostris\u003c/em\u003e). Conserv Genet Resour 4:1059\u0026ndash;1061\u003c/li\u003e\n\u003cli\u003eChen W, Hong W, Chen S (2015) Population genetic structure and demographic history of the mudskipper \u003cem\u003eBoleophthalmus pectinirostris\u003c/em\u003e on the northwestern Pacific coast. Environ Biol Fish 98:845\u0026ndash;856\u003c/li\u003e\n\u003cli\u003eChen S, Zhou Y, Chen Y, Gu J (2018) fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34:i884\u0026ndash;i890\u003c/li\u003e\n\u003cli\u003eBian C, Huang Y, Li RH, Xu PW, You XX, Lv YY et al. (2024) Genomics comparisons of three chromosome-level mudskipper genome assemblies reveal molecular clues for water-to-land evolution and adaptation. J Adv Res 58:93\u0026ndash;104\u003c/li\u003e\n\u003cli\u003eLi H, Durbin R (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25:1754\u0026ndash;1760\u003c/li\u003e\n\u003cli\u003eLi H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N et al. (2009) The Sequence Alignment/Map format and SAMtools. Bioinformatics 25:2078\u0026ndash;2079\u003c/li\u003e\n\u003cli\u003eMcKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A et al. (2010) The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res 20:1297\u0026ndash;1303\u003c/li\u003e\n\u003cli\u003eWang K, Li M, Hakonarson H (2010) ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res 38:e164\u003c/li\u003e\n\u003cli\u003eJones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C et al. (2014) InterProScan 5: genome-scale protein function classification. Bioinformatics 30:1236\u0026ndash;1240\u003c/li\u003e\n\u003cli\u003eFelsenstein J (1993) PHYLIP (phylogeny inference package), version 3.5 c. University of Washington, Seattle.\u003c/li\u003e\n\u003cli\u003ePurcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D et al. (2007) PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet 81:559\u0026ndash;575\u003c/li\u003e\n\u003cli\u003eAlexander DH, Novembre J, Lange K (2009) Fast model-based estimation of ancestry in unrelated individuals. Genome Res 19:1655\u0026ndash;1664\u003c/li\u003e\n\u003cli\u003eZhang C, Dong SS, Xu JY, He WM, Yang TL (2019) PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics 35:1786\u0026ndash;1788\u003c/li\u003e\n\u003cli\u003eDanecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA et al. (2011) The variant call format and VCFtools. Bioinformatics 27:2156\u0026ndash;2158\u003c/li\u003e\n\u003cli\u003eLi H, Durbin R (2011) Inference of human population history from individual whole-genome sequences. Nature 475:493\u0026ndash;496\u003c/li\u003e\n\u003cli\u003ePickrell JK, Pritchard JK (2012) Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet 8:e1002967\u003c/li\u003e\n\u003cli\u003eSzpiech ZA, Hernandez RD (2014) selscan: An Efficient Multithreaded Program to Perform EHH-Based Scans for Positive Selection. Mol Biol Evol 31:2824\u0026ndash;2827\u003c/li\u003e\n\u003cli\u003eCantalapiedra CP, Hern\u0026aacute;ndez-Plaza A, Letunic I, Bork P, Huerta-Cepas J (2021) eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale. Mol Biol Evol 38:5825\u0026ndash;5829\u003c/li\u003e\n\u003cli\u003eAshburner M, Ball C, Blake J (2000) Gene Ontology: tool for the unification of biology. Nat Genet 25:25\u0026ndash;29\u003c/li\u003e\n\u003cli\u003eKanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M (2023) KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res 51:D587\u0026ndash;D592\u003c/li\u003e\n\u003cli\u003eLangley CH, Stevens K, Cardeno C, Lee YCG, Schrider DR, Pool JE et al. (2012) Genomic variation in natural populations of Drosophila melanogaster. Genetics 192:533\u0026ndash;598\u003c/li\u003e\n\u003cli\u003eZhou J, He W, Wang J (2023) The pan-plastome of tartary buckwheat (\u003cem\u003eFagopyrum tataricum\u003c/em\u003e): key insights into genetic diversity and the history of lineage divergence. BMC Plant Biol 23:212\u003c/li\u003e\n\u003cli\u003eMa DN, Lai ZF, Ding QS, Zhang K, Chang KZ, Li SH et al. (2022) Identification, characterization and function of orphan genes among the current Cucurbitaceae genomes. Front Plant Sci 13:872137\u003c/li\u003e\n\u003cli\u003eChen H, Polgar G, Yin W, Fu CZ (2014) Cryptic species and evolutionary history of \u003cem\u003eBoleophthalmus pectinirostris\u003c/em\u003e complex along the northwestern Pacific coast. Acta Hydrobiol Sin 38:75\u0026ndash;86\u003c/li\u003e\n\u003cli\u003eYang F, He LJ, Lei GC, Zhang AB (2012) Genetic diversity and DNA barcoding of mudskipper common species along southeast coasts of China. Chin J Ecol 31:676\u0026ndash;683\u003c/li\u003e\n\u003cli\u003eMcKeown NJ, Arkhipkin AI, Shaw PW (2019) Genetic analysis reveals historical and contemporary population dynamics in the longfin squid Doryteuthis gahi: implications for cephalopod management and conservation. ICES J Mar Sci 76:1019-1027\u003c/li\u003e\n\u003cli\u003eWang D, Gao L, Tian H, et al. (2020) Population genetics and sympatric divergence of the freshwater gudgeon, Gobiobotia filifer, in the Yangtze River inferred from mitochondrial DNA. Ecol Evol 10:50\u0026ndash;58\u003c/li\u003e\n\u003cli\u003eWu B, Chen X, Yu M, Ren J, Hu J, Shao C et al. (2022) Chromosome-level genome and population genomic analysis provide insights into the evolution and environmental adaptation of Jinjiang oyster \u003cem\u003eCrassostrea ariakensis\u003c/em\u003e. Mol Ecol Resour 22:1529\u0026ndash;1544\u003c/li\u003e\n\u003cli\u003eBernos TA, Avlija\u0026scaron; S, Hill J, Morissette O, Ricciardi A, Mandrak NE et al. (2022) Genetic diversity and structure of a recent fish invasion: Tench (\u003cem\u003eTinca tinca\u003c/em\u003e) in eastern North America. Evol Appl 16:173\u0026ndash;188\u003c/li\u003e\n\u003cli\u003eGao K, He Z, Xiong J (2024) Population structure and adaptability analysis of \u003cem\u003eSchizothorax o\u0026apos;connori\u003c/em\u003e based on whole-genome resequencing. BMC Genomics 25:145\u003c/li\u003e\n\u003cli\u003eKon T, Pei L, Ichikawa R (2021) Whole-genome resequencing of large yellow croaker (\u003cem\u003eLarimichthys crocea\u003c/em\u003e) reveals the population structure and signatures of environmental adaptation. Sci Rep 11:11235\u003c/li\u003e\n\u003cli\u003eZhang G, Fang X, Guo X (2012) The oyster genome reveals stress adaptation and complexity of shell formation. Nature 490:49\u0026ndash;54\u003c/li\u003e\n\u003cli\u003eGoh KM, Gonz\u0026aacute;lez-Siso MI, Sani RK (2023) Genomics of extreme environments: unveiling the secrets of survival. Sci Rep 13:21441\u003c/li\u003e\n\u003cli\u003eZhang X, Yan X (2014) Spatiotemporal change in geographical distribution of global climate types in the context of climate warming. Clim Dyn 43:595\u0026ndash;605\u003c/li\u003e\n\u003cli\u003eYang D, Li W, Xiang P, Ge T, Li H, Zhang Y (2024) Rhein promotes skin wound healing by activating the PI3K/AKT signaling pathway. Open Med 19:20241116\u003c/li\u003e\n\u003cli\u003eMartini M, De Santis MC, Braccini L, Gulluni F, Hirsch E (2014) PI3K/AKT signaling pathway and cancer: an updated review. Ann Med 46:372\u0026ndash;383\u003c/li\u003e\n\u003cli\u003eShin JM, Kim MY, Sohn KC, Jung SY, Lee HE, Lim JW et al. (2014) Nrf2 negatively regulates melanogenesis by modulating PI3K/Akt signaling. PLoS One 9:e96035\u003c/li\u003e\n\u003cli\u003eTang H, Liu J, Wang Z, Zhang L, Yang M, Huang J et al. (2023) Genome-wide association study (GWAS) analysis of black color trait in the leopard coral grouper (\u003cem\u003ePlectropomus leopardus\u003c/em\u003e) using whole genome resequencing. Comp Biochem Physiol Part D Genomics Proteomics 48:101138\u003c/li\u003e\n\u003cli\u003eCoelho SG, Hearing VJ (2010) UVA tanning is involved in the increased incidence of skin cancers in fair-skinned young women. Pigment Cell Melanoma Res 23:57\u0026ndash;63\u003c/li\u003e\n\u003cli\u003eQian W, Deng L, Lu D, Xu S (2013) Genome-wide landscapes of human local adaptation in Asia. PLoS One 8:e54224\u003c/li\u003e\n\u003cli\u003eNie MM, Ni WL, Wang LH, Gao Q, Liu D, Tian F et al. (2022) Insights into miRNA-mRNA regulatory mechanisms of cold adaptation in \u003cem\u003eGymnocypris eckloni\u003c/em\u003e: ubiquitin-mediated proteolysis is pivotal for adaptive energy metabolism. Front Genet 13:903995\u003c/li\u003e\n\u003cli\u003eZhou X, Cao J, Xie J, Tong W, Jia B, Fu J (2025) Effect of Dync1h1 on phototransduction protein transport and the development and maintenance of photoreceptor cells in zebrafish. Invest Ophthalmol Vis Sci 66:38\u003c/li\u003e\n\u003cli\u003eHaase J, Jones AKC, McVeigh CJ, Brown E, Clarke G, Ahnert-Hilger G (2021) Sex and brain region-specific regulation of serotonin transporter activity in synaptosomes in guanine nucleotide-binding protein G(q) alpha knockout mice. J Neurochem 159:156\u0026ndash;171\u003c/li\u003e\n\u003cli\u003eWang C, Yang L, Xiao TY, Li JH, Liu QL, Xiong ST (2022) Identification and expression analysis of zebrafish gnaq in the hypothalamic-pituitary-gonadal axis. Front Genet 13:1015796\u003c/li\u003e\n\u003cli\u003ePridgeon JW, Klesius PH (2013) G-protein coupled receptor 18 (GPR18) in channel catfish: expression analysis and efficacy as immunostimulant against \u003cem\u003eAeromonas hydrophila\u003c/em\u003e infection. Fish Shellfish Immunol 35:1070\u0026ndash;1078\u003c/li\u003e\n\u003cli\u003eKiryu I, K\u0026ouml;llner B, Kuroda A, Ototake M, Dijkstra JM (2003) A new putative G-protein coupled receptor gene associated with the immune system of rainbow trout (\u003cem\u003eOncorhynchus mykiss\u003c/em\u003e). Fish Shellfish Immunol 15:117\u0026ndash;127\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Boleophthalmus pectinirostris, population genetic divergence, whole-genome resequencing, adaptive evolution, environmental adaptation","lastPublishedDoi":"10.21203/rs.3.rs-6768090/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-6768090/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"Boleophthalmus pectinirostris is an amphibious fish widely distributed in the coastal intertidal zone, which is a good object for evolutionary biology research. With the increasing application of whole-genome resequencing technology, genetic evolution analysis and gene screening can better clarify the genetic basis and molecular mechanisms of environmental adaptation. In the present study, B. pectinirostris samples were collected across five coastal regions of China, and their population genetic viarations, demographic history alteration, and environmental adaptation mechanisms were analyzed. After alignment with the reference genome, a total of 33,016,712 high-quality SNPs were obtained for 68 specimens. The results showed that population Qingdao (QD) exhibited strong genetic heterogeneity with other four populations, and its low genetic diversity may be related to its long-term population decline. The significant changes in the effective population size of different populations over history were revealed, with notable differences between northern and southern groups. In addition, Selective sweep analysis uncovered northern-adapted genes in B. pectinirostris, linked to energy metabolism and stress response under colder environments. The results of this study contribute to a better understanding of the adaptive evolution of intertidal fishes, particularly the mechanisms associated with temperature, seasonal fluctuations and other adaptations. These findings will provide an important scientific basis for the conservation and sustainable management of this species.","manuscriptTitle":"Evolutionary Clues of North-South Divergence: Genomic Adaptation of the Intertidal Mudskipper Boleophthalmus pectinirostris","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-06-16 10:00:31","doi":"10.21203/rs.3.rs-6768090/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"63cbe67e-c79d-4cb8-a89d-9fbedbfc5469","owner":[],"postedDate":"June 16th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":49909794,"name":"Biological sciences/Evolution/Population genetics/Genetic variation"},{"id":49909795,"name":"Biological sciences/Genetics/Genome/Genetic variation"}],"tags":[],"updatedAt":"2025-09-23T11:20:38+00:00","versionOfRecord":[],"versionCreatedAt":"2025-06-16 10:00:31","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-6768090","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-6768090","identity":"rs-6768090","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}