In silico bioinformatic analyses for 2B-type enzyme definition

The in silico identification of four 2B-type restriction enzymes (BcgI, BsaXI, AlfI, and CspCI; mat.compl1, mat.compl2) was done on the genomes of three aquaculture important species: whiteleg shrimp Penaeus (Litopenaeus) vannamei (GenBank accession no.: ASM378908v1; ^{Zhang et al., 2019}), Pacific oyster Crassostrea gigas (GenBank accession no.: cgigas_uk_roslin_v1 NC_047559-NC_047568; ^{Peñaloza et al., 2021}), and almaco jack Seriola rivoliana (GenBank accession no.: GCA_002994505.1 ASM299450v1; ^{Seetharam, 2018}). The genomes of the species were obtained from the GenBank database (National Center for Biotechnology Information https://www.ncbi.nlm.nih.gov/). They were downloaded to the local server at the Genomics and Bioinformatics Laboratory of Centro de Investigaciones Biológicas del Noroeste, S.C. (CIBNOR).

The number of recognition sites of each enzyme on the genomic sequencing on sense- and anti-sense directions were counted with search and line-count commands (mat.compl3). The next step consisted of estimating the remaining sites depending on the number of selective bases included in the adaptors.

The proportion of variable sites relative to the number of recognition sites in the whiteleg shrimp was estimated. A previously generated database of genomic analysis of 88 individuals of several origins was used (^{Perez-Enriquez et al., 2018}). This database was depurated only to contain fragments with the recognition sites of the enzyme BcgI (zgrep command) and aligned to the complete genome of ^{Zhang et al. (2019}) following the procedure in Perez-Enriquez et al. (2018).

2bRAD technique standardization

The standardization of the 2bRAD technique was done on a reduced number of samples of the species of interest to verify the number of markers to be obtained with the selected enzyme and with adaptors with four selective bases. This step was divided into three sections: sample obtaining, library construction and sequencing, and data analysis.

Sample obtaining and DNA extraction. Tissue and fin samples (depending on the species) of each aquaculture species (shrimp, oyster, almaco jack, n = 8 each) were obtained. DNA extractions were done with the Qiagen DNA extraction kit (QIAgen), following manufacturer instructions. DNA concentration of each sample was measured in a Nanodrop (Thermo Scientific) spectrophotometer and adjusted to 25 ng/µl. Pools of 4 DNA samples each were made: two from oyster (Cg-pool-1 and Cg-pool-2), two from almaco jack (Sr-pool-1 and Sr-pool-2), and four for shrimp (Lv-sel2-pool-1, Lv-sel2-pool-2, Lv-sel4-pool-1, and Lv-sel4-pool-2).

Library construction. 2bRAD libraries were prepared following the protocol of ^{Matz et al. (2018}) with modifications (see detailed protocol in mat.compl4). After the DNA was digested with the enzyme BcgI (defined in the previous step), the adaptors were ligated (adaptors without selective bases for the two oyster pools, adaptors without selective bases for the two almaco jack pools, one adaptor with two selective bases for two shrimp pools and two adaptors with two selective bases for the other two shrimp pools), PCR amplification and barcode adding, libraries pooling and purification of a 180 bp fragment.

Sequencing and data analysis.The library was prepared for sequencing (mat.compl4) in a MiSeq (Illumina) using a MiSeq Reagent Kit v3 (150-cycle) sequencing cell at CIBNOR. Sequencing reads quality was analyzed with the software FastQC version 0.11.2 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The Trimmomatic software (^{Bolger et al., 2014}) was used to delete low-quality and contaminating sequences (adaptors) and separate the sequences for each species according to the assigned barcodes. The sequences were aligned with BWA (^{Li & Durbin, 2009}) to each species’ reference genome, and the alignment percentage was calculated with Samtools (^{Li et al., 2009}). The bam files of each pool were used for SNP discovery with samtools, Bcftools, and Vcftools (Li et al., 2009; ^{Danecek et al., 2011}). For oyster and almaco jack, in silico analyses were done to estimate the number of remaining SNPs in the case of using adaptors with selective bases.

2bRAD technique validation

The validation of the 2bRAD technique and the SNP panel obtained for the genetic characterization and parentage analyses were done by analyzing 96 samples of almaco jackS. rivoliana, oysterC. gigas, red snapperLutjanus peru, rose spotted snapperL. guttatus(n = 8) and totoabaTotoaba macdonaldi(mat.compl5). Three known families of cultivated oyster (n = 24, 6 broodstock and 18 progenies), two almaco jack spawning groups from two breeding tanks (n = 45, 9 broodstock and 13 progenies from one culture tank, and eight broodstock and 15 progenies from the second tank) were used. In addition, wild samples from other aquaculture-important species, such as red snapperLutjanus peru(n = 8) from the state of Sinaloa (^{Reguera-Rouzaud et al., 2021}),L. guttatus(n = 8) from Panamá (^{Rivera-Lucero, 2020}), and totoabaTotoaba macdonaldi(n = 8) from the Gulf of California (^{Valenzuela-Quiñonez et al., 2016}), were included.

Library construction and sequencing were done with the procedure described in the supplementary material (mat.compl4), using the two adaptors with two selective bases each. The supplementary material (mat.compl6) contains the adaptor sequences and primers used for library construction and the position in a 96-well plate of the individual identification barcodes. Sequencing was done in a MiSeq (Illumina) using a MiSeq Reagent Kit v3 (150-cycle) sequencing cell. Sequence files were analyzed as described above throughout the demultiplexing, alignment, and SNP discovery process. In the case ofLutjanusandTotoaba, the sequence analysis was done with Stacks (^{Catchen et al., 2011}) because there are no complete genomes of these species in the GenBank. Stacks software allows the creation of pseudo-reference genomes using the sequences of each species. The SNP catalog creation step was using thedenovo_map.plfunction, with the following parameters: minimum depth to form a stack m = 3, maximum number of differences among stacks M = 3, and maximum number of differences among putative loci n =3.

SNP genotypes vcf files of each species were filtered with vcftools (^{Danecek et al., 2011}) with the following parameters: removal ofindels, retention of two alleles, SNPs represented in at least 80% of the individuals (max-miss= 0.8), and a minimum read quality of Q = 30. For oyster, themax-missparameter was 0.9, and only one SNP per fragment was retained. Loci with two alleles in at least one lot or varying relative to the reference genome were considered polymorphic. Thus, it was not necessary to calculate the minimum allele frequency of the minor allele (MAF). The filtered SNP files were converted into GenePop format with the software PGDSpider2 (^{Lischer & Excoffier, 2012}). They were imported from GenAlEx (^{Peakall & Smouse, 2012}) to analyze Hardy-Weinberg equilibrium and genetic variability (number of alleles per locus, heterozygosity) by species and lot (i.e., broodstock and progenies in oyster and almaco jack). Deviating loci from Hardy-Weinberg equilibrium (p< 0.05) were excluded. In oyster and almaco jack, parentage analyses were done with the software Colony (^{Jones & Wang, 2010}) and Cervus (^{Kalinowski et al., 2007}).

To implement a standardized genotyping analysis strategy for the subject species and the obtained SNPs, the number of individuals (Nej) that can be analyzed in a single sequencing run of a MiSeq (Illumina) platform was estimated with the following equation:

Where:

Nlps: Number of reads per sequencing cell; Mer: Reads depletion percentage derived from the PhiX control and low-quality reads; Prof: Desired read sequencing depth; Nsr: Approximate number of recognition sites of the 2B-type enzyme in the species’ genome after using selective bases.

In silico bioinformatic analysis for 2B-type enzyme selection

The average number of recognition sites for each enzyme among species corresponded, in ascending order, toCspCI, BcgI, AlfI, andBsaXI with 53,328, 116,354, 118,873, and 278,549, respectively (Table 1). A wide heterogeneity in the number of recognition sites per million bases for each enzyme and species is observed, suggesting that there is no random distribution of the nucleotidic bases in these species. For instance, the number of recognition sites per million bases of BcgI is 30% higher in shrimp than in oysters and 27% than in almaco jack; however, the contrary is true forBsaXI, AlfI, andCspCI. In the 2bRAD study of ^{Hernandez-Castro et al. (2017}) with the 702 Mb genome ofRhodnius prolixus(vector insect the Chagas disease), 146, 292, 100 sites per million bases for each of the enzymesBcgI, AlfI, andCspCI, were respectively obtained. This result indicates that there is not an associated number of recognition sites dependent on the taxonomy level.

Under the consideration that the genome size is similar between almaco jack and oyster and that the BcgI enzyme results in a similar number of recognition sites in both species (between 65 K and 70 K, respectively, Table 1), BcgI was the selected enzyme for subsequent analyses. Even though the CspCI enzyme would also be an adequate enzyme due to a more restricted variability, the random variation in the number of bases on both sides of the sequence (New England Biolabs, pers. comm.) derives in an undesired uncertainty during the digestion step.

The average reduction in the number of recognition sites as a function of the adaptors with one to four selective bases is relatively homogenous, independently of the species and the enzyme, with 22.9 %, 5.4 %, 1.3 %, and 0.3 %, respectively. In the case ofBsaXI, the reduction with five selective bases is 0.1 %. However, when the number of remaining sites, as a function of the selective bases, the number is higher when the bases are A or T (data not shown). This result implies that, depending on the application for which the genetic characterization will be applied, an adequate number of selective bases can be estimated. That is to say that if the study implies an association analysis of markers to a relevant genetic trait (known as GWAS), the number of bases could be 0, 1, or 2, but if the study is about parentage testing to simultaneously analyze a large number of individuals, 3 or 4 selective bases would be recommended.

The analysis with the shrimp database designed by ^{Perez-Enriquez et al. (2018}) to determine proportion variable sites (SNPs) within the fragments withBcgI recognition sites indicated that, in a random sample of 50 fragments, 65% showed at least one SNP. Thus, the in silico analysis indicates that using theBcgI enzyme without selective bases would result in 214,518 recognition sites and, thus, approximately 139,437 SNPs. This number is an overestimation compared to the 94,113 SNPs (after filtering) reported by ^{Lyu et al. (2021}). If, during the digestion step, the adaptors are added with two or four selective bases, the expected number of SNPs would be 7,530 and 418, respectively. Under the same parameters, usingBcgI and adaptors without selective bases in oyster and almaco jack would result in 42,000 to 45,000 SNPs for each species.

2bRAD technique standardization

BcgI digested libraries of shrimp, oyster, and almaco Jack pools resulted in a highly variable number of sequencing reads (Table 2), probably due to the preparation process rather than the species’ genome nature. The mean sequencing reads alignment against the shrimp’s, and almaco jack’s reference genomes (68 % and 95 %, respectively) is coincident with expectations according to the size of the genomes of those species in the GenBank (^{Seetharam, 2018}; ^{Zhang et al., 2019}). For oysters, a species for which the reference genome is complete, there were no differences in the alignment percentages between pools (^{Peñaloza et al., 2021}).

In shrimp, the number of filtered SNPs obtained using four selective bases in two adaptors is 2,874 (Table 2). The number of SNPs in oyster and almaco jack was similar before filtering (approx. 35,000), which were reduced by 25 % and 60 % after filtering, respectively (Table 2). The lower reduction in oysters compared to almaco jack might be due to the high heterozygosity reported in oysters (^{Peñaloza et al., 2021}). The in silico analysis for these two species indicates that using two adaptors with two selective bases each, 300 to 400 SNPs would be obtained. Depending on the objective of a low-density SNP panel (hundreds of SNPs), it was determined that the 2bRAD technique with four selective bases is adequate for oyster and almaco jack. Nevertheless, for shrimp, this technique is not viable due to the high number of markers (2,874); an alternative strategy for the development of a low-density panel (100 to 200 SNPs) in this species is the technique known as GT-Seq (Genotyping-in-thousands by sequencing; ^{Campbell et al., 2015}) that uses the genomic information and previously identified species-specific SNPs.

2bRAD technique validation

Sequencing of five species (almaco Jack, oyster, red snapper, rose spotted snapper, and totoaba) resulted in an average of 165,428 reads per individual and a high alignment percentage against references genomes (almaco jack and oyster) (Table 3). In almaco jack, 3,320 and 160 SNPs were obtained after alignment and filtering, respectively; in oysters, the number of SNPs was 3,731 and 114, respectively (Table 3). The comparative analysis of replicates (almaco jack, red snapper, and oyster) showed a low genotyping rate in the first two (0.6 % and 1.03 %, respectively). In oysters, the genotyping rate after SNP depuration was 2.1 %, although relatively low, which may have implications for the precision of parentage assignment analyses.

The analysis of deviation from Hardy-Weinberg equilibrium indicated that several loci were significant (p< 0.05), which were excluded for subsequent analyses (almaco jack: n = 1; oyster n = 3; totoaba n= 7). In red snapper and rose spotted snapper, the loci potentially in disequilibrium were not excluded due to a low sample size. Nevertheless, the presence of loci under selective pressure cannot be ruled out, as reported in the catfishRhamdia quelen(^{Ríos et al., 2020}), for which further analyses with adequate sample sizes for each species will be required.

The genetic variability in almaco JackS. rivolianaindicates a high number of alleles per locus (Na) in the broodstock. However, both the effective number of alleles (Ne) and heterozygosity (Ho, He) were medium to low (Table 4). Nevertheless, the negative value of F indicates a lack of inbreeding in this group. There are no data on the genetic variability in cultivated stocks ofS. rivoliana; however, in the closely related speciesS. lalandi, the genetic variability of cultivated populations estimated with microsatellites indicated medium He values (0.652 − 0.810) (^{Martínez Matus, 2016}).

When the variability between broodstock and progeny of S. rivoliana is compared, a reduction in Na, Ho, and He of approximately 7 % is observed, with a minimum variation in n_e and F (Table 4). This genetic variability reduction from one generation to the other is commonly seen in cultivated stocks because not all broodstock individuals participate during the spawning events (e.g., ^{Perez-Enriquez et al., 1999}). The parentage analysis in the almaco jack cultivation tanks based on 159 SNPs allowed the inference of the familiar structure of full- and half-sibs with high certainty. (the non-exclusion probability of parent pairs was 9 × 10^-9). In addition, it can be observed that, among ten females, three were represented in the progenies, and of nine males (one not genetically analyzed), six participated in the spawning events (Table 5). This result indicates a familiar structure in which more males than females are represented in the following generation, as in other species [e.g.,Seriola lalandi(^{Martínez Matus, 2016}),Lutjanus guttatus(^{Perez-Enriquez et al., 2020})]. The sequences of the 159 SNPs panel for S. rivoliana are shown in the supplementary material (mat.compl.7).

In oysters, following the criteria of ^{Perez-Enriquez et al. (2024}), the genetic variability showed medium to high values in the number of alleles per locus, but medium to low (and lower in progeny than in broodstock) in He (Table 4). However, there was an increment in Na from one generation to the other due to alleles present in the progeny but not in the broodstock; this can be a problem associated with a low sequencing depth in some of the parents in the broodstock or particular mutations, which has been reported in C. gigas as a necessary process to counteract against the effect of deleterious alleles (^{Durland et al., 2021}). The SNP panel had a non-exclusion probability for the parent pair of 1 × 10^-15. The parentage assessment in oysters with Colony using 114 SNPs showed an agreement of the genetic grouping at the fullsib level in most individuals as expected (Table 5). The sequences of the 114 SNPs panel forC. gigasare shown in the supplementary material (mat.compl.8).

The variability values for the other three species (red snapper, rose spotted snapper, totoaba) showed relatively high for n_e but medium for He (Table 4), which can be due to a small sample size and a lack of a reference genome. The availability of reference genomes for these species will be relevant for a better description of the diversity. Even though ^{Chambers et al. (2023}) recognize that the lack of reference genomes in non-model species and the small size of the 2bRAD restriction fragments are problematic for genotyping, the phylogenetic reconstruction and the proportion of missing genotypes derived from the technique generate comparable results, and higher repeatability than other ddRAD approaches. The genetic characterization of non-model aquatic species by 2bRAD with one adaptor with one selective base was reported by ^{Barbanti et al. (2020}). These authors reported that the number of samples for sequencing is balanced between the number of loci, the number of reads per individual derived from the sequencing platform (in that case, HiSeq, Illumina), and the available budget at a suggested depth coverage of 20×.

The present study demonstrated the feasibility of reducing the genome size of several aquaculture species by using two adaptors with two selective bases in each to obtain a sufficient number of loci in a MiSeq (Illumina) sequencing platform. In this regard, the estimated number of individuals that can be analyzed in a single sequencing run is approximatelyNej= 270 (see eq. 1), considering a sequencing cell of 25 million reads (Nlps), 20% of reads depletion (Mer), a sequencing depth of 20× (Prof), and an estimated number of recognition sites of the 2B enzyme due to the use of selective bases of 3700 (Nsr) (Table 3). A costing analysis of such characteristics to determine the applicability in aquaculture companies will be required; nevertheless, this estimation will still have high variability depending on the conditions of the sequencing laboratory and the acquisition costs of the supplies.