Method and device for tracing origin of manis pentadactyla based on multiplex PCR primer set

Abstract

The application discloses a Chinese pangolin geographical tracing method based on a multiplex PCR primer group and a device thereof, and comprises the following steps: step S1: obtaining whole genome sequencing data and geographical coordinate information of a source place of a reference sample; step S2: comparing each reference sample sequencing data in the reference sample data set to a published Chinese pangolin reference genome assembled to a chromosome level; step S3: eliminating samples with a genetic relationship in the reference sample data set; step S4: performing PCA analysis on a high-quality SNPs data set; step S5: screening an optimal combination by taking a plurality of gradient small-scale site combinations; step S6: generating a candidate primer pair of a target site; step S7: obtaining genotype information of an unknown source sample to be identified; and step S8: performing tracing analysis to predict the latitude and longitude information of the sample to be identified. In the variation detection and geographical tracing process, the method is relatively simple to operate.

Description

Technical Field

[0001] This application relates to the field of source tracing and identification technology for confiscated pangolins, and in particular to a method and apparatus for geographical source tracing of the Chinese pangolin based on multiplex PCR primer sets. Background Technology

[0002] Chinese pangolin ( Manis pentadactyla The Chinese pangolin (Pangolinus chinensis) is a mammal belonging to the order Pholidota and the family Manidae. Its geographical distribution is wide, ranging from Nepal in the southern Himalayas in the west to Taiwan in eastern my country, and extending from south of the Yangtze River to the northern part of the Indochina Peninsula. Currently, the Chinese pangolin is listed as a Class I protected wild animal in my country and is listed as critically endangered on the IUCN Red List of Threatened Species.

[0003] Precisely combating illegal trade at its source is the most direct and effective means of protecting the Chinese pangolin, leading to the need to develop traceability methods. Existing geographical traceability methods for the Chinese pangolin include trade record analysis and molecular genetic analysis. Molecular genetic analysis establishes a relationship between pangolin genetic information and geographical regions or coordinates, thereby predicting the origin of the sample based on specific molecular marker information.

[0004] Existing population studies have revealed a distinct population structure in the Chinese pangolin, with each population possessing a relatively independent geographical distribution area, providing a core foundation for its genetic origin tracing research. Furthermore, the research method of screening highly differentiated SNP markers among a small number of populations for origin tracing has been successfully applied in areas such as tracing the illegal trade in tree pangolins and tracking bird migrations, providing important insights for the genetic origin tracing research of the Chinese pangolin.

[0005] Identification of samples from the illegal trade of Chinese pangolins: Due to the large number of samples to be identified in smuggling cases, the lack of trade records, or the mixing of various pangolin scales, existing identification methods for identifying such scales have the following problems: (1) Limitations of trade record analysis method: This method is highly dependent on the completeness and accuracy of trade records. However, due to the long transportation routes, many transshipment links, and the concealment of trade information, it is difficult to guarantee the completeness and accuracy of trade records in illegal trade. Therefore, the results obtained by trade record analysis method are unreliable. At the same time, trade record analysis method is difficult to accurately locate the source of illegal trade samples. It can only point to a general area, which affects the accurate identification of the sample varieties. Although trade record analysis method is simple and direct, it is difficult to guarantee the reliability and accuracy of the results.

[0006] (2) Limitations of molecular genetic analysis methods: This method includes various analytical approaches. One type of method determines the population of the sample based on its genetic information and determines its geographical origin based on the distribution area of ​​the corresponding population. However, this method requires a large number of loci during detection, and the genetic information loci extracted from the sample cannot be completely matched with the population loci, resulting in low accuracy of geographical tracing results and low cost-effectiveness in practical applications. Another type of method, which screens highly differentiated sites among populations at the whole-genome level, requires prior population structure analysis to provide a population label for each sample. The population structure analysis, in turn, requires multiple analytical methods to cross-validate in order to determine the population corresponding to the sample to be tested. This method has a complex analysis process, is time-consuming, and the accurate interpretation of the results is complex and requires a lot of professional knowledge. Another type of method uses whole-genome variation site information and deep learning frameworks for geographic origin tracing. This method requires whole-genome sequencing data and a large amount of computing resources, takes a long time, and is costly and time-consuming in practical applications.

[0007] While various methods in molecular genetic analysis can be used to measure the accuracy and reliability of their results in multiple ways, they are inefficient and costly in practical applications. Neither can simultaneously meet the core law enforcement requirements of high accuracy, low time consumption, and low cost. As a result, grassroots law enforcement agencies face a double burden of manpower and finances when identifying large-scale smuggled armor plates, making it difficult to effectively combat illegal trade.

[0008] The information disclosed in the background section is intended only to enhance the understanding of the overall background of the invention and should not be construed as an admission or in any way implying that such information constitutes prior art known to those skilled in the art. Summary of the Invention

[0009] This application addresses the aforementioned technical problems by providing a method and apparatus for geographical tracing of the Chinese pangolin based on multiplex PCR primer sets.

[0010] This application provides a method for geographical tracing of the Chinese pangolin based on multiplex PCR primer sets, including the following steps: Step S1: Obtain the whole genome sequencing data and source geographical coordinates of the reference samples, and store the whole genome sequencing data of each sample and the source geographical coordinates one by one as the reference sample dataset; Step S2: Align the sequencing data of each reference sample in the reference sample dataset to the published Chinese pangolin reference genome assembled at the chromosome level. Use picard to remove redundant sequences generated by PCR during sequencing. Use Genome Analysis Toolkit v4.0 to obtain the SNP dataset. Use the HaplotypeCaller module to identify the SNPs of each individual corresponding to each sample. Merge the SNPs of each individual into a file to form a vcf file. Use vcftools in the VariantFiltration module to remove SNPs with minor allele frequencies less than 0.05. Remove SNPs at the population level where the deletion site at the same position covers more than 20% of the samples to ensure that the number of alleles in the population is 2. Remove SNPs with sequencing depths distributed at all sites that are less than 2.5% or greater than 97.5%. Complete the filtering of the obtained SNPs, remove false positive sites and multiple allele sites, and obtain the reference sample SNPs dataset. Step S3: Based on KING analysis, perform phylogenetic analysis on each sample in the reference sample SNPs dataset to obtain multiple sets of the same individual groups and individual datasets. Only retain the information of individuals with specific geographical origins or good sequencing quality in the same individual groups. After removing redundant data in the same individual groups, a clean dataset is obtained. Filter the clean dataset to obtain a high-quality SNPs dataset. Step S4: Use the smartPCA program in the Eigensoft v.4.2 software package to perform PCA analysis on the high-quality SNPs dataset to obtain the principal component values, the variance values ​​corresponding to each principal component value, and the contribution values ​​of each SNP to the principal components. Through the Tracy-Widom test, determine the principal components that have a high degree of data interpretation and are true in the high-quality SNPs dataset. Sort the SNP sites in the high-quality SNPs dataset in descending order of contribution. Select SNP sites and merge the selected SNP sites as the geographic origin SNP molecular marker dataset. The selection of SNP sites is as follows: determine the proportion of sites to be selected in each principal component based on the variance ratio of the principal components. Then, multiply the total number of sites to be selected by the corresponding proportion to obtain how many sites to select in each principal component. Select the corresponding number of top-contributing SNP sites under each principal component based on the number of sites to be selected. Step S5: Take several small-scale locus combinations of gradients, and input the genotype information of these locus combinations in the reference samples and the latitude and longitude coordinates of the geographical origin of each reference sample into the geo-origination software OriGen. Use leave-one-out cross-validation to evaluate the error distance of each group of loci for geo-origination. That is, in each round, the latitude and longitude coordinates of the geographical origin of a reference sample are set to unknown, and the predicted coordinates are obtained by the geo-origination software OriGen. The distance between the predicted coordinates and the true coordinates is calculated using the R package geodist, which is the error distance of geo-origination. The geo-origination error distance of all samples in each group is counted. The average geo-origination error distance of each group and the number of individuals with an error distance greater than 500km are compared. The group of loci in the geo-origination SNP molecular marker dataset corresponding to at least one of the above two indicators being the smallest is taken as the best SNP molecular marker set. Step S6: Use Primer3 to generate candidate primer pairs for the target sites, and then use MFEprimer-3.0 as the core quality control module. After multiple rounds of iterative screening and global checks, a pool of multiplex PCR primers that do not interfere with each other in the same reaction system and can be amplified efficiently and specifically is obtained. Step S7: Extract total DNA from the sample of unknown origin to be identified, use a customized multiplex PCR primer set to target and capture the sample and construct a second-generation sequencing library, use a high-throughput second-generation sequencing platform to sequence the library to obtain a second-generation sequencing dataset, and perform genotyping data processing on the second-generation sequencing dataset to obtain the genotype information of the optimal SNP molecular marker set for the sample to be tested. Step S8: Use the R package OriGen for source tracing analysis. Input the genotype information of the best SNP molecular marker set obtained in S7 into the OriGen software to obtain the predicted longitude and predicted latitude information of the geographical location of the sample to be identified, as well as the predicted posterior probability of each sample, as the prediction result output.

[0011] Preferably, the whole genome sequencing data and geographical coordinates of the reference sample are sourced from published literature and newly collected samples.

[0012] Preferably, the same group of individuals is screened and filtered according to the following steps: In PLINK v.2.0, the “make-king” command is used to iterate through and calculate the KING-robust coefficients between any two individuals in the reference sample SNPs dataset. Any two individuals with KING-robust coefficients greater than 0.354 are constructed into a group of identical individuals. In each group of identical individuals, the individual with the highest sequencing coverage is selected for subsequent analysis, and the remaining identical individuals are filtered out from the dataset.

[0013] Preferably, step S5 further includes: site combination simplification operation; The site combination simplification process includes the following steps: Step S511: Use the ancestor estimation function of the software HITSNP to construct training data based on the simulated hybrid population of purebred samples, and select the predicted probabilities output by LR-I, SVM-I, and RF-I in Classifier I to train ClassifierII; Step S512: Use Classifier II to distinguish whether the samples in the best SNP molecular marker set are purebreds or hybrids and infer the number of ancestors. Then, based on the threshold that the LR-I prediction probability is greater than 0.05, lock the specific ancestral species of the hybrid sample. Step S513: Using the minimum feature SNP subset screening function, with the feature SNP as the initial set, the sequence forward selection algorithm combined with the BagNB-III classifier is used to iteratively screen according to the quantitative gradient. The classification accuracy is used as the core evaluation criterion to obtain the minimum feature SNP subset that can accurately distinguish purebred populations as the best SNP molecular marker set for use in step S6.

[0014] Preferably, the process of constructing the next-generation sequencing library in step S7 includes: two rounds of PCR amplification and purification, followed by detection of the length distribution of library fragments to obtain a next-generation sequencing library that meets all quality control standards.

[0015] Preferably, the data processing of the second-generation sequencing dataset in step S7 includes the following steps: Step S721: The raw sequencing data is evaluated for quality using FastQC. Trimmomatic is used to remove adapter contamination sequences, low-quality bases, and reads that are too short to obtain high-quality sequencing data. Step S722: Use Burrows-Wheeler Aligner to align high-quality sequencing data to the Chinese pangolin reference genome. Use Samtools to sort, index, and convert the alignment results. Use Genome Analysis Toolkit v4.0 to obtain SNP data. Finally, use vcftools to obtain the genotype information of the optimal SNP molecular marker set for the sample to be tested.

[0016] Another aspect of this application provides an apparatus for operating the method described above, comprising: The sample information collection module is used to acquire the whole genome sequencing data and source geographical coordinates of the reference samples, and to store the whole genome sequencing data of each reference sample as a reference sample dataset after matching the source geographical coordinates one by one. The variant site detection and filtering module is used to align the sequencing data of each reference sample in the reference sample dataset to the published Chinese pangolin reference genome assembled at the chromosome level. It uses picard to remove redundant sequences generated by PCR during sequencing, uses Genome Analysis Toolkit v4.0 to obtain the SNP dataset, uses the HaplotypeCaller module to identify the SNPs corresponding to each individual sample, merges the SNPs of each individual into a single file to form a vcf file, uses vcftools in the VariantFiltration module to remove SNPs with a minor allele frequency less than 0.05, removes SNPs where the deletion site at the same position covers more than 20% of the samples at the population level, ensures that the allele count in the population is 2, and removes SNPs with sequencing depths below 2.5% or above 97.5% of all sites. This completes the filtering of the obtained SNPs, removing false positive sites and multiple allele sites, resulting in the reference sample SNPs dataset. The kinship identification module is used to perform kinship analysis on each sample in the reference sample SNPs dataset based on KING analysis, and obtain multiple sets of the same individual groups and individual datasets. Only individuals with specific geographical origins or good sequencing coverage in the same individual groups are retained. After removing redundant data in the same individual groups, a clean dataset is obtained. After filtering the clean dataset, a high-quality SNPs dataset is obtained. The PCA module is used to perform PCA analysis on high-quality SNP datasets using the smartPCA program in the Eigensoft v.4.2 software package. It obtains the principal component values ​​for each sample, the variance values ​​corresponding to each principal component value, and the contribution values ​​of each SNP to the principal components, using the Tracy-Widom algorithm. The process involves verifying and identifying principal components with high explanatory power and true accuracy in high-quality SNP datasets. SNP sites in these datasets are then sorted in descending order of contribution. Selected SNP sites are merged to form a geographic origin SNP molecular marker dataset. Each group of sites in the geographic origin SNP molecular marker dataset is compared against the mean error distance and the number of samples with an error distance greater than 500 km in the geographic source dataset. The dataset with the smallest mean error distance and the smallest number of samples with an error distance greater than 500 km is selected as the optimal site set. SNP site selection is based on the proportion of sites to be selected in each principal component according to the principal component variance ratio. The total number of sites to be selected is multiplied by the corresponding proportion to obtain the number of sites to be selected in each principal component. The corresponding number of SNP sites are then selected based on this number. The module for obtaining the optimal SNP molecular marker set takes several small-scale locus combinations at different gradients. The genotype information of these locus combinations from the reference samples and the latitude and longitude coordinates of the geographical origin of each reference sample are input into the geo-origination software OriGen. Leave-one-out cross-validation is used to evaluate the error distance of each group of loci for geo-origination. That is, in each round, the latitude and longitude coordinates of the geographical origin of a reference sample are set to unknown. The geo-origination software OriGen predicts the coordinates, and the R package geodist is used to calculate the distance between the predicted coordinates and the true coordinates, i.e., the error distance of geo-origination. The geo-origination error distance of all samples with different locus combinations is counted. The average geo-origination error distance of each group and the number of individuals with an error distance greater than 500km are compared. The group of loci in the geo-origination SNP molecular marker dataset with at least one of the above two indicators being the smallest is taken as the optimal SNP molecular marker set. The multiplex PCR primer pool acquisition module is used to generate candidate primer pairs for target sites using Primer3. Then, MFEprimer-3.0 is used as the core quality control module. Through multiple rounds of iterative screening and global checks, multiplex PCR primer pools that do not interfere with each other in the same reaction system and can be amplified efficiently and specifically are obtained. The high-quality sequencing data acquisition and genotyping module is used to extract total DNA from samples of unknown origin to be identified, use a customized multiplex PCR primer set to target and capture quality control samples and construct a second-generation sequencing library for the unknown sample to be identified, use a high-throughput sequencing platform to sequence the second-generation sequencing library to obtain a second-generation sequencing dataset, process the second-generation sequencing dataset to obtain high-quality sequencing data, and obtain the genotype information of the optimal SNP molecular marker set for the sample to be tested through genotyping data processing. The prediction result output module is used to perform source tracing analysis using the R package OriGen. It inputs the genotype information of the best SNP molecular marker set of the sample to be identified into the OriGen software to obtain the predicted longitude and predicted latitude information of the geographical location of the sample to be identified, as well as the predicted posterior probability of each sample, as the prediction result output. Genotype information of the target SNP molecular marker set of the reference sample was obtained by using vcftools to extract the genotype information of the optimal SNP molecular marker set from the whole genome SNP dataset of the reference sample.

[0017] The beneficial effects that this application can produce include: 1) Achieved low-cost and simple-to-operate geographical traceability of the Chinese pangolin: This method does not require whole-genome sequencing of the genome information of Chinese pangolin samples. It only requires multiplex PCR amplification of the target fragment to obtain the genotype information of the required loci, effectively reducing the detection cost. At the same time, because this method requires fewer loci, it is simple to operate in the process of variation detection and geographical tracing, requires less computing power and consumes less time during identification. 2) Site screening is simple Currently, the methods for selecting loci in the geographical tracing of the Chinese pangolin are generally quite complex. For example, some methods require first analyzing the population structure of the Chinese pangolin using multiple methods to provide a population label for each sample before subsequent locus selection can proceed. This method does not require such pre-analysis as population structure analysis; once the SNP dataset is obtained, target loci can be selected directly.

