Fine mapping of anti-sporozoite gene in carassius auratus based on linkage disequilibrium analysis

Abstract

The application discloses a method for fine positioning of anti-sporozoan genes of Carassius auratus based on linkage disequilibrium analysis, relates to the technical field of fine positioning of fish genes, and comprises the following steps: collecting whole blood samples of a Carassius auratus population, extracting genomic nucleic acid sequences, constructing a whole genome marker site set with a preset marker density gradient, dividing haplotype blocks to generate an initial haplotype data set, obtaining a population linkage disequilibrium distribution map through sliding window comparison, selecting a region with a decay rate lower than a standard value as a candidate correlation section and extracting a genotype coding sequence, matching and calculating correlation strength values with sporozoan infection survival phenotype data, reconstructing haplotype phase information of the candidate section and iteratively calculating until the threshold is met if the preset threshold is not reached. The method is suitable for genetic structure differences of genomes, optimizes correlation determination processes, and improves the accuracy of anti-sporozoan gene positioning and the reliability of section screening.

Description

Technical Field

[0001] This invention belongs to the field of fish gene fine mapping technology, specifically a method for fine mapping of the antisporidian gene in Chuzhou crucian carp based on linkage disequilibrium analysis. Background Technology

[0002] Current gene mapping studies on the resistance to sporozoites in Chuzhou crucian carp mostly employ conventional linkage disequilibrium analysis techniques. After collecting whole blood samples from Chuzhou crucian carp populations and extracting genomic nucleic acids, a set of whole-genome marker sites is constructed at a fixed or uniform density. Haplotype blocks are divided for each marker site to generate initial haplotype data. A linkage disequilibrium distribution map of the population is obtained through sliding window alignment. Regions with low decay rates are selected as candidate association segments. Genotype coding sequences are extracted and matched with post-sporozoite infection survival phenotype data. The association strength is calculated in a single step, and the association results are determined. This technique is a commonly used implementation method for mapping disease resistance genes in fish at present.

[0003] Constructing whole-genome marker sites at fixed or uniform densities cannot adapt to the differentiated characteristics of linkage disequilibrium decay rates in different regions of the genome. The distribution of marker sites does not match the actual genetic structure of the genome sufficiently, leading to biases in the selection of candidate associated regions. Haplotype phase information is only determined in the initial stage. When the association strength value does not reach the preset standard, without a targeted haplotype phase reconstruction process, relying solely on a single association calculation cannot achieve accurate gene association determination. The fine-mapping resolution of the Chuzhou crucian carp antisporidian gene is significantly limited.

[0004] To address the issue of insufficient density of genome-wide marker sites and their adaptation to genome linkage disequilibrium features, the construction logic of the marker site set was optimized. To address the lack of an iterative process for haplotype phase reconstruction when association strength values ​​are insufficient, the closed-loop execution path for linkage disequilibrium analysis and association strength calculation was improved, enabling precise and detailed localization of the antisporidian gene in Chuzhou crucian carp. Summary of the Invention

[0005] This invention aims to solve at least one of the technical problems existing in the prior art;

[0006] Therefore, this invention proposes a method for fine mapping of antisporidian genes in Chuzhou crucian carp based on linkage disequilibrium analysis, including:

[0007] Whole blood samples were collected from a group of crucian carp in Chuzhou and genomic nucleic acid sequences were extracted. A set of marker sites covering the entire genome was constructed based on a pre-set marker density gradient.

[0008] The set of marked sites is divided into haplotype blocks to generate an initial haplotype data set containing information on multiple site combinations;

[0009] The initial haplotype data set is input into the linkage disequilibrium calculation process, and the population-level linkage disequilibrium distribution map is obtained by sliding window comparison.

[0010] Regions in the linkage disequilibrium distribution map with decay rates lower than the standard value were selected as candidate association segments, and the genotype coding sequences of Chuzhou crucian carp individuals were extracted from the candidate association segments.

[0011] Phenotypic record data of survival status after infection with *Spirometra* are obtained, and the genotype coding sequence is matched and associated with the phenotypic record data to generate an association strength value;

[0012] Determine whether the correlation strength value reaches a preset correlation threshold. If it does not, reconstruct the haplotype phase information of the candidate correlation segment and recalculate the correlation strength value until the correlation strength value meets the preset correlation threshold.

[0013] Furthermore, the construction of a set of marker sites covering the entire genome based on a preset marker density gradient includes:

[0014] Three levels of marker interval distance parameters were set for low density, medium density, and high density, and virtual probes were designed for the reference genome of the Chuzhou crucian carp population, respectively.

[0015] Based on the virtual probe design results, polymorphic sites with a minimum allele frequency exceeding a specific threshold in the Chuzhou crucian carp population were screened to form a primary marker pool;

[0016] Linkage equilibrium pre-check is performed on polymorphic sites in the primary marker pool to remove redundant markers that have strong linkage disequilibrium relationships with other sites, and retain independently distributed marker sites.

[0017] The independently distributed marker sites are sorted and integrated according to their chromosome positions to generate the marker site set, and the marker site set is divided into several consecutively arranged marker window units.

[0018] Further, the set of marked sites is divided into haplotype blocks to generate an initial haplotype data set containing multi-site combination information, including:

[0019] Read all genotype call data within the labeled window unit, and count the allele combination patterns appearing in each labeled window unit;

[0020] For each of the aforementioned allele combination patterns, its frequency of occurrence in the Chuzhou crucian carp population is calculated, and rare combination patterns with a frequency of occurrence of less than one percent of the population size are filtered out.

[0021] The filtered allele combination patterns are spliced ​​together according to the physical location on the chromosome to construct a haplotype fragment spanning multiple marker window units;

[0022] Each haplotype fragment is assigned a unique fragment identifier, and the marker site number and corresponding base type contained in the haplotype fragment are recorded to form the initial haplotype data set.

[0023] Furthermore, the step of obtaining the population-level linkage disequilibrium distribution map through sliding window comparison includes:

[0024] Set a fixed-width alignment window, and move the alignment window sequentially along the genome coordinates, with each movement being half the width of the alignment window;

