Seven Uranoscopidae species along the coast of China were collected ( Table 1 , Figure 1 ). All samples were cryopreserved and transported to the Marine Fishery Resource and Biodiversity Laboratory of Zhejiang Ocean University. Preliminary morphological identification was conducted following the methods outlined by Chen et al (Kawai et al., 2003; Chen and Zhang, 2016; Wu and Zhong, 2021). Approximately 1 g of muscle tissue was extracted from each sample for DNA extraction. All animal specimens were collected in compliance with legal regulations and adhered to the Animal Care and Use Ethics policies of Zhejiang Ocean University. Every effort was made to minimize invasiveness during the collection process.

Phenol/Chloroform extraction method was used to extract the DNA from muscle tissue. The extracted DNA concentration was assessed using NanoDrop 2000 (Thermo Fisher Scientific Inc, USA) and verified through 1% agarose gel electrophoresis. After ultrasonic fragmentation, libraries were constructed with insert fragment sizes of approximately 350 bp for sequencing. All libraries were sequenced on the Illumina NovaSeq 6000 platform according to the manufacturer′s protocol.

FASTP (Chen et al., 2018) with default parameters was utilized for raw data filtering and quality control, resulting in the acquisition of clean data. Following the step, a total of 396.89 Gb (averaging 100× coverage of the genome) of high-quality data were retained for further analysis. To assess sequencing quality, metrics such as quality values Q20 and Q30, along with GC content, were calculated. GCE v1.0.0 (Liu et al., 2013) was used to estimate the genome characteristics of each sample, with K-mer size of 17. The outcomes of the K-mer analysis were leveraged to estimate genome size, heterozygosity, and repeat ratio. The algorithm used for determining genome size is as follows: genome size = K-mer_num/peak_depth, where K-mer_num represents the total number of K-mers, and peak_depth is the expected value of the K-mer depth.

SOAPdenovo2 (Luo et al., 2012) was employed to individually assemble the clean reads into distinct contigs and scaffolds. The assembly genome sequences were evaluated by BUSCO using actinopterygii_odb10 database (Seppey et al., 2019). Simultaneously, the conserved gene sequences of each genome were identified. For SNP calling, the clean reads were mapped to X. elongatus genome using BWA-mem method (Li, 2014) and the variance were identified by GATK (v 4.1.9.0) (DePristo et al., 2011) with parameters of “QD < 2.0, MQ < 40.0, FS > 60.0, SOR > 3.0, MQRankSum < -12.5, ReadPosRankSum < -8.0, DP < 8.0”. The raw variances were filtered using vcf2phylip (v2.0) (Ortiz, 2019). SNPs with integrity of 100% were analyzed in the followed study. The identification of potential microsatellite motifs was conducted using the Perl script “misa.pl” from the MISA software (Beier et al., 2017). The search parameters were configured to detect dinucleotide, trinucleotide, tetranucleotide, pentanucleotide, and hexanucleotide microsatellite motifs with a range of five-to-fifteen repeats.

The clean data were employed for assembling the mitochondrial genome using NOVOplasty v4.2 (Dierckxsens et al., 2017). The resulting mitochondrial genome sequences were annotated and visualized as circular maps using OrganellarGenomeDRAW (Greiner et al., 2019). As an outgroup, mitochondrial genome sequences of Gymnammodytes semisquamatus (16,484 bp, MT410950), Ammodytes personatus (16,527 bp, KF672362), Parapercis cylindrica (16,195 bp, AP006022), Trachinus draco (16,810 bp, AP006024), and Parapercis sexfasciata (16,497 bp, AP006816) were retrieved from NCBI. RAxML v8.2.12 (Stamatakis, 2014) was employed to reconstruct the maximum-likelihood (ML) tree using the optimal GTRGAMMAI substitution model and 1,000 bootstrap replicates. For constructing the mitochondrial genome based phylogenetic tree, all genes were stringed together to form a supergene as input. For constructed the core-genome origined conserved gene tree, all BUSCO identified genes were clustered using OrthMCL to find the single copy gene families. For each gene families, a tree was constructed and the final tree was generated by merging tree of each gene family using ASTRAL (Sayyari and Mirarab, 2016).