[0018] 3) High traceability accuracy The optimal locus set obtained by this method had an average geographical origination error distance of 163.92 km and 6 individuals with an error distance greater than 500 km. In contrast, the origination method for whole-genome SNPs using the deep learning framework Locator had an average error distance of 175.43 km and 9 individuals with an error distance greater than 500 km. Considering all three indicators, the geographical origination performance of the small-scale SNPs screened for the Chinese pangolin is close to that of the whole-genome SNPs. Attached Figure Description

[0019] Figure 1 This is a graph showing the error distance results corresponding to different gradient site sets selected in the embodiments of this application; Figure 1 The horizontal axis represents the number of SNP sites to be selected in step S5. Figure 1 The left vertical axis represents the average error distance between the predicted coordinates and the true coordinates obtained by using the leave-one-out method with the selected locations. Figure 1 The right vertical axis represents the number of individuals whose error distance between the predicted coordinates and the true coordinates obtained by using the leave-one-out method at the selected sites is greater than 500 km. Figure 2 This is a distribution map of the geographic origination error distance results for the optimal SNP combination in the embodiments of this application.

[0020] Figure 3 A flowchart illustrating the method provided in this application; Figure 4 This application provides a schematic diagram of the module connections of the device; Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0022] Therefore, the following detailed description of the embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0023] Technical means not detailed in this application and not used to solve the technical problems of this application are all set according to common general knowledge in the field, and multiple common general knowledge setting methods can be implemented.