[0025] Within the comparison window, select a baseline marker site and calculate the linkage disequilibrium statistic between the baseline marker site and all other marker sites.

[0026] Summarize all the chain imbalance statistics calculated within the comparison window and plot the chain imbalance decay curve as a function of physical distance;

[0027] Identify the coordinate positions where the slope of the chain imbalance decay curve changes abruptly, mark the coordinate positions as block boundary points, and divide independent chain imbalance blocks based on the block boundary points;

[0028] The haplotype diversity index within each independent chain imbalance block is calculated, and blocks with a haplotype diversity index greater than a preset diversity index are included in the chain imbalance distribution map.

[0029] Furthermore, regions in the chain imbalance distribution map with decay rates lower than a standard value are selected as candidate association segments, including:

[0030] Traverse each independent chain imbalance block in the chain imbalance distribution map and extract the decay rate value corresponding to the independent chain imbalance block;

[0031] The attenuation rate value is compared with a preset standard attenuation rate value to filter out slow attenuation blocks whose attenuation rate value is less than the standard attenuation rate value.

[0032] Analyze the annotation information of the slow decay blocks on the genome and exclude the slow decay blocks located in repetitive sequence regions or gene desert regions;

[0033] Gene function enrichment analysis was performed on the remaining slow decay blocks, and the slow decay blocks containing pathway genes related to immune response or cellular defense were retained.

[0034] The slow decaying blocks after screening and enrichment analysis are defined as candidate associated segments, and the start and end positions of the candidate associated segments are recorded.

[0035] Furthermore, the genotype coding sequence of individual Chuzhou crucian carp is extracted from the candidate associated region, including:

[0036] Based on the start and end positions of the candidate associated segments, retrieve haplotype fragments containing the candidate associated segments from the initial haplotype data set;

[0037] For each individual Chuzhou crucian carp, the base information of all polymorphic sites on the corresponding haplotype fragment is extracted;

[0038] The base information is converted into a digital code form, wherein a homozygous reference base is coded as zero, a heterozygous base is coded as one, and a homozygous variant base is coded as two.

[0039] The numerical codes are arranged in chromosome order to generate a genotype coding sequence representing each individual Chuzhou crucian carp in the candidate associated region;

[0040] A genotype data matrix containing the genotype coding sequences of all the individuals of the Chuzhou crucian carp was constructed, wherein each row of the genotype data matrix corresponds to one individual and each column corresponds to one marker site.

[0041] Furthermore, the acquisition of phenotypic record data on the survival status after infection with *Spirogyrfur* includes:

[0042] The experimental group of Chuzhou crucian carp was exposed to a water environment containing spores of the spirochete, and a fixed infection duration was set.

[0043] After the infection period ended, each individual in the Chuzhou crucian carp population was examined, and the infection site and infection intensity level were recorded.

[0044] Based on the pathological examination results, individuals in which no sporozoan cysts were detected were identified as having a resistant phenotype, while individuals in which sporozoan cysts were detected were identified as having a susceptible phenotype.

[0045] Assign numerical identifier three to the resistance phenotype and numerical identifier four to the susceptibility phenotype, and generate phenotypic numerical labels for each individual.

[0046] The phenotypic numerical labels of all individuals are aggregated to form a phenotypic record data vector corresponding to the row number of the genotype data matrix.

[0047] Further, the genotype coding sequence is matched and associated with the phenotypic record data to generate an association strength value, including:

[0048] The genotype data matrix and the phenotype record data vector are aligned row by row to ensure that the genotype coding sequence of each individual corresponds to its phenotype numerical label;

[0049] The rank-based nonparametric test method was used to calculate the significance of the difference between the two phenotype groups for each marker site within the candidate associated region.

[0050] By integrating the significance of differences at all marker sites, the overall association score of the candidate associated regions is calculated.

[0051] The association score is compared with the empirical distribution obtained from random permutation simulation to calculate the corrected association strength value;

[0052] Output the correlation strength value and its corresponding chromosome interval information.

[0053] Furthermore, the step of reconstructing the haplotype phase information of the candidate association segment and recalculating the association strength value if the condition is not met includes:

[0054] When the correlation strength value is less than the preset correlation threshold, the single phase reconstruction process is initiated.

[0055] Hidden Markov models were used to perform phase inference on the genotype coding sequences within the candidate associated regions to determine the haplotype origin on each chromosome;

[0056] Update the haplotype phase labeling of the corresponding region in the initial haplotype data set according to the inference results, and generate the reconstructed haplotype data;

[0057] Based on the reconstructed haplotype data, the genotype coding sequence of the Chuzhou crucian carp individual is extracted again, and the step of matching and associating the genotype coding sequence with the phenotypic record data is repeated.

[0058] The reconstruction and association calculation process is executed cyclically, and the association strength value is updated after each cycle until the association strength value exceeds the preset association threshold or the maximum number of iterations is reached.

[0059] Furthermore, it also includes:

[0060] After reaching the preset association threshold, the current candidate association segment is locked as the antisporidian gene localization region.

[0061] Extract the location information of all recombination breakpoints within the antisporidian gene localization region and analyze the distribution density of recombination breakpoints;

[0062] Identify the sub-interval with the lowest density of recombination breakpoints, and determine the sub-interval as the minimum containment region of the candidate gene;

[0063] Output the physical coordinates of the minimum contained region and the annotation list of all known genes within the region to complete the fine localization of the antisporidian gene of Chuzhou crucian carp.

[0064] Compared with the prior art, the beneficial effects of the present invention are:

[0065] A set of marker sites covering the entire genome was constructed based on a preset marker density gradient, replacing the construction method of marker sites with fixed or uniform density. The distribution of marker sites is adapted to the linkage disequilibrium decay characteristics of different regions of the genome. The haplotype block division results are consistent with the actual linkage genetic state of the Chuzhou crucian carp population genome. The linkage disequilibrium distribution map generated by sliding window alignment is closer to the real genetic distribution characteristics of the population. The selection range of candidate associated segments has improved the fit with the genetic segments related to disease resistance in the genome. The extraction targets of genotype coding sequences are more consistent with the genetic regions associated with antisporidian trait.

