A total of 180 samples from the L. vannamei broodstock of Renhai No. 1 (RH), Kehai No. 1 (KH), Benchmark Genetics (BMK), Charoen Pokphand (CP), Shrimp Improvement Systems (SIS), and Top Aquaculture Technology (TA) were collected from Hairen Aquatic Seed Industry Technology Co., Ltd (Hebei, China) (Supplementary Table 1). The muscle samples of each individual were collected for DNA extraction and genome resequencing. Genomic DNA was extracted using the TIANamp Marine Animal DNA Kits (TIANGEN Biotech Co., Ltd. Beijing, China). Paired-end sequencing libraries, with an insert size of 250 to 350 bp, were constructed by the VAHTS Universal Plus DNA Library Prep Kit for MGI (Vazyme Biotech Co., Ltd., Nanjing, China) in their lab, and sequenced on an MGI DNBSEQ-T7 system (BGI Genomics Co., Ltd., Shenzhen, China). All raw sequencing data were subject to quality control¹. Low-quality bases and reads were trimmed by Trimmomatic with parameters “LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36” (Bolger et al., 2014).

The reference genome was downloaded from the National Center for Biotechnology Information (NCBI) with the accession of GCA_003789085.1 (Zhang X. et al., 2019). We anchored the contigs into scaffolds by the guidance of a genetic linkage maps (Yu et al., 2015; Jones et al., 2017), and all trimmed reads were aligned to this genomic assembly using the BWA (version 2.3.4.1), with the default parameters (Li and Durbin, 2009). The Binary Alignment Map (BAM) files were imported to the samtools (v0.1.19) for reference index building, format conversion, and reads sorting (Li et al., 2009). The Picard² (v1.92) and the sambamba (v0.8.0) were used to assign read group information (Tarasov et al., 2015). PCR and optical duplicates were marked and filtered using GATK (v4.1). Variants and haplotypes were identified using the HaplotypeCaller algorithm in Genomic Varient Call Format (GVCF) mode (Mckenna et al., 2010). The GenotypeGVCFs were used subsequently for joint genotyping.

To keep the most reliable SNPs for subsequent analysis, variant sites were marked and filtered by the vcffilter with the following criteria: low quality score (GQ < 20); low quality by depth score (QD < 10); high Fisher strand score (FS > 10); low mapping quality (MQ < 40); low read position rank sum score (ReadPosRankSum < 8); high strand odds ratio (SOR > 4); and low mapping quality rank-sum score (MQRankSum < 12.5). The SNP loci with multiple alleles, low minor allele frequency (MAF < 0.05), and missing genotypes (max-missing < 1) were also removed using the VCFtools (v0.1.16) (Danecek et al., 2011). The potential effect of each SNP was annotated using the SNPEff (Cingolani et al., 2012).

Principal component analysis (PCA) was performed using the PLINK (Purcell et al., 2007). Ancestry estimation and population structure were analyzed using the ADMIXTURE (Alexander et al., 2009). Pairwise fixation index (F_st) was calculated by VCFtools with parameters –fst-window-size 50,000 and –fst-window-step 10,000 (Danecek et al., 2011). The genomic observed heterozygosity, expected heterozygosity, and inbreeding coefficient were evaluated by the VCFtools with parameter –het (Danecek et al., 2011). For the phylogenetic analysis, representative SNP markers were extracted using the VCFtools with parameter –thin 10,000 to mitigate the effects of LD, homozygous sites were removed from the alignment file. Maximum likelihood phylogenetic inference was carried out by IQ-TREE (Nguyen et al., 2015).

The SNPs with specific alleles in one population were defined as population-specific SNPs (ps-SNPs). If the frequency of the population-specific allele was greater than 50% in a population, the SNPs were defined as common population-specific SNPs (cps-SNPs). The LD decay was estimated for each population using the PopLDdecay (Zhang C. et al., 2019). Genomic inbreeding coefficients F_ROH were calculated according to the following equation, with a minimum window of 100 kb (McQuillan et al., 2008):

The pairwise F_st and the ratio of genetic polymorphisms (θ_π) were calculated using a 50 kb window and a 10 kb step between each breed. Genomic regions with signals of differences in the spectrum of genetic polymorphisms were identified as regions that are potentially under selection. Empirical cut-offs for F_st and θ_π ratio were set as the top 2% largest and the top 5% largest or smallest, respectively (Li L. et al., 2018; Wang et al., 2021). Candidate genes under selection were defined as within or overlapping with the regions showing signals of selection. Candidate genes were characterized with Gene Ontology (GO) enrichment analysis implemented in the EnrichPipeline based on Fisher’s exact test (Huang et al., 2009). The GO terms with a P-value < 0.05 were considered significantly enriched.