[0024] See Figures 1-3 The geographical origin tracing method for the Chinese pangolin based on multiplex PCR primer sets provided in this application includes the following steps: Step S1: Obtain the whole genome sequencing data and source geographical coordinates of the reference samples, and store the whole genome sequencing data of each sample and the source geographical coordinates one by one as the reference sample dataset; Whole-genome sequencing data of samples covering the distribution range of the Chinese pangolin were obtained. The source of the gene sequencing data was the download of sequencing data of samples from published literature. After downloading, the MD5 value was checked to ensure the integrity of the downloaded files. For newly collected samples, tissue from the samples is used to extract DNA. PCR amplification and bidirectional Sanger sequencing are performed on specific mitochondrial gene fragments. The DNA is then compared with the GenBank database to confirm that the sample is from the Chinese pangolin. Library construction and sequencing are then performed using next-generation sequencing methods to obtain the sequencing data. The specific sequencing method used can follow existing methods or, if necessary, be outsourced to a sequencing company. The geographical coordinates of the sample's origin are also obtained.

[0025] While acquiring the whole genome sequencing data of each reference sample, the geographical origin coordinates of each sample are also acquired. The whole genome sequencing data of each reference sample are then associated with the geographical origin coordinates of that sample and stored accordingly.

[0026] Step S2: Detection of variant sites Sequencing data from each reference sample in the reference sample dataset were aligned to the published chromosome-level assembly of the Chinese pangolin reference genome. Redundant sequences generated during sequencing due to PCR were removed using picard. Single nucleotide polymorphism (SNP) datasets were obtained using Genome Analysis Toolkit (GATK) v4.0. The HaplotypeCaller module was used to identify the SNPs for each individual sample. The SNPs for each individual were merged into a single file to form a vcf file. The vcftools software was used to remove SNPs with a minor allele frequency less than 0.05, SNPs with deletions at the same position covering more than 20% of samples at the population level (ensuring a population allele count of 2), and SNPs with sequencing depths below 2.5% or above 97.5% of all sites. This filtering process removed false positives and multiple alleles, resulting in the reference sample SNP dataset.

[0027] Step S3: Kinship Testing Kinship analysis was performed on the reference sample SNPs dataset based on KING analysis. In PLINK v.2.0, the "make-king" command was used to calculate the KING-robust coefficients between any two individuals in the reference sample SNPs dataset. This value reflects the proportion of identical SNPs between individuals. Individual pairs with a KING-robust coefficient greater than 0.354 were considered identical individuals. Individuals with whole genome sequences corresponding to KING-robust coefficients greater than 0.354 were grouped into identical individual groups. Only individuals with specific geographic origins or good sequencing quality within these groups were retained. After removing redundant data from identical individual groups, a cleaned dataset was obtained. This cleaned dataset was then filtered, and the resulting high-quality SNPs dataset was used for subsequent analysis. A geographic origin dataset was also obtained for analyzing error distance.

[0028] Step S4: Screening of SNP molecular markers for geographical origin tracing of the Chinese pangolin based on principal component analysis (PCA) First, the high-quality SNPs dataset obtained in step S3 is subjected to linkage disequilibrium removal. The smartPCA program in the Eigensoft v.4.2 software package is used to perform PCA analysis on the high-quality SNPs dataset. After the analysis is completed, the principal component value, the variance value corresponding to each principal component value, and the contribution value of each SNP to the principal component are obtained.