[0066] When the association strength value does not reach the preset association threshold, the haplotype phase information of the candidate association segment is reconstructed and the association strength value is recalculated. The haplotype phase information is dynamically adjusted according to the association determination result. The matching association logic between the genotype coding sequence and the phenotypic record data of the survival status after infection with spores is continuously optimized. The calculation result of the association strength value gradually conforms to the real association relationship between the genome and the phenotype. The gene localization process of linkage disequilibrium analysis forms a closed-loop verification mode. The determination result of the antispore gene association segment is more in line with the actual association characteristics of the genetics and phenotype of the Chuzhou crucian carp population. Attached Figure Description

[0067] Figure 1 This is a flowchart illustrating the steps of the method for fine mapping of antisporidian genes in Chuzhou crucian carp based on linkage disequilibrium analysis as described in this invention.

[0068] Figure 2 A flowchart for obtaining linkage disequilibrium distribution maps by sliding window alignment;

[0069] Figure 3 A heatmap of the genotype coding matrix for candidate segments of Chuzhou crucian carp;

[0070] Figure 4 A comprehensive analysis chart of multiple indicators for candidate related sections of Chuzhou crucian carp;

[0071] Figure 5 The curve for iterative optimization of phase reconstruction in single-type phase. Detailed Implementation

[0072] The technical solution of the present invention will be clearly and completely described below with reference to the embodiments. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0073] See Figure 1 This invention provides a method for fine mapping of antisporidian genes in Chuzhou crucian carp based on linkage disequilibrium analysis. The overall implementation scheme of this method is as follows:

[0074] Whole blood samples were collected from a population of crucian carp in Chuzhou, and genomic nucleic acid sequences were extracted. A set of marker sites covering the entire genome was constructed based on a preset marker density gradient. The marker site set was divided into haplotype blocks to generate an initial haplotype dataset containing information on multiple site combinations. The initial haplotype dataset was input into a linkage disequilibrium calculation process, and a population-level linkage disequilibrium distribution map was obtained through sliding window alignment. Regions in the linkage disequilibrium distribution map with decay rates lower than a standard value were selected as candidate association segments, and genotype coding sequences of individual crucian carp in Chuzhou were extracted from these candidate association segments. Phenotypic records of survival status after infection with *Pseudomonas aeruginosa* were obtained, and the genotype coding sequences were matched and associated with the phenotypic records to generate association strength values. It was determined whether the association strength value reached a preset association threshold. If not, the haplotype phase information of the candidate association segments was reconstructed, and the association strength value was recalculated until the association strength value met the preset association threshold.

[0075] In one embodiment of the present invention, three levels of marker spacing parameters—low density, medium density, and high density—are set, and virtual probes are designed for the reference genome of a Chuzhou crucian carp population. Based on the virtual probe design results, polymorphic sites with a minimum allele frequency exceeding a specific threshold in the Chuzhou crucian carp population are screened to form a primary marker pool. Linkage equilibrium pre-screening is performed on the polymorphic sites in the primary marker pool, removing redundant markers with strong linkage disequilibrium relationships with other sites and retaining independently distributed marker sites. The independently distributed marker sites are sorted and integrated according to their chromosomal positions to generate a marker site set, which is then divided into several consecutively arranged marker window units.

[0076] In the specific implementation, the method constructs a set of marker sites covering the entire genome based on a preset marker density gradient. The preset marker density gradient includes three levels: low density, medium density, and high density, and virtual probes are designed for the reference genome of the Chuzhou crucian carp population. In the specific implementation, the marker interval parameter is systematically set using a formula. This formula is expressed as:

[0077]

[0078] in: Representing the Theoretical average distance between markers at each density level, parameters It is the baseline physical distance constant. These are integer level values ​​assigned to low-density, medium-density, and high-density levels, respectively. In some embodiments, the level value for the low-density level... Set as medium density level rating Set as The level value of high-density hierarchy Set as In some embodiments, the reference physical distance constant The values ​​can be adjusted based on the genome size of the Chuzhou crucian carp and the total number of expected markers. Based on the marker interval distance parameters at different levels, regularly spaced virtual probes are designed on the reference genome sequence of the Chuzhou crucian carp population, with each virtual probe corresponding to a potential marker site.

[0079] In the specific implementation, based on the virtual probe design results, genotyping was performed on samples from the Chuzhou crucian carp population to obtain polymorphic information at the corresponding position of each virtual probe. In the specific implementation, polymorphic loci with a minimum allele frequency exceeding a specific threshold were screened, and all selected loci formed a primary marker pool. In the specific implementation, the specific threshold was set to... That is, only polymorphic sites with a minimum allele frequency greater than 5% in the Chuzhou crucian carp population are preserved. Optionally, the specific threshold can be adjusted according to the actual level of genetic diversity in the population. to Adjustments are made between them. After the primary marker pool is constructed, linkage equilibrium pre-checks are performed on polymorphic sites in the primary marker pool.

[0080] In practice, linkage equilibrium pre-detection calculates the linkage disequilibrium statistic between any two polymorphic sites in the primary marker pool. The linkage disequilibrium statistic is implemented using... The value is used for measurement. In practice, if the difference between any two polymorphic sites... If the value is greater than a preset linkage threshold, then a strong linkage disequilibrium relationship is determined between the two polymorphic sites. In specific implementation, the preset linkage threshold is set to... When pairwise polymorphic sites with strong linkage disequilibrium are identified, one polymorphic site is randomly removed from each pair as a redundant marker, leaving a set of independently distributed marker sites. Optionally, the rule for removing redundant markers can be to retain polymorphic sites with the lowest allele frequency. The purpose of removing redundant markers is to reduce collinearity among markers, ensuring that the final set of independently distributed marker sites is as evenly distributed and independent as possible across the genome.

[0081] In the specific implementation, independently distributed marker loci are sorted and integrated according to their chromosome numbers and physical location coordinates on the Chuzhou crucian carp reference genome. After sorting and integration, an ordered list of marker loci is generated, which represents the set of marker loci covering the entire genome. In the specific implementation, the marker loci set is further divided into several consecutively arranged marker window units. It can be understood that when dividing the marker window units, each marker window unit is typically set to contain a fixed number of independently distributed marker loci. In the specific implementation, the number of independently distributed marker loci contained in each marker window unit is set to... Each marker window cell allows for partial overlap of independently distributed marker sites to ensure the continuity of the analysis. Ultimately, the ordered set of marker sites and the marker window cells they form provide the basic structure for subsequent haplotype block partitioning.