An average of 164 million 300 bp paired-end reads were generated for each library, representing ∼18.7 × coverage. After quality control, approximately 80% of reads could be aligned to the reference genome (∼12.34 × mean unique coverage) and used for variant calling. A total of 37,308,098 SNPs were identified in the six populations. The number of SNPs in the chromosomes ranged from 182,599 to 1,274,977, and the density was between 20.60 SNPs/Kb and 24.58 SNPs/kb. Of these, 4,509,526 reliable SNPs were kept after the removal of low-quality samples and those within repetitive regions of missing genotypes. The transitions/transversions ratio in this biallelic SNPs set was found to be 1.63, with 1,117,185,030 transitions and 72,199,828 transversions, respectively. The SNPs were unevenly distributed across the genome, with over 10.56% (1,460,366) located in the intergenic regions, about 55.43% (7,662,564) located in the introns, and only 2.283% (238,482) located in the exons (Table 1). More than 73.89% of the SNPs (221,021) were silent mutations, the remaining 25.82 (77,219) and.30% of SNPs (886) could result in missense and non-sense mutations, respectively (Table 1).

The population structure was inferred based on the genome-wide biallelic SNPs. As shown by the PCA, the first three components cumulatively explained 24.56% of the variance (Figure 1A and Supplementary Figures 1, 2). Samples from the six populations were grouped into three major clusters, which largely restored the sampling sources. Obvious stratifications were observed among BMK, SIS, and the Asian breeds (RH, KH, TA, and CP). Although RH, KH, TA, and CP were distributed closely on the first two components, they possessed distinct population structures (Figure 1 and Supplementary Figure 2). The pairwise comparison revealed that the F_st values among the six populations ranged from 0.0005 to 0.13415, and that BMK had the largest genetic distances from the others. The KH and TA share a relatively higher similarity in genetic structure (Figure 1B). Populations defined by estimated genetic ancestries also confirmed the distribution of samples, where six ancestral populations were inferred to be the optimal scenario by cross-validation (Figure 1C and Supplementary Figure 3). The phylogenetic relationships were reconstructed via the maximum likelihood inference. Except for two samples of TA clustered together with CP, samples from each population formed monophyletic clades, which is consistent with the population structure inferred by PCA and ancestral estimation (Figure 1D).

As shown in Figure 2A, population-specific SNPs were identified in each breed. The BMK and SIS possessed 83,273 and 15,980 cps-SNPs, respectively, which are much larger than that of the other four breeds. The RH, KH, TA, and CP only had 3,153, 1,889, 485, and 836 cps-SNPs, respectively (Table 2). We evaluated the genetic diversity by the observed heterozygosity, expected heterozygosity, and inbreeding coefficient (Supplementary Data Sheet 1). These six populations have an average genomic observed heterozygosity of 0.122, 0.137, 0.144, 0.147, 0.149; and 0.156 in BMK, CP, TA, RH, KH, and SIS. The extent of genome-wide LD decay was measured against the physical distance in each population (Figure 2B). Asian breeds (RH, KH, TA, and CP) typically exhibited an overall lower level of LD than the others. Their LD levels (r²) fell below 0.1 at less than 15Kb, whereas, the physical distances reached over 70 kb in BMK and SIS. Analysis of runs of homozygosity (ROH) revealed that BMK and CP reserved higher levels of inbreeding than the others (Figure 2C and Supplementary Figure 4). Homozygous segments larger than 1Mb were mostly identified in these two populations. Despite that the genome-wide distribution of ROH is distinct in each population, peaks in specific regions of chromosomes 10, 15, and 42 were identified in each population (Figure 2D).

The pairwise F_st and θ_π ratio were calculated to scan the genomic regions with genetic signals of divergence. The top 5% of the windows, with high values of F_st and differentiation of polymorphism frequency spectrum, were identified as regions that are potentially under selection (Figure 3A and Supplementary Figure 5). Genome-wide selective sweeps in BMK were illustrated in Figure 3A. We detected 271, 440, 257, 269, 431, and 429 candidate genes under a strong selection from 166,258 genomic windows in RH, BMK, KH, TA, SIS, and CP, respectively (Supplementary Data Sheet 2). Interestingly, 206 genes were under selection in more than one population (Supplementary Data Sheet 2). Function enrichment analysis of these candidate genes provided 631 significant enriched categories (Supplementary Data Sheet 3). For example, genes under selection in BMK were enriched in response to stimulus, carbohydrate transport, and several biological processes (Figure 3B). Selective genes in RH were mainly enriched in trialkyl sulfonium hydrolase activity, protein modification by a small protein conjugation, and phospholipid binding. Selective genes in KH were mainly enriched in the structural constituent of ribosome, structure molecule activity, and cellular macromolecule metabolic process. Selective genes in TA were enriched in the pyruvate kinase activity, peptidase inhibitor activity, and ATP metabolic process. Selective genes in SIS were enriched in the ribonuclease H2 complex, GTP cyclohydrolase I activity, and cellular macromolecule metabolic process. Selective genes in CP were mainly enriched in the structural molecule activity, regulation of ATPase activity, and 5′-deoxynucleotidase activity (please check Supplementary Data Sheet 3 for detailed in section “Results”). Notably, B3GT5, CDC42, PXDN, and several genes that play important roles in cellular macromolecule metabolic process, proteolysis, and responses to stimulus were identified under selection in more than one population (Figure 4).