PSMC (Li and Durbin, 2011) was employed for inferring the effective population size of each individual. Specifically, clean reads were aligned to the assembled genome sequence using BWA-mem method. Samtools (v0.1.19) (Danecek et al., 2021) was used to handing the mapped bam file using parameter of “-bF 12”. Then, bcftools and vcftools were used to convert the sorted bam files into “fq.gz” file, and “fq2psmcfa” script in PSMC was used to convert the “fq.gz” file into psmcfa file with parameter of “-q20”. For running PSMC, the generation interval (g) was set to 2 years, and the mutation rate (μ) was set at 2.5×10^-8 (Song et al., 2022), with the remaining parameters set to their default values.

After conducting the data quality control step, a total of 396.69 Gb of high-quality sequencing data was successfully obtained from seven samples, with each sample’s data ranging from 41.77 to 74.46 Gb. All samples exhibited Q30 values surpassing 91.96%, indicating an exceptionally high overall sequencing quality. The sequencing depth for each sample ranged from 64× to 122×, averaging at 87×, meeting the prerequisites for subsequent genome prediction ( Supplementary Table S1 ).

K-mer analysis was employed to evaluate genome size, heterozygosity rates, and the proportion of repetitive sequences for different samples ( Table 2 ). The assessment revealed that all samples were diploid, with genome sizes ranging from 536 Mb to 650 Mb. Notably, U. japonicus displayed the largest genome at 652.49 Mb, while U. bicinctus exhibited the smallest genome at 536 Mb. In comparison to cyprinid fish, the Uranoscopidae fishes showcased relatively low proportions of repetitive sequences, ranging from 20.09% to 36.64%. Despite all seven fish species belonging to the Uranoscopidae family, there were considerable differences in heterozygosity rates among individuals, none exceeding 1%. U. oligolepis recorded the highest heterozygosity rate at 0.88%, while I. pollicaris exhibited the lowest at 0.41%. This implies significant diversity in effective population size among these fish species under natural conditions.

We employed the SOAPdenovo2 software for assembling clean data, yielding genome sizes ranging from 533 Mb to 669 Mb ( Table 3 ). The assembled genome sequences for all samples accounted for 85.54% to 114.09% of the genome sizes assessed by K-mer method, indicating relative completeness in all assembled genome sequences. Due to the use of next-generation sequencing (NGS) data for assembly, the scaffold N50 for all seven sequences was relatively low, ranging from 7.3 Kb to 62.5 Kb. Besides, the genome continuity was suboptimal, with varying numbers of scaffolds between 172,912 and 837,792. Evaluating genome completeness using BUSCO revealed significant differences in completeness among different assembly sequences, ranging from 14.6% to 86.3% ( Supplementary Table S2 ). Among the seven samples, U. sp. exhibited the highest number of contigs with the lowest contig N50, resulting in low complete BUSCOs. Conversely, X. elongatus demonstrated the best assembly performance, boasting the highest scaffold N50 at 62,497 bp and the longest contig reached 1.14 Mb. Meanwhile, BUSCO evaluation indicated an impressive completeness of 86.30%, with 85.5% of BUSCO genes being complete and single-copy. This is reasonable, as X. elongatus genome exhibited very low heterozygosity (0.41%) and relatively low repetitive sequence proportion (20.09%).

Utilizing the assembled genome, we performed the prediction of simple sequence repeats (SSR). The predicted SSR counts for the seven species closely aligned, ranging from 248,600 to 294,869 ( Table 4 ). When assessed in relation to the assembled sequence lengths, the SSR density exhibited striking similarity across different species, with each Mb sequence containing 403-510 SSRs. The SSR composition patterns are highly similar across all seven species ( Figure 2A ). In general, SSRs with higher unit repetition numbers constitute a smaller proportion of the overall composition. This trend is observed in Trinucleotide, Tetranucleotide, Pentanucleotide, and Hexanucleotide SSRs. As for Dinucleotide SSRs, with the exception of X. elongatus, both the six-unit and ten-unit repetitions are notably higher compared to SSRs with other repetition numbers in the remaining six species.