[0082] In one embodiment of the present invention, all genotype call data within a marker window unit are read, and the allele combination patterns appearing in each marker window unit are statistically analyzed. For each allele combination pattern, its frequency of occurrence in the Chuzhou crucian carp population is calculated, and rare combination patterns with a frequency less than one percent of the population size are filtered out. The filtered allele combination patterns are spliced ​​according to their chromosomal physical locations to construct haplotype fragments spanning multiple marker window units. A unique fragment identifier is assigned to each haplotype fragment, and the marker site number and corresponding base type contained in the haplotype fragment are recorded to form an initial haplotype dataset. See also... Figure 2 A fixed-width alignment window is set, and the window is moved sequentially along the genomic coordinates, with each movement being half the width of the alignment window. Within the alignment window, a baseline marker locus is selected, and linkage disequilibrium statistics are calculated between the baseline marker locus and all other marker loci. The linkage disequilibrium statistics calculated across all alignment windows are summarized, and a linkage disequilibrium decay curve is plotted as a function of physical distance. The coordinate locations where the slope of the linkage disequilibrium decay curve abruptly changes are identified, and these locations are marked as block boundary points. Independent linkage disequilibrium blocks are then defined based on these boundary points. The haplotype diversity index within each independent linkage disequilibrium block is calculated, and blocks with a haplotype diversity index greater than a preset diversity threshold are included in the linkage disequilibrium distribution map.

[0083] In the specific implementation, the method reads all genotype call data within the marker window unit. This genotype call data originates from sequencing and genotyping results of samples from the Chuzhou crucian carp population. Allele combination patterns appearing within each marker window unit are statistically analyzed; each allele combination pattern represents a specific arrangement of alleles at multiple polymorphic loci within the marker window unit. Specifically, for each allele combination pattern, its frequency in the Chuzhou crucian carp population is calculated. The frequency equals the number of individuals carrying that allele combination pattern divided by the total number of individuals in the Chuzhou crucian carp population. In the specific implementation, rare allele combination patterns with a frequency below one percent of the population size are filtered out, and the remaining allele combination patterns constitute a set of high-frequency allele combination patterns. Optionally, the filtering threshold can be adjusted to 0.5% or 2% depending on the population size. The filtered allele combination patterns are then spliced ​​according to their chromosomal physical locations. High-frequency allele combination patterns in adjacent marker window units are compared and connected within individuals to construct haplotype fragments spanning multiple marker window units. In practice, each haplotype fragment is assigned a unique fragment identifier, and the encoding rule for the fragment identifier can combine chromosome number and starting physical location. In practice, the marker site numbers and corresponding base types contained in the haplotype fragments are recorded, and the information of all haplotype fragments is summarized to form the initial haplotype dataset.

[0084] In practice, the generation process of the linkage disequilibrium distribution map uses a fixed-width alignment window, which slides continuously across the genome to cover the target region. Specifically, the alignment window is moved sequentially along the genomic coordinates, with each movement being half the width of the alignment window. This overlapping sliding method ensures that any location on the genome is covered by at least one complete alignment window. Within the alignment window, a baseline marker site is selected, typically a polymorphic site located in the middle of the window. Linkage disequilibrium statistics are calculated between the baseline marker site and all other marker sites within the same alignment window. These statistics can be calculated and recorded using D' values ​​or r² values. In some embodiments, multiple linkage disequilibrium statistics are calculated within a single alignment window. All linkage disequilibrium statistics calculated within the alignment windows are then summarized, grouped and averaged according to the physical distance between the corresponding marker site pairs. Plot a linkage disequilibrium decay curve as a function of physical distance, with the horizontal axis representing physical distance and the vertical axis representing the average linkage disequilibrium statistic. Identify the coordinate locations where the slope of the linkage disequilibrium decay curve abruptly changes. These abrupt changes can be determined by calculating the first derivative of the curve and finding its local extrema. Mark these coordinate locations as block boundary points. Independent linkage disequilibrium blocks can be delineated based on these block boundary points. Each independent linkage disequilibrium block represents a segment of the genome where a high level of linkage disequilibrium is maintained between marker sites. Calculate the haplotype diversity index within each independent linkage disequilibrium block. The haplotype diversity index measures the richness of haplotypes within the block. The formula for calculating the haplotype diversity index is:

[0085]

[0086] in: Represents the haplotype diversity index. This represents the total number of different haplotype fragments observed within that independent linkage disequilibrium block. Representing the The frequency of haplotype fragments appearing in the *Carassius chuzhouensis* population within this independent linkage disequilibrium block is measured. In some embodiments, a preset diversity threshold is set to 0.5. It is understood that blocks with haplotype diversity indices greater than the preset diversity threshold are included in the linkage disequilibrium distribution map, as blocks with haplotype diversity indices above the threshold generally possess richer genetic polymorphism information. Optionally, the preset diversity threshold can be adjusted between 0.3 and 0.7 depending on the research objectives. The final linkage disequilibrium distribution map is visualized with the genome coordinates on the horizontal axis and the block-based linkage disequilibrium level and haplotype diversity index on the vertical axis.