[0029] The Tracy-Widom test was used to identify principal components that provide significant and accurate explanations for the high-quality SNPs dataset. Within these principal components, SNP sites in the high-quality SNPs dataset were sorted in descending order of contribution. Selected SNP sites were then merged to form a geographic origin SNP molecular marker dataset. The proportion of sites selected in each principal component represents the proportion of variance values ​​corresponding to that principal component.

[0030] The error analysis results of the optimal site set obtained in this embodiment relative to the real geographic coordinates from the IUCN Chinese pangolin distribution area dataset are as follows: Figure 2 As shown in the figure, the color of the circle coordinates and the length of the line segment represent the magnitude of the error distance. The line segment connects the predicted coordinates and the true coordinates, and its length is the magnitude of the error distance. The color of the circle changes from blue to red, indicating that the error distance increases from small to large. The size of the circle represents the number of samples corresponding to that sample point. The figure demonstrates that the method provided in this application can accurately identify the coordinates of most unknown samples.

[0031] Step S5: Selection of the optimal combination of SNP molecular markers Following the method in step S4, several small-scale locus combinations of gradients are selected. The genotype information of these locus combinations in the reference samples and the latitude and longitude coordinates of the geographical origin of each reference sample are input into the geo-origination software OriGen. Leave-one-out cross-validation is used to evaluate the error distance of each group of loci for geo-origination. That is, in each round, the latitude and longitude coordinates of the geographical origin of a reference sample are set to unknown. The predicted coordinates are obtained by the geo-origination software OriGen. The distance between the predicted coordinates and the true coordinates is calculated using the R package geodist, which is the error distance of geo-origination.

[0032] The geographic origination error distance of all samples in each group was calculated. The average geographic origination error distance and the number of individuals with an error distance greater than 500 km were compared in each group. The group of sites in the geographic origination SNP molecular marker dataset corresponding to at least one of the above two indicators being the smallest was taken as the best SNP molecular marker set to be adopted in the final selection.

[0033] Step S6: Design of the optimal SNP molecular marker set for multiplex PCR primers For the optimal SNP molecular marker set selected in step S5, candidate primer pairs for the target sites are first generated using Primer3. Then, MFEprimer-3.0 is used as the core quality control module. After multiple rounds of iterative screening and global checks, a multiplex PCR primer pool with no mutual interference in the same reaction system and capable of efficient and specific amplification is obtained in step S7.

[0034] Step S7: Genotyping of unknown samples Total DNA was extracted from samples of unknown origin to be identified. A customized multiplex PCR primer set was used to target and capture the DNA from quality control samples, and next-generation sequencing libraries for the unknown samples were constructed. The library construction process included two rounds of PCR amplification and purification, followed by fragment length distribution analysis to obtain next-generation sequencing libraries that met all quality control standards. These libraries were then sequenced using a high-throughput sequencing platform to obtain the next-generation sequencing dataset.

[0035] After sequencing, the raw sequencing data was assessed for quality using FastQC. Trimmomatic was used to remove adapter contamination sequences, low-quality bases, and excessively short reads. Subsequently, the high-quality sequencing data was aligned to the Chinese pangolin reference genome (from CNGBdb, BioProject: PRJCA020583) using Burrows-Wheeler Aligner. Samtools was then used to sort, index, and convert the alignment results to obtain high-quality sequencing data.

[0036] To address the technical characteristics of multiplex PCR sequencing, residual multiplex primer sequences in the aligned sequences were specifically removed to eliminate their interference with variant detection. Based on this, the sequencing depth and coverage of each target region were statistically analyzed to evaluate the validity of the sequencing data.

[0037] GATK, Samtools, and VarScan were used in combination to detect single nucleotide polymorphisms (SNPs) in the samples. The obtained SNPs were stored in the dataset, and the genotype information of the target SNP sites was finally obtained.

[0038] Step S8: Geographical Origin Tracing of Unknown Samples Source tracing analysis was performed using the R package OriGen. Step S7 obtained the genotype information of the target SNPs in the sample to be identified. Simultaneously, the genotype information of the target SNPs obtained in step S5 of the reference sample and the geographical origin latitude and longitude coordinates were integrated as input files to obtain the predicted longitude and latitude of the sample's origin, as well as the predicted posterior probability for each sample, which were output as the prediction results. The reliability of the predicted location was determined based on the predicted posterior probability. The genotype information of the target SNPs was obtained by extracting the genotype information of the optimal SNP molecular marker set from the reference sample SNP dataset using vcftools. Example

[0039] This section provides a detailed explanation of the technical solution and its resulting technical effects using specific experimental data and operational details.

[0040] Step S1: Collection of whole genome sequencing data and geographical origin coordinates of the samples A total of 302 Chinese pangolin samples with sequencing data were collected. The number of samples from published articles were as follows: 22 from Hu et al. (2020), 25 from Hu et al. (2026), 15 from Wang et al. (2022), 56 from Wei et al. (2024), 35 from Zhai et al. (2025), and 6 from Koju et al. (2025). Sequencing data from the published articles were downloaded from NCBI. MD5 values ​​were checked to ensure the integrity of the downloaded files. The geographical coordinates of the source of each sequencing data point in the published literature were obtained. The literature sources for the above data are as follows: Hu, JY, Hao, ZQ, Frantz, L., Wu, SF, Chen, W., Jiang, YF, Wu, H., Kuang, WM, Li, HP, Zhang, YP, & Yu, L. (2020). Genomicconsequences of population decline in critically endangered pangolins and their demographic histories. National Science Review, 7(4), 798-814. https: / / doi.org / 10.1093 / nsr / nwaa031 Hu, JY, Jiang, Y., Li, S., Ying, TT, & Yu, L. (2026). Pangolinmuseum genomics reveal the dynamic genetic consequences and extinction risk of the critically endangered Chinese pangolin. Molecular Biology and Evolution, 43(5). https: / / doi.org / 10.1093 / molbev / msag099 Koju, NP, Zeng, Z., Zhang, G., Huang, X., Yao, Z., Wang, X.,Hawkins, MTR, Hinckley, A., Flores, MFC, Guo, C., Li, J., Maharjan,D., Huang, L., Yu, W., Leng, L., He, K. A., & Feijóa, Y., Hua, Y. (2025).Revalidation of Manis aurita, based on integrative morphological and genomicevidence. bioRxiv, 2025.2007.2005.663294. https: / / doi.org / 10.1101 / 2025.07.05.663294 Wang, Q., Lan, T., Li, H., Sahu, SK, Shi, M., Zhu, Y., Han, L.,Yang, S., Li, Q., Zhang, L., Deng, Z., Liu, H.,&Hua, Y. (2022). Whole-genomeresequencing of Chinese pangolins reveals a population structure and providesinsights into their conservation. Communications Biology, 5(1), 821. https: / / doi.org / 10.1038 / s42003-022-03757-3 Wei, SC, Fan, HZ, Zhou, WL, Huang, GP, Hua, Y., Wu, SB, Wei, X., Chen, YT, Tan, XY,&Wei, FW (2024). Conservationgenomics of the critically endangered Chinese pangolin. SCIENCE CHINA-LIFESCIENCES, 67(10), 2051-2061. https: / / doi.org / 10.1007 / s11427-023-2540-y Zhai, TY, Wang, JC, Zhan, GQ, Hu, JY, & Yu, L. (2025).Population genomic landscapes and insights for conservation of the critically endangered island-endemic Chinese pangolin in Taiwan [; Early Access].SCIENCE CHINA-LIFE SCIENCES, 16. https: / / doi.org / 10.1007 / s11427-024-2904-6 The newly collected samples were analyzed using the MagPure Forensic DNA Extraction Kit (Catalog No.: D6359D) from Guangzhou Magen Biotechnology Co., Ltd. After species identification as the Chinese pangolin, BGI Genomics and Berry Genomics Sequencing Co., Ltd. were commissioned to perform genome resequencing on the obtained samples. Using BGI's DNBseq sequencing platform and Berry's Illumina platform, 150 bp insert sequencing libraries were constructed, and paired-end sequencing was performed to obtain sequencing data from the Chinese pangolin scales. The geographical coordinates of the origin of the newly collected samples identified as belonging to the Chinese pangolin were also obtained.