To explore the evolutionary relationships among different species, we used zebrafish as an outgroup and constructed an evolutionary tree based on conserved genes in the genomes. The conserved genes identified by the BUSCO prediction model for each species were utilized as input, resulting in the identification of a total of 90 shared single-copy genes. Subsequently, we employed the maximum likelihood method to construct an evolutionary tree comprising 8 species ( Figure 3 ). The results reveal that I. pollicaris and X. elongatus both belong to the Ammodytidae family, clustering together on the evolutionary tree. The genus Uranaosopus belongs to the Uranoscopodae family, with its 5 species forming a cluster on a branch. The bootstrap support rate between the two families is 100%. Within the genus Uranaosopus, U. japonicus is evolutionarily closer to I. pollicaris and X. elongatus, with a bootstrap support rate of 100%. The support rate among other species within the genus Uranaosopus also exceeds 76%.

To further validate the evolutionary relationships within the nuclear genome of the seven species, we utilized whole-genome genetic variance to construct an additional phylogenetic tree. The X. elongatus genome was chosen as the reference genome due to its superior integrity and continuity, evidenced by the highest BUSCO completeness value, contig N50, and scaffold N50, as well as the lowest scaffold number among the seven assembled genomes. A total of 3,756,984 SNPs, constituting 0.69% of the assembled X. elongatus contig sequence, were identified across the seven species, with the number of single species SNP variants amount ranging from 2.19 million to 2.49 million ( Supplementary Table S3 ). After filtering variant loci with integrity below 100%, 1,771,650 SNPs were retained for construction the phylogenetic tree. Notably, the resulting tree demonstrated a high level of consistency with the tree derived from common conserved genes, thus reaffirming the reliability of the nuclear genome phylogenetic relationships ( Figure 4 ).

The assembled total lengths of the mitogenomes for U. tosae, U. japonicus, U. sp., U. bicinctus, U. oligolepis, I. pollicaris, and X. elongatus were 17,532 bp, 18,446 bp, 15,763 bp, 14,966 bp, 16,675 bp, 16,456 bp, and 16,546 bp, respectively ( Figure 2B ). Each mitochondrial genome encompasses a variable number of genes, ranging from 36 to 39 ( Table 5 ). Typically, except for one NADH dehydrogenase coding gene (ND6) and eight tRNA genes (Glu, Ala, Asn, Cys, Tyr, Ser, Gln, and Pro) were encoded on the L-strand, most mitochondrial genes were encoded on the H-strand. All 13 protein-coding genes, comprising 7 NADH dehydrogenases, 3 cytochrome c oxidases, 1 cytochrome b, and 2 ATP synthases, were identified in these mitochondrial sequences. The number of transfer RNA (tRNA) and ribosomal RNA (rRNA) genes varies among the seven sequences, with most species having 23 tRNA genes and 1 rRNA gene.

To further unravel the evolutionary relationships among different species within the genus Uranaosopus, we constructed a species phylogenetic tree using all conserved genes in the mitochondrial genome. As outgroups, we selected five species from the Perciformes order, including two from the Ammodytidae family (Gymnammodytes semisquamatus and Ammodytes personatus), two from the Pinguipedidae family (Parapercis cylindrica and Parapercis sexfasciata), and one from the Trachinidae family (Trachinus draco). At the suborder level, all species, except T. draco, belong to Uranoscopoidei, with T. draco serving as the outgroup for the entire tree.

The mitochondrial genome-based gene tree exhibits strong consistency with the nuclear-genome-conserved-genes based tree. I. pollicaris and X. elongatus, along with U. japonicus and U. bicinctus, maintain consistent positions of two trees. However, there is inconsistency in the relative relationships among the other three species in genus Uranoscopus. Although the bootstrap support rates among these three species exceed 50%, they are slightly lower compared to other consistent branches ( Figure 5 ).

We conducted a population history analysis of various species within the Uranoscopidae family using PSMC ( Figure 6 ). The results reveal that following the Last Interglacial Period (LIP), all species underwent a rapid expansion in population size. By the time of the Last Glacial Period (LGP), the effective populations size of U. bicinctus and U. japonicus started to decline, experiencing a bottleneck effect. U. oligolepis began to undergo a bottleneck effect towards the end of the Last Glacial Period (LGP). Notably, based on the speculative results, U. tosae, X. elongatus, U. sp., and I. pollicaris have remained unaffected by bottleneck events since the LIP. Furthermore, significant differences in effective population sizes are observed among these species. Although U. sp. and I. pollicaris underwent expansion in one branch after LIP, they still maintained relatively small populations. In comparison, U. bicinctus, U. oligolepis, and U. japonicus consistently maintained higher population levels.