[0087] In one embodiment of the invention, each independent linkage disequilibrium block in the linkage disequilibrium distribution map is traversed, and the decay rate value corresponding to the independent linkage disequilibrium block is extracted. The decay rate value is compared with a pre-set standard decay rate value, and slow decay blocks with decay rate values ​​less than the standard decay rate value are screened out. The annotation information of the slow decay blocks on the genome is analyzed, and slow decay blocks located in repetitive sequence regions or gene desert regions are excluded. Gene function enrichment analysis is performed on the remaining slow decay blocks, and slow decay blocks containing pathway genes related to immune response or cell defense are retained. The slow decay blocks after screening and enrichment analysis are defined as candidate associated segments, and the start and end positions of the candidate associated segments are recorded. Based on the start and end positions of the candidate associated segments, haplotype fragments containing candidate associated segments are retrieved from the initial haplotype dataset. For each Chuzhou crucian carp individual, the base information of all polymorphic sites on its corresponding haplotype fragment is extracted. The base information is converted into a digital encoding form, where homozygous reference bases are encoded as zero, heterozygous bases are encoded as one, and homozygous variant bases are encoded as two. The numerical codes are arranged according to chromosome order to generate the genotype coding sequence representing each Chuzhou crucian carp individual within the candidate associated region. A genotype data matrix containing the genotype coding sequences of all Chuzhou crucian carp individuals is constructed, with each row of the genotype data matrix corresponding to one individual and each column corresponding to one marker locus.

[0088] In practice, the implementation method traverses each independent linkage disequilibrium block in the linkage disequilibrium distribution map. Each independent linkage disequilibrium block has clearly defined genome start and end coordinates. In practice, the decay rate value corresponding to each independent linkage disequilibrium block is extracted. The decay rate value describes how quickly the linkage disequilibrium statistic within that independent linkage disequilibrium block decreases with increasing physical distance. The decay rate value is calculated. The formula is:

[0089]

[0090] in: The average chain imbalance statistic near the starting position of the block. The average chain imbalance statistic near the end position of the block. The starting physical location of the block. The physical location where the block terminates. In some embodiments, the standard decay rate value... It can be set based on the median of the genomic background distribution. In practice, the decay rate value is compared with the pre-set standard decay rate value, and slow decay blocks with a decay rate value less than the standard decay rate value are screened out. Slow decay blocks mean that the linkage disequilibrium state inside them changes slowly in terms of physical distance.

[0091] In the specific implementation, the annotation information of slow decay blocks on the genome is analyzed. The annotation information comes from the repetitive sequence and gene annotation files of the Chuzhou crucian carp reference genome. In the specific implementation, slow decay blocks located in repetitive sequence regions or gene desert regions are excluded, as these slow decay blocks generally do not have the potential to encode functional genes. Optionally, repetitive sequence regions are determined by comparing with known transposon sequence databases. Gene function enrichment analysis is performed on the remaining slow decay blocks. The gene function enrichment analysis takes all annotated genes within the slow decay blocks as input and performs a hypergeometric test on a gene database of pathways related to immune response or cell defense. In some embodiments, the gene database of pathways related to immune response or cell defense includes relevant pathway entries in the KEGG or GO database. Slow decay blocks containing pathway genes related to immune response or cell defense are retained, while slow decay blocks that do not contain any known functional genes or contain genes but are unrelated to the preset pathway are excluded. It can be understood that the slow decay blocks after screening and enrichment analysis are defined as candidate associated segments. In practice, the start and end positions of candidate associated regions are recorded, and these positions are stored in the form of physical coordinates of the reference genome (e.g., base pair positions). Detailed information about the candidate associated regions is provided in Table 1.

[0092] Table 1: Example Information Table of Candidate Related Segments

[0093]

[0094] In the specific implementation, haplotype fragments containing candidate associated regions are retrieved from the initial haplotype dataset based on the start and end positions of the candidate associated regions. The initial haplotype dataset stores the chromosomal coordinate range of each haplotype fragment. For each individual Chuzhou crucian carp, base information for all polymorphic sites on its corresponding haplotype fragment is extracted. The extracted base information includes four types: adenine (A), thymine (T), cytosine (C), and guanine (G). In the specific implementation, the base information is converted into a digital encoding form, and the digital encoding conversion rules are uniform and preset. In the specific implementation, homozygous reference bases are encoded as zero, heterozygous bases as one, and homozygous variant bases as two. Optionally, the reference base refers to the base type published in the Chuzhou crucian carp reference genome at that site. The digital codes are arranged according to chromosomal order, and the arrangement order is based on the physical coordinates of the polymorphic site on the reference genome from smallest to largest. Generate the genotype coding sequence representing each individual Chuzhou crucian carp within the candidate associated region. An individual's genotype coding sequence is a string consisting of the numbers 0, 1, and 2. Construct a genotype data matrix containing the genotype coding sequences of all Chuzhou crucian carp individuals. The genotype data matrix is ​​a two-dimensional array. In the specific implementation, each row of the genotype data matrix corresponds to one individual, and each column corresponds to one marker locus.

[0095] See Figure 3 This is a heatmap of the genotype coding matrix for candidate segments of crucian carp from Chuzhou. Genotype distribution exhibits population polymorphism, with 0 / 1 / 2 differences among individuals at the same marker locus, consistent with natural population genetic characteristics. High-frequency homozygous variations are observed at some loci (concentrated in yellow), possibly related to selection pressure for sporozoite resistance phenotypes. Overall, there are no obvious complete deletions / fixations in individuals or loci, indicating good preservation of genetic diversity in the candidate segments. This matrix visually displays genotype diversity within candidate segments, helping to determine if selection sweeps or population stratification exist. This matrix is ​​the core input data for subsequent association analysis with sporozoite resistance / susceptibility phenotypes. Haplotypes can be further reconstructed based on this coding matrix to finely locate sporozoite-resistant functional genes. It also allows for rapid identification of genotype deletions, misclassifications, or abnormal individuals.

[0096] In one embodiment of the invention, a population of Chuzhou crucian carp was exposed to a water environment containing sporozoan spores for a fixed duration of infection. After the infection period ended, each individual in the Chuzhou crucian carp population was examined, and the infection site and infection intensity level were recorded. Based on the pathological examination results, individuals without detected sporozoan cysts were identified as having a resistant phenotype, and individuals with detected sporozoan cysts were identified as having a susceptible phenotype. A numerical identifier three was assigned to the resistant phenotype, and a numerical identifier four was assigned to the susceptible phenotype, generating a phenotypic numerical label for each individual. The phenotypic numerical labels of all individuals were summarized to form a phenotypic record data vector corresponding to the row number of the genotype data matrix. The genotype data matrix and the phenotypic record data vector were aligned row by row to ensure that the genotype coding sequence of each individual corresponds to its phenotypic numerical label. A rank-based nonparametric test method was used to calculate the significance of the difference between the two phenotypic groups for each marker locus within the candidate association segment. The significance results of the differences of all marker loci were integrated to calculate the overall association score of the candidate association segment. The association score is compared with the empirical distribution obtained from random permutation simulation to calculate the corrected association strength value. The association strength value and its corresponding chromosome interval information are output.

[0097] In the implementation, the experimental group of Chuzhou crucian carp was exposed to a water environment containing sporozoan spores. The concentration of sporozoan spores was standardized, and a fixed infection duration was set to ensure that some individuals within the group could produce distinguishable infection results. After the infection duration ended, each individual in the Chuzhou crucian carp group was examined under a dissecting microscope, and the infection site and infection intensity level were recorded. The infection intensity level was classified based on the number of sporozoan cysts observed. Based on the pathological examination results, individuals without detected sporozoan cysts were identified as having a resistant phenotype, and individuals with detected sporozoan cysts were identified as having a susceptible phenotype. In the implementation, a numerical identifier of three was assigned to the resistant phenotype, and a numerical identifier of four was assigned to the susceptible phenotype, generating a phenotypic numerical label for each individual. The phenotypic numerical labels of all individuals were summarized to form a phenotypic record data vector corresponding to the row number of the genotype data matrix. An example of a phenotypic record data vector is shown in Table 2.

[0098] Table 2: Example of individual phenotypic records of Chuzhou crucian carp

[0099]

[0100] In practice, the genotype data matrix and phenotype record data vector are aligned row-wise to ensure that each individual's genotype coding sequence corresponds to its phenotype numerical label. This alignment is performed based on the unique identifier of the individual's ID. In practice, a rank-based nonparametric test is used to calculate the significance of each marker locus within the candidate association segment between the two phenotype groups: a resistance phenotype group and a susceptible phenotype group. In some embodiments, the calculated significance of each marker locus is expressed as a p-value. The significance of all marker loci is then integrated; this integration can be achieved by combining the p-values ​​of all marker loci using statistical methods. Finally, the overall association score for the candidate association segment is calculated. The formula is:

[0101]

[0102] in: This represents the total number of marker sites contained within the candidate associated region. Representing the The p-values ​​for statistical significance were calculated from each marker locus. This can be understood as association score values. This is a statistic that comprehensively reflects the strength of the association between a region and a phenotype. The association score is compared with an empirical distribution obtained from a random permutation simulation, which is constructed by repeatedly randomly shuffling the correspondence between the phenotype record data vector and the genotype data matrix and recalculating the association score. The corrected association strength value is then calculated; this corrected association strength value can be an empirical p-value based on the empirical distribution.

[0103] In some embodiments, the correlation strength value Through the formula:

[0104]

[0105] in: This represents the empirical p-value obtained through random permutation simulation. Optionally, the correction process can also employ a multiple test correction method. The output includes the association strength value and its corresponding chromosomal interval information, including the chromosome number of the candidate associated region and its start and end physical locations. It can be understood that the association strength value is used to quantify the strength of statistical association evidence between a specific genomic region and the antisporidian phenotype.

[0106] See Figure 4This is a multi-indicator comprehensive analysis chart of candidate association segments for crucian carp from Chuzhou, used to comprehensively evaluate the association potential of antisporidian genes in different chromosomal candidate segments. The association score represents the overall association strength between the candidate segment and the sporidian resistance phenotype; a higher score indicates a stronger co-segregation signal between genetic variation and the resistance phenotype in that region. -log10 (corrected p-value) represents the statistical significance of the association results; a higher value indicates a lower probability of false positives and more reliable results. Gene density represents the distribution density of functional genes within the candidate segment, reflecting the coding potential and biological functional richness of the region. Integrating association strength, statistical significance, and functional gene density avoids bias caused by a single indicator, allowing for more precise identification of functional gene intervals. The chart visually displays the comprehensive potential of each candidate segment, providing a clear priority for subsequent fine-tuning.

[0107] In one embodiment of the present invention, when the association strength value is less than a preset association threshold, a haplotype phase reconstruction process is initiated. A Hidden Markov Model is used to perform phase inference on the genotype coding sequences within the candidate association segment to determine the haplotype origin on each chromosome. Based on the inference results, the haplotype phase labels of the corresponding regions in the initial haplotype data set are updated to generate reconstructed haplotype data. Based on the reconstructed haplotype data, the genotype coding sequences of individual Chuzhou crucian carp are re-extracted, and the step of matching and associating the genotype coding sequences with phenotypic record data is repeated. The reconstruction and association calculation process is executed iteratively, and the association strength value is updated after each iteration until the association strength value exceeds the preset association threshold or the maximum number of iterations is reached. After reaching the preset association threshold, the current candidate association segment is locked as the antisporidian gene localization region. The location information of all recombination breakpoints within the antisporidian gene localization region is extracted, and the recombination breakpoint distribution density is analyzed. The sub-interval with the lowest recombination breakpoint distribution density is identified, and this sub-interval is determined as the minimum containment region of the candidate gene. Output the physical coordinates of the smallest contained region and a list of annotations for all known genes within the region, thus completing the fine localization of the antisporidian gene in Chuzhou crucian carp.

[0108] In practice, when the association strength value is less than a preset association threshold, the haplotype phase reconstruction process is initiated. In this process, the preset association threshold is a predefined statistical significance standard used to determine the reliability of the genotype-phenotype association. A Hidden Markov Model (HMM) is used to infer the phase of the genotype coding sequences within the candidate association segments. The HMM utilizes the allele transmission patterns and linkage models in the population to estimate the most probable haplotype combinations and determine the haplotype origin on each chromosome, where haplotype origin refers to the allele's affiliation on the parent chromosome. In some embodiments, the parameters of the HMM are trained using an expectation-maximization algorithm. Based on the inference results, the haplotype phase labels of the corresponding regions in the initial haplotype data set are updated, generating reconstructed haplotype data. The reconstructed haplotype data provides more accurate allele linkage phase information. Based on the reconstructed haplotype data, the genotype coding sequences of individual Chuzhou crucian carp are re-extracted, converting the updated haplotype information into a digital encoding form. The step of matching and associating the genotype coding sequences with the phenotypic record data is repeated, and the association test is performed again to generate new association strength values. The reconstruction and association calculation process is executed cyclically. In each iteration, the newly generated reconstructed haplotype data is used for association analysis, and the association strength value is updated after each iteration until the association strength value exceeds a preset association threshold or the maximum number of iterations is reached. The maximum number of iterations is used to prevent infinite loops. Optionally, the maximum number of iterations can be set to 10.

[0109] In practice, after reaching a preset association threshold, the current candidate association segment is identified as the antisporidian gene localization region. This region is a genomic interval determined after multiple haplotype phase reconstructions and association calculations. The location information of all recombination breakpoints within the antisporidian gene localization region is extracted. This location information is identified by analyzing the boundaries between haplotype blocks in the population. The recombination breakpoint distribution density is analyzed, reflecting the hotspot intensity of historical recombination events per unit genome length. In practice, the recombination breakpoint distribution density... Through the formula:

[0110]

[0111] in: This represents the total number of recombination breakpoints detected within the antisporidian gene localization region. The physical span length (in megabases) represents the region where the antisporidian gene is located. The sub-region with the lowest recombination breakpoint density is identified. This sub-region indicates fewer recombination events in the genome's genetic history and a more complete haplotype structure. It can be understood that the sub-region with the lowest recombination breakpoint density is determined as the minimum containment region of the candidate gene. The minimum containment region is a more refined, presumed region containing the target gene within the antisporidian gene location region. In some embodiments, the minimum containment region is determined by sliding a fixed-width window within the antisporidian gene location region, calculating the recombination breakpoint density within each window, and selecting the window with the lowest density value as the minimum containment region. Optionally, the window width can be adjusted according to the size of the initially located region. The physical coordinates of the minimum containment region and a list of annotations for all known genes within the region are output. The physical coordinates include chromosome number, start and end positions, and the annotation list is derived from the bioinformatics annotation of the Chuzhou crucian carp reference genome. The fine mapping of the Chuzhou crucian carp antisporidian gene is completed, resulting in a clearly defined genomic coordinate range and a list of genes it contains.

[0112] See Figure 5 This is a haplotype phase reconstruction iterative optimization curve, which intuitively demonstrates the efficiency and effectiveness of the optimization process. The association strength value represents the association strength between the candidate segment and the antisporidian phenotype, which gradually increases with each iteration. The preset association threshold is the qualified threshold for association strength set by the project (value is 10). Reaching this threshold locks the candidate segment. The green solid line (square) represents the percentage increase in association strength after each iteration relative to the previous step, reflecting the optimization efficiency. The improvement rate reaches its peak (18.1%) in the initial iteration (step 1→2), and then gradually decreases, which conforms to the law of "diminishing marginal returns", indicating that the early haplotype phase reconstruction has the most significant effect on improving the association strength. When iterating to step 4, the association strength first exceeds the preset threshold (10), meeting the project termination condition, proving that the optimization process is efficient and feasible. The association strength continues to rise steadily with iteration, verifying that haplotype phase reconstruction can effectively correct genotype phase errors and enhance phenotype-genotype association signals.

[0113] The above embodiments are only used to illustrate the technical methods of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical methods of the present invention without departing from the spirit and scope of the technical methods of the present invention.

Claims

1. A method for fine mapping of antisporidian genes in Chuzhou crucian carp based on linkage disequilibrium analysis, characterized in that, The method includes: Whole blood samples were collected from a group of crucian carp in Chuzhou and genomic nucleic acid sequences were extracted. A set of marker sites covering the entire genome was constructed based on a pre-set marker density gradient. The set of marked sites is divided into haplotype blocks to generate an initial haplotype data set containing information on multiple site combinations; The initial haplotype data set is input into the linkage disequilibrium calculation process, and the population-level linkage disequilibrium distribution map is obtained by sliding window comparison. Regions in the linkage disequilibrium distribution map with decay rates lower than the standard value were selected as candidate association segments, and the genotype coding sequences of Chuzhou crucian carp individuals were extracted from the candidate association segments. Phenotypic record data of survival status after infection with *Spirometra* are obtained, and the genotype coding sequence is matched and associated with the phenotypic record data to generate an association strength value; Determine whether the correlation strength value reaches a preset correlation threshold. If it does not, reconstruct the haplotype phase information of the candidate correlation segment and recalculate the correlation strength value until the correlation strength value meets the preset correlation threshold.

2. The method for fine mapping of the antisporidian gene in Chuzhou crucian carp based on linkage disequilibrium analysis according to claim 1, characterized in that, The construction of a set of marker sites covering the entire genome based on a preset marker density gradient includes: Three levels of marker interval distance parameters were set for low density, medium density, and high density, and virtual probes were designed for the reference genome of the Chuzhou crucian carp population, respectively. Based on the virtual probe design results, polymorphic sites with a minimum allele frequency exceeding a specific threshold in the Chuzhou crucian carp population were screened to form a primary marker pool; Linkage equilibrium pre-check is performed on polymorphic sites in the primary marker pool to remove redundant markers that have strong linkage disequilibrium relationships with other sites, and retain independently distributed marker sites. The independently distributed marker sites are sorted and integrated according to their chromosome positions to generate the marker site set, and the marker site set is divided into several consecutively arranged marker window units.

3. The method for fine mapping of the antisporidian gene in Chuzhou crucian carp based on linkage disequilibrium analysis according to claim 2, characterized in that, The set of marked loci is divided into haplotype blocks to generate an initial haplotype data set containing multi-locus combination information, including: Read all genotype call data within the labeled window unit, and count the allele combination patterns appearing in each labeled window unit; For each of the aforementioned allele combination patterns, its frequency of occurrence in the Chuzhou crucian carp population is calculated, and rare combination patterns with a frequency of occurrence of less than one percent of the population size are filtered out. The filtered allele combination patterns are spliced ​​together according to the physical location on the chromosome to construct a haplotype fragment spanning multiple marker window units; Each haplotype fragment is assigned a unique fragment identifier, and the marker site number and corresponding base type contained in the haplotype fragment are recorded to form the initial haplotype data set.

4. The method for fine mapping of the antisporidian gene in Chuzhou crucian carp based on linkage disequilibrium analysis according to claim 3, characterized in that, The method of obtaining the population-level linkage disequilibrium distribution map through sliding window comparison includes: Set a fixed-width alignment window, and move the alignment window sequentially along the genome coordinates, with each movement being half the width of the alignment window; Within the comparison window, select a baseline marker site and calculate the linkage disequilibrium statistic between the baseline marker site and all other marker sites. Summarize all the chain imbalance statistics calculated within the comparison window and plot the chain imbalance decay curve as a function of physical distance; Identify the coordinate positions where the slope of the chain imbalance decay curve changes abruptly, mark the coordinate positions as block boundary points, and divide independent chain imbalance blocks based on the block boundary points; The haplotype diversity index within each independent chain imbalance block is calculated, and blocks with a haplotype diversity index greater than a preset diversity index are included in the chain imbalance distribution map.

5. The method for fine mapping of the antisporidian gene in Chuzhou crucian carp based on linkage disequilibrium analysis according to claim 4, characterized in that, Regions in the chain imbalance distribution map with decay rates lower than a standard value are selected as candidate association segments, including: Traverse each independent chain imbalance block in the chain imbalance distribution map and extract the decay rate value corresponding to the independent chain imbalance block; The attenuation rate value is compared with a preset standard attenuation rate value to filter out slow attenuation blocks whose attenuation rate value is less than the standard attenuation rate value. Analyze the annotation information of the slow decay blocks on the genome and exclude the slow decay blocks located in repetitive sequence regions or gene desert regions; Gene function enrichment analysis was performed on the remaining slow decay blocks, and the slow decay blocks containing pathway genes related to immune response or cellular defense were retained. The slow decaying blocks after screening and enrichment analysis are defined as candidate associated segments, and the start and end positions of the candidate associated segments are recorded.

6. The method for fine mapping of the antisporidian gene in Chuzhou crucian carp based on linkage disequilibrium analysis according to claim 5, characterized in that, And extract the genotype coding sequence of individual Chuzhou crucian carp from the candidate associated regions, including: Based on the start and end positions of the candidate associated segments, retrieve haplotype fragments containing the candidate associated segments from the initial haplotype data set; For each individual Chuzhou crucian carp, the base information of all polymorphic sites on the corresponding haplotype fragment is extracted; The base information is converted into a digital code form, wherein a homozygous reference base is coded as zero, a heterozygous base is coded as one, and a homozygous variant base is coded as two. The numerical codes are arranged in chromosome order to generate a genotype coding sequence representing each individual Chuzhou crucian carp in the candidate associated region; A genotype data matrix containing the genotype coding sequences of all the individuals of the Chuzhou crucian carp was constructed, wherein each row of the genotype data matrix corresponds to one individual and each column corresponds to one marker site.

7. The method for fine mapping of the antisporidian gene in Chuzhou crucian carp based on linkage disequilibrium analysis according to claim 6, characterized in that, The acquisition of phenotypic data on the survival status after infection with *Spirogyrfur* includes: The experimental group of Chuzhou crucian carp was exposed to a water environment containing spores of the spirochete, and a fixed infection duration was set. After the infection period ended, each individual in the Chuzhou crucian carp population was examined, and the infection site and infection intensity level were recorded. Based on the pathological examination results, individuals in which no sporozoan cysts were detected were identified as having a resistant phenotype, while individuals in which sporozoan cysts were detected were identified as having a susceptible phenotype. Assign numerical identifier three to the resistance phenotype and numerical identifier four to the susceptibility phenotype, and generate phenotypic numerical labels for each individual. The phenotypic numerical labels of all individuals are aggregated to form a phenotypic record data vector corresponding to the row number of the genotype data matrix.

8. The method for fine mapping of the antisporidian gene in Chuzhou crucian carp based on linkage disequilibrium analysis according to claim 7, characterized in that, Matching and associating the genotype coding sequence with the phenotypic record data to generate an association strength value includes: The genotype data matrix and the phenotype record data vector are aligned row by row to ensure that the genotype coding sequence of each individual corresponds to its phenotype numerical label; The rank-based nonparametric test method was used to calculate the significance of the difference between the two phenotype groups for each marker site within the candidate associated region. By integrating the significance of differences at all marker sites, the overall association score of the candidate associated regions is calculated. The association score is compared with the empirical distribution obtained from random permutation simulation to calculate the corrected association strength value; Output the correlation strength value and its corresponding chromosome interval information.

9. The method for fine mapping of the antisporidian gene in Chuzhou crucian carp based on linkage disequilibrium analysis according to claim 8, characterized in that, If the threshold is not reached, the haplotype phase information of the candidate association segment is reconstructed and the association strength value is recalculated, including: When the correlation strength value is less than the preset correlation threshold, the single phase reconstruction process is initiated. Hidden Markov models were used to perform phase inference on the genotype coding sequences within the candidate associated regions to determine the haplotype origin on each chromosome; Update the haplotype phase labeling of the corresponding region in the initial haplotype data set according to the inference results, and generate the reconstructed haplotype data; Based on the reconstructed haplotype data, the genotype coding sequence of the Chuzhou crucian carp individual is extracted again, and the step of matching and associating the genotype coding sequence with the phenotypic record data is repeated. The reconstruction and association calculation process is executed cyclically, and the association strength value is updated after each cycle until the association strength value exceeds the preset association threshold or the maximum number of iterations is reached.

10. The method for fine mapping of the antisporidian gene in Chuzhou crucian carp based on linkage disequilibrium analysis according to claim 9, characterized in that, Also includes: After reaching the preset association threshold, the current candidate association segment is locked as the antisporidian gene localization region. Extract the location information of all recombination breakpoints within the antisporidian gene localization region and analyze the distribution density of recombination breakpoints; Identify the sub-interval with the lowest density of recombination breakpoints, and determine the sub-interval as the minimum containment region of the candidate gene; Output the physical coordinates of the minimum contained region and the annotation list of all known genes within the region to complete the fine localization of the antisporidian gene of Chuzhou crucian carp.

Images