[0041] The raw data was filtered using the FASTP software to remove adapter sequences and low-quality sequences (more than 40% of the bases with a quality value of Q≤15 and more than 5 N bases) to obtain high-quality clean reads for subsequent analysis.

[0042] Of the newly collected samples, 105 had a clearly defined geographical origin, covering the distribution area of ​​the Chinese pangolin. A clearly defined geographical origin is defined as a description of the origin at the city level or below, involving six countries. Among these, 33 published samples have clearly defined latitude and longitude coordinates of their origin. The sequencing data of each reference sample were stored in a one-to-one correspondence with the geographical coordinate information of its origin.

[0043] Step S2: Detection and Filtering of Variant Sites First, using the arn algorithm in BWA software (Li & Durbin, 2009), high-quality clean reads were aligned to the published Chinese pangolin reference genome assembled at the chromosome level (from CNGBdb, BioProject:PRJCA020583). Redundant sequences generated during sequencing due to PCR were removed using picard, and alignment-completed BAM files were generated using samtools (Li et al., 2009).

[0044] Single nucleotide polymorphism (SNP) datasets were obtained using Genome Analysis Toolkit (GATK) v4.0 (Mckenna et al., 2010). The HaplotypeCaller module was used to identify SNPs for each individual, generating a gvcf file for each individual. The CombineGVCFs module was used to merge the gvcf files of all individuals into a vcf file. A high-quality SNP dataset was obtained using the VariantFiltration module. Vcftools (Danecek et al., 2011) was used to remove SNPs with minor allele frequencies less than 0.05, SNPs with deletions at the same location covering more than 20% of samples at the population level, ensuring an allele count of 2 in the population, and SNPs with sequencing depths below 2.5% or above 97.5% across all sites were removed.

[0045] Step S3: Kinship Testing In PLINK v.2.0 (Purcell et al., 2007), the "make-king" command was used to calculate the KING-robust coefficients between individuals for kinship identification. A total of 39 pairs of samples had KING-robust coefficients greater than 0.354, and were considered to be samples of the same individual, involving a total of 64 samples. Samples without geographical origin or with low sequencing quality were removed, resulting in 34 samples being discarded. After removing identical samples within the same individual group, a total of 268 Chinese pangolin samples remained as a clean dataset, of which 96 samples had a clear geographical origin. After re-filtering and cleaning the dataset, a Chinese pangolin dataset containing 12,733,581 bp of high-quality autosomal SNPs was obtained.

[0046] In addition, a geographic origin dataset containing 9,036,169 bp of high-quality autosomal SNPs with clear geographic origins was obtained, which was used to evaluate the error distance in subsequent analyses.

[0047] Step S4: Screening of SNP molecular markers for geographical origin tracing of the Chinese pangolin based on principal component analysis (PCA) First, principal component analysis (PCA) was performed using the smartPCA program in the Eigensoft v.4.2 software package (N. Patterson et al., 2006). At the same time, by setting the parameter snpweightoutname, the contribution of each SNP to all output principal components was output.

[0048] Then, the Tracy-Widom test was used to identify principal components that provided significant and accurate explanations for the high-quality SNPs dataset. Within these principal components, each SNP in the high-quality SNPs dataset was sorted in descending order of contribution. Selected SNPs were then merged to form a geographic origin SNP molecular marker dataset of highly differentiated sites. The proportion of sites selected in each principal component represents the proportion of the variance value corresponding to that principal component. The selection of SNPs involved determining the proportion of sites to be selected in each principal component based on the variance ratio of the principal components. Then, the total number of sites to be selected was multiplied by the corresponding proportion to obtain the number of sites to be selected from each principal component. Finally, the corresponding number of SNPs was selected based on this number of sites. A gradient was taken every 50 loci, and the above method was used to screen five small-scale loci with a total number of 50 to 250 loci. Then, the genotype of each locus corresponding to the reference sample and the latitude and longitude coordinates of the reference sample's origin were input into the geo-originating software OriGen (Ranola et al., 2014), and the leave-one-out method was used to evaluate the error distance of each locus for geo-originating purposes.

[0049] Step S5: Using the geo-originating software OriGen, in each round, the latitude and longitude coordinates of the geographical origin of a reference sample in the geo-originating SNP molecular marker dataset are set to unknown. The predicted coordinates are obtained by the geo-originating software OriGen. The distance between the predicted coordinates and the true coordinates is calculated using the R package geodist, which is the error distance of geo-originating.

[0050] The geographic origination error distance of all samples in each group of the geographic origination SNP molecular marker dataset is statistically analyzed. The average geographic origination error distance of each group and the number of individuals with an error distance greater than 500 km are compared. The group of sites in the geographic origination SNP molecular marker dataset that have at least one of the above two indicators as the minimum are taken as the optimal SNP molecular marker set.

[0051] In this embodiment, comparing these three indicators, it was found that the 150 sites selected by the PCA method performed best in geographic tracing, with the lowest average value of 162.76 km among all groups. The number of individuals with an error distance greater than 500 km was only 5, indicating that the prediction results had fewer extreme error distances. The results obtained in this embodiment are as follows: Figure 1 As shown in the figure, among the three methods in this embodiment, the PCA method performs best under all point gradients; and within the PCA method, the 150-point geolocation method performs best.

[0052] Step S6: Design of the optimal SNP molecular marker combination multiplex PCR primer set First, candidate primer pairs were generated for each target sequence in the optimal SNP molecular marker set using Primer3. MFEprimer-3.0 was used as the core quality control module to comprehensively screen the candidate primer pairs based on their specificity, secondary structure (hairpins, self-dimers), SNP distribution at binding sites, and non-specific amplification performance. Specificity indicators ensured that the target fragment corresponding to the primer could be correlated with a unique fragment on the reference genome, i.e., they did not possess non-specific amplification performance. Secondary structure indicators themselves could not form secondary structures such as hairpins or self-dimers. Regarding the SNP distribution at binding sites, both the primer itself and the 10-20 bp regions flanking the primer binding region were avoided, preventing the presence of SNPs and indels.

[0053] After screening primers that meet the requirements of singleton PCR using the decision module, cross-compatibility verification was further performed with the selected primers to exclude combinations with cross-dimers and combined non-specific amplification. After multiple rounds of iterative screening and global final inspection, 6 sites that did not meet the primer design requirements were excluded, and finally a pool of 144 pairs of multiplex PCR primers with no mutual interference in the same reaction system and capable of efficient and specific amplification were obtained.

[0054] Step S7: Genotyping of unknown samples First, the MagPure Forensic DNA Extraction Kit (Catalog No.: D6359D) from Guangzhou Magen Biotechnology Co., Ltd. was used to extract total DNA from the unknown sample.

[0055] Subsequently, customized multiplex PCR primers were used for targeted capture and construction of next-generation sequencing libraries for unknown samples. First, genomic DNA samples underwent dual quality control: agarose gel electrophoresis was used to assess DNA degradation and contamination by RNA, proteins, and other impurities; Qubit quantitative PCR was used to accurately determine DNA concentration, retaining only OD. 260 / 280 Samples with a value between 1.8 and 2.0 and a total DNA volume ≥40 ng were used for subsequent library construction. Customization here refers to the primers required for genotyping at 144 SNP loci.

[0056] The library construction process consists of two rounds of PCR amplification and purification: the first round uses a primer pool (Primer pool T1) containing multiple upstream and downstream primers to simultaneously and specifically amplify multiple target regions on the genome (reaction time 2 h), and the amplification products are purified by magnetic beads in the first round (35 min) to remove free primers and non-specific fragments. A second round of PCR reaction was then performed (25 min), with Illumina paired-end sequencing adapters (PE1.0 / PE 2.0) and sample-specific tag sequences introduced at both ends of the amplicon. After a second round of magnetic bead purification (35 min), library construction was completed (commissioned to a gene company).

[0057] Library quality control included: quantifying library concentration using the Qubit dsDNA HS Assay Kit, requiring a acceptable library concentration in the range of 10–50 ng / μl; and detecting the fragment length distribution of the library using the Agilent 2100 Bioanalyzer system (Agilent DNA 1000 Kit), confirming that the fragments were concentrated in the 300–450 bp range. Finally, a next-generation sequencing library meeting all quality control standards was obtained and sequenced using a high-throughput sequencing platform.

[0058] After sequencing, this embodiment employs a standardized bioinformatics workflow to systematically analyze the multiplex PCR targeted sequencing data: After obtaining the raw sequencing data, data quality assessment was first conducted using FastQC (Andrews, 2014), and Trimmomatic (Bolger et al., 2014) was used to remove adapter contamination sequences, low-quality bases, and excessively short reads to obtain high-quality sequencing data (clean reads); subsequently, the high-quality sequencing data was aligned to the Chinese pangolin reference genome (from CNGBdb, BioProject: PRJCA020583) using Burrows-Wheeler Aligner (BWA) (Li & Durbin, 2009), and the alignment results were sorted, indexed, and format converted using Samtools (Li et al., 2009).

[0059] To address the technical limitations of multiplex PCR sequencing, residual multiplex primer sequences in the aligned sequences were specifically removed to eliminate their interference with variant detection. Based on this, the sequencing depth and coverage of each target region were statistically analyzed to evaluate the validity of the sequencing data. Simultaneously, three tools—GATK (Mckenna et al., 2010), Samtools (Li et al., 2009), and VarScan (Koboldt et al., 2012)—were used in combination to detect single nucleotide polymorphisms (SNPs) in the samples. The obtained SNP data were stored in the dataset, ultimately yielding the genotype information of the target SNP loci. SNP detection involved using these three tools to identify SNP variants from the alignment results, and then extracting the required SNP genotype information based on the locus coordinates.

[0060] Step S8: Geographical Origin Tracing of Unknown Samples In step S7, the genotype information of the target SNPs in the sample to be identified is obtained. Simultaneously, the genotype information of the target SNPs in the reference sample and the geographical origin latitude and longitude coordinates of the reference sample are integrated as input files. Source tracing analysis is performed using the R package OriGen. The genotype information of the target SNPs is obtained in step S5, and the genotype information of the optimal SNP molecular marker set is extracted from the reference sample SNP dataset using vcftools.

[0061] OriGen is an R package for geolocation. It divides the target study area into several raster pixels, constructs an allele frequency surface for each SNP, calculates the posterior probability based on Bayes' theorem, locates the geographic origin of the target individual, and obtains the raster pixel with the highest posterior probability, which is the predicted geographic origin of the individual.

[0062] Use the FitOriGenModelFindUnknowns function in this package to perform geographic prediction of the sample to be identified, where the parameter RhoParameter is obtained by the function FindRhoParameterCrossValidation.

[0063] After the analysis is completed, the results include two columns: predicted longitude (predLong) and predicted latitude (predLat). The results also report the predicted posterior probability (highprob) for each sample. This probability can be used to determine the reliability of the prediction; the higher the value, the more reliable the prediction. If the predicted posterior probability is below 0.1, the result is considered invalid.

[0064] In this embodiment, the results obtained after identifying the sample are: 1. Genotype information corresponding to 144 loci; 2. Latitude and longitude coordinates obtained from geographical tracing prediction and the predicted posterior probability.

[0065] Another aspect of this application provides an apparatus for implementing the above method, comprising: a sample information collection module for acquiring whole genome sequencing data and source geographical coordinate information of reference samples, and storing the whole genome sequencing data of each reference sample in a one-to-one correspondence with the source geographical coordinate information as a reference sample dataset; The variant site detection and filtering module is used to align the sequencing data of each reference sample in the reference sample dataset to the published Chinese pangolin reference genome assembled at the chromosome level. It uses picard to remove redundant sequences generated by PCR during sequencing, uses Genome Analysis Toolkit v4.0 to obtain the SNP dataset, uses the HaplotypeCaller module to identify the SNPs corresponding to each individual sample, merges the SNPs of each individual into a single file to form a vcf file, uses vcftools in the VariantFiltration module to remove SNPs with a minor allele frequency less than 0.05, removes SNPs where the deletion site at the same position covers more than 20% of the samples at the population level, ensures that the allele count in the population is 2, and removes SNPs with sequencing depths below 2.5% or above 97.5% of all sites. This completes the filtering of the obtained SNPs, removing false positive sites and multiple allele sites, resulting in the reference sample SNPs dataset. The kinship identification module is used to perform kinship analysis on each sample in the reference sample SNPs dataset based on KING analysis, and obtain multiple sets of the same individual groups and individual datasets. Only individuals with specific geographical origins or good sequencing coverage in the same individual groups are retained. After removing redundant data in the same individual groups, a clean dataset is obtained. After filtering the clean dataset, a high-quality SNPs dataset is obtained. The PCA module is used to perform PCA analysis on high-quality SNP datasets using the smartPCA program in the Eigensoft v.4.2 software package. It obtains the principal component values ​​for each sample, the variance values ​​corresponding to each principal component value, and the contribution values ​​of each SNP to the principal components, using the Tracy-Widom algorithm. The process involves verifying and identifying principal components with high explanatory power and true accuracy in high-quality SNP datasets. SNP sites in these datasets are then sorted in descending order of contribution. Selected SNP sites are merged to form a geographic origin SNP molecular marker dataset. Each group of sites in the geographic origin SNP molecular marker dataset is compared against the mean error distance and the number of samples with an error distance greater than 500 km in the geographic source dataset. The dataset with the smallest mean error distance and the smallest number of samples with an error distance greater than 500 km is selected as the optimal site set. SNP site selection is based on the proportion of sites to be selected in each principal component according to the principal component variance ratio. The total number of sites to be selected is multiplied by the corresponding proportion to obtain the number of sites to be selected in each principal component. The corresponding number of SNP sites are then selected based on this number. The module for obtaining the optimal SNP molecular marker set takes several small-scale locus combinations at different gradients. The genotype information of these locus combinations from the reference samples and the latitude and longitude coordinates of the geographical origin of each reference sample are input into the geo-origination software OriGen. Leave-one-out cross-validation is used to evaluate the error distance of each group of loci for geo-origination. That is, in each round, the latitude and longitude coordinates of the geographical origin of a reference sample are set to unknown. The geo-origination software OriGen predicts the coordinates, and the R package geodist is used to calculate the distance between the predicted coordinates and the true coordinates, i.e., the error distance of geo-origination. The geo-origination error distance of all samples with different locus combinations is counted. The average geo-origination error distance of each group and the number of individuals with an error distance greater than 500km are compared. The group of loci in the geo-origination SNP molecular marker dataset with at least one of the above two indicators being the smallest is taken as the optimal SNP molecular marker set. The multiplex PCR primer pool acquisition module is used to generate candidate primer pairs for target sites using Primer3. Then, MFEprimer-3.0 is used as the core quality control module. Through multiple rounds of iterative screening and global checks, multiplex PCR primer pools that do not interfere with each other in the same reaction system and can be amplified efficiently and specifically are obtained. The high-quality sequencing data acquisition and genotyping module is used to extract total DNA from samples of unknown origin to be identified, use a customized multiplex PCR primer set to target and capture quality control samples and construct a second-generation sequencing library for the unknown sample to be identified, use a high-throughput sequencing platform to sequence the second-generation sequencing library to obtain a second-generation sequencing dataset, process the second-generation sequencing dataset to obtain high-quality sequencing data, and obtain the genotype information of the optimal SNP molecular marker set for the sample to be tested through genotyping data processing. The prediction result output module is used to perform source tracing analysis using the R package OriGen. It inputs the genotype information of the best SNP molecular marker set of the sample to be identified into the OriGen software to obtain the predicted longitude and predicted latitude information of the geographical location of the sample to be identified, as well as the predicted posterior probability of each sample, as the prediction result output. Genotype information of the target SNP molecular marker set of the reference sample was obtained by using vcftools to extract the genotype information of the optimal SNP molecular marker set from the whole genome SNP dataset of the reference sample.

[0066] Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for geographical tracing of the Chinese pangolin based on multiplex PCR primer sets, characterized in that, Includes the following steps: Step S1: Obtain the whole genome sequencing data and source geographical coordinates of the reference samples, and store the whole genome sequencing data of each sample and the source geographical coordinates one by one as the reference sample dataset; Step S2: Align the sequencing data of each reference sample in the reference sample dataset to the published Chinese pangolin reference genome assembled at the chromosome level. Use picard to remove redundant sequences generated by PCR during sequencing. Use Genome Analysis Toolkit v4.0 to obtain the SNP dataset. Use the HaplotypeCaller module to identify the SNPs of each individual corresponding to each sample. Merge the SNPs of each individual into a file to form a vcf file. Use vcftools in the VariantFiltration module to remove SNPs with minor allele frequencies less than 0.

05. Remove SNPs at the population level where the deletion site at the same position covers more than 20% of the samples to ensure that the number of alleles in the population is 2. Remove SNPs with sequencing depths distributed at all sites that are less than 2.5% or greater than 97.5%. Complete the filtering of the obtained SNPs, remove false positive sites and multiple allele sites, and obtain the reference sample SNPs dataset. Step S3: Based on KING analysis, perform phylogenetic analysis on each sample in the reference sample SNPs dataset to obtain multiple sets of the same individual groups and individual datasets. Only retain the information of individuals with specific geographical origins or good sequencing quality in the same individual groups. After removing redundant data in the same individual groups, a clean dataset is obtained. Filter the clean dataset to obtain a high-quality SNPs dataset. Step S4: Use the smartPCA program in the Eigensoft v.4.2 software package to perform PCA analysis on the high-quality SNPs dataset to obtain the principal component values, the variance values ​​corresponding to each principal component value, and the contribution values ​​of each SNP to the principal components. Through the Tracy-Widom test, determine the principal components that have a high degree of data interpretation and are true in the high-quality SNPs dataset. Sort the SNP sites in the high-quality SNPs dataset in descending order of contribution. Select SNP sites and merge the selected SNP sites as the geographic origin SNP molecular marker dataset. The selection of SNP sites is as follows: determine the proportion of sites to be selected in each principal component based on the variance ratio of the principal components. Then, multiply the total number of sites to be selected by the corresponding proportion to obtain how many sites to select in each principal component. Select the corresponding number of top-contributing SNP sites under each principal component based on the number of sites to be selected. Step S5: Take several small-scale locus combinations of gradients, and input the genotype information of these locus combinations in the reference samples and the latitude and longitude coordinates of the geographical origin of each reference sample into the geo-origination software OriGen. Use leave-one-out cross-validation to evaluate the error distance of each group of loci for geo-origination. That is, in each round, the latitude and longitude coordinates of the geographical origin of a reference sample are set to unknown, and the predicted coordinates are obtained by the geo-origination software OriGen. The distance between the predicted coordinates and the true coordinates is calculated using the R package geodist, which is the error distance of geo-origination. The geo-origination error distance of all samples in each group is counted. The average geo-origination error distance of each group and the number of individuals with an error distance greater than 500km are compared. The group of loci in the geo-origination SNP molecular marker dataset corresponding to at least one of the above two indicators being the smallest is taken as the best SNP molecular marker set. Step S6: Use Primer3 to generate candidate primer pairs for the target sites, and then use MFEprimer-3.0 as the core quality control module. After multiple rounds of iterative screening and global checks, a pool of multiplex PCR primers that do not interfere with each other in the same reaction system and can be amplified efficiently and specifically is obtained. Step S7: Extract total DNA from the sample of unknown origin to be identified, use a customized multiplex PCR primer set to target and capture the sample and construct a second-generation sequencing library, use a high-throughput second-generation sequencing platform to sequence the library to obtain a second-generation sequencing dataset, and perform genotyping data processing on the second-generation sequencing dataset to obtain the genotype information of the optimal SNP molecular marker set for the sample to be tested. Step S8: Use the R package OriGen for source tracing analysis. Input the genotype information of the best SNP molecular marker set obtained in S7 into the OriGen software to obtain the predicted longitude and predicted latitude information of the geographical location of the sample to be identified, as well as the predicted posterior probability of each sample, as the prediction result output.

2. The geographical origin tracing method for the Chinese pangolin based on multiplex PCR primer sets according to claim 1, characterized in that, Reference samples' whole-genome sequencing data and geographic coordinates of origin are sourced from published literature and newly collected samples.

3. The geographical origin tracing method for the Chinese pangolin based on multiplex PCR primer sets according to claim 1, characterized in that, The same group of individuals is filtered and selected using the following steps: In PLINK v.2.0, the "make-king" command is used to iterate through and calculate the KING-robust coefficients between any two individuals in the reference sample SNPs dataset. Any two individuals with KING-robust coefficients greater than 0.354 are constructed into a group of identical individuals. Within each group of identical individuals, the individual with the highest sequencing coverage is selected for subsequent analysis, and the remaining identical individuals are filtered out from the dataset.

4. The geographical origin tracing method for the Chinese pangolin based on multiplex PCR primer sets according to claim 1, characterized in that, Step S5 also includes: site combination simplification operation; The site combination simplification process includes the following steps: Step S511: Use the ancestor estimation function of the software HITSNP to construct training data based on the simulated hybrid population of purebred samples, and select the predicted probabilities output by LR-I, SVM-I, and RF-I in Classifier I to train Classifier II; Step S512: Use Classifier II to distinguish whether the samples in the best SNP molecular marker set are purebreds or hybrids and infer the number of ancestors. Then, based on the threshold that the LR-I prediction probability is greater than 0.05, lock the specific ancestral species of the hybrid sample. Step S513: Using the minimum feature SNP subset screening function, with the feature SNP as the initial set, the sequence forward selection algorithm combined with the BagNB-III classifier is used to iteratively screen according to the quantitative gradient. The classification accuracy is used as the core evaluation criterion to obtain the minimum feature SNP subset that can accurately distinguish purebred populations as the best SNP molecular marker set for use in step S6.

5. The geographical origin tracing method for the Chinese pangolin based on multiplex PCR primer sets according to claim 1, characterized in that, The process of constructing the next-generation sequencing library in step S7 includes: two rounds of PCR amplification and purification, followed by detection of the length distribution of library fragments to obtain a next-generation sequencing library that meets all quality control standards.

6. The geographical origin tracing method for the Chinese pangolin based on multiplex PCR primer sets according to claim 1, characterized in that, The data processing of the second-generation sequencing dataset in step S7 includes the following steps: Step S721: The raw sequencing data is evaluated for quality using FastQC. Trimmomatic is used to remove adapter contamination sequences, low-quality bases, and reads that are too short to obtain high-quality sequencing data. Step S722: Use Burrows-Wheeler Aligner to align high-quality sequencing data to the Chinese pangolin reference genome. Use Samtools to sort, index, and convert the alignment results. Use Genome Analysis Toolkit v4.0 to obtain SNP data. Finally, use vcftools to obtain the genotype information of the optimal SNP molecular marker set for the sample to be tested.

7. An apparatus for performing the method as described in any one of claims 1 to 6, characterized in that, include: The sample information collection module is used to acquire the whole genome sequencing data and source geographical coordinates of the reference samples, and to store the whole genome sequencing data of each reference sample as a reference sample dataset after matching the source geographical coordinates one by one. The variant site detection and filtering module is used to align the sequencing data of each reference sample in the reference sample dataset to the published Chinese pangolin reference genome assembled at the chromosome level. It uses picard to remove redundant sequences generated by PCR during sequencing, uses Genome Analysis Toolkit v4.0 to obtain the SNP dataset, uses the HaplotypeCaller module to identify the SNPs corresponding to each individual sample, merges the SNPs of each individual into a single file to form a vcf file, uses vcftools in the VariantFiltration module to remove SNPs with a minor allele frequency less than 0.05, removes SNPs where the deletion site at the same position covers more than 20% of the samples at the population level, ensures that the allele count in the population is 2, and removes SNPs with sequencing depths below 2.5% or above 97.5% of all sites. This completes the filtering of the obtained SNPs, removing false positive sites and multiple allele sites, resulting in the reference sample SNPs dataset. The kinship identification module is used to perform kinship analysis on each sample in the reference sample SNPs dataset based on KING analysis, and obtain multiple sets of the same individual groups and individual datasets. Only individuals with specific geographical origins or good sequencing coverage in the same individual groups are retained. After removing redundant data in the same individual groups, a clean dataset is obtained. After filtering the clean dataset, a high-quality SNPs dataset is obtained. The PCA module is used to perform PCA analysis on high-quality SNP datasets using the smartPCA program in the Eigensoft v.4.2 software package. It obtains the principal component values ​​for each sample, the variance values ​​corresponding to each principal component value, and the contribution values ​​of each SNP to the principal components, using the Tracy-Widom algorithm. The process involves verifying and identifying principal components with high explanatory power and true accuracy in high-quality SNP datasets. SNP sites in these datasets are then sorted in descending order of contribution. Selected SNP sites are merged to form a geographic origin SNP molecular marker dataset. Each group of sites in the geographic origin SNP molecular marker dataset is compared against the mean error distance and the number of samples with an error distance greater than 500 km in the geographic source dataset. The dataset with the smallest mean error distance and the smallest number of samples with an error distance greater than 500 km is selected as the optimal site set. SNP site selection is based on the proportion of sites to be selected in each principal component according to the principal component variance ratio. The total number of sites to be selected is multiplied by the corresponding proportion to obtain the number of sites to be selected in each principal component. The corresponding number of SNP sites are then selected based on this number. The module for obtaining the optimal SNP molecular marker set takes several small-scale locus combinations at different gradients. The genotype information of these locus combinations from the reference samples and the latitude and longitude coordinates of the geographical origin of each reference sample are input into the geo-origination software OriGen. Leave-one-out cross-validation is used to evaluate the error distance of each group of loci for geo-origination. That is, in each round, the latitude and longitude coordinates of the geographical origin of a reference sample are set to unknown. The geo-origination software OriGen predicts the coordinates, and the R package geodist is used to calculate the distance between the predicted coordinates and the true coordinates, i.e., the error distance of geo-origination. The geo-origination error distance of all samples with different locus combinations is counted. The average geo-origination error distance of each group and the number of individuals with an error distance greater than 500km are compared. The group of loci in the geo-origination SNP molecular marker dataset with at least one of the above two indicators being the smallest is taken as the optimal SNP molecular marker set. The multiplex PCR primer pool acquisition module is used to generate candidate primer pairs for target sites using Primer3. Then, MFEprimer-3.0 is used as the core quality control module. Through multiple rounds of iterative screening and global checks, multiplex PCR primer pools that do not interfere with each other in the same reaction system and can be amplified efficiently and specifically are obtained. The high-quality sequencing data acquisition and genotyping module is used to extract total DNA from samples of unknown origin to be identified, use a customized multiplex PCR primer set to target and capture quality control samples and construct a second-generation sequencing library for the unknown sample to be identified, use a high-throughput sequencing platform to sequence the second-generation sequencing library to obtain a second-generation sequencing dataset, process the second-generation sequencing dataset to obtain high-quality sequencing data, and obtain the genotype information of the optimal SNP molecular marker set for the sample to be tested through genotyping data processing. The prediction result output module is used to perform source tracing analysis using the R package OriGen. It inputs the genotype information of the best SNP molecular marker set of the sample to be identified into the OriGen software to obtain the predicted longitude and predicted latitude information of the geographical location of the sample to be identified, as well as the predicted posterior probability of each sample, as the prediction result output. Genotype information of the target SNP molecular marker set of the reference sample was obtained by using vcftools to extract the genotype information of the optimal SNP molecular marker set from the whole genome SNP dataset of the reference sample.

Images