Goat genetic diversity and lambing gene mining method based on whole genome analysis

Abstract

The invention discloses a method for analyzing goat genetic diversity and lambing gene mining based on a whole genome, and relates to the field of genetic breeding. According to the technical scheme, the method comprises the steps of whole genome sequencing, whole genome variation data analysis, principal component analysis, population neighbor tree construction, genetic diversity analysis, phenotype data processing, variance analysis and whole genome association analysis on extracted DNA. The method can be used for identifying variation sites obviously related to the lambing number character of the Dazu black goats and breeding goats with high-reproduction characters.

Description

Technical Field

[0001] The present invention relates to the field of genetic breeding, and in particular to a method for goat genetic diversity and multi-lamb gene mining based on whole genome analysis. Background Art

[0002] The Dazu Black Goat is a superior goat breed, cultivated through long-term natural selection and artificial breeding, bred for both leather and meat. It boasts a pure black coat, stable genetics, a large build, high reproductive performance, strong disease resistance, and excellent meat quality. It is also known for its exceptional multiparity, with a lambing rate of 2.18 for first-time and 2.52 for second-time lambs, and a lamb survival rate of ≥90%. Therefore, analyzing the high fertility traits of the Dazu Black Goat and establishing a conservation and breeding system for the Dazu Black Goat are crucial. Fertility is a key factor in the development of the meat goat industry, and the number of lambs born by a ewe is an important indicator of fertility. However, existing technologies for the lambing trait in goats at the genetic level have yet to demonstrate its effectiveness. Summary of the Invention

[0003] To address the above technical issues, the present invention provides a method for whole-genome analysis of goat genetic diversity and multiplicity gene mining. This method can be used to identify variant sites significantly associated with litter size in Dazu black goats and to breed goats with high fertility traits.

[0004] The technical solution adopted by the present invention to solve its technical problems is: a method for analyzing goat genetic diversity and multi-lamb gene mining based on whole genome analysis, including whole genome sequencing of goat blood DNA, whole genome variation data analysis, principal component analysis and population neighbor tree construction, phenotypic data processing and variance analysis and whole genome association analysis.

[0005] Preferably, the whole genome variation data analysis includes raw data quality control, sequence alignment, duplication removal, chromosome splitting, chromosome merging, variation detection, variation identification, preliminary filtering and SNP site filtering quality control.

[0006] Preferably, the phenotypic data processing and variance analysis, the lambing number phenotypes sorted out include the average number of lambs, the average number of lambs per parity, the average number of lambs born in the first three parities and the average number of lambs born in the first four parities;

[0007] The fixed effects were mainly the year and season of birth of the ewes. The general linear model in SAS software was used to perform variance analysis by least squares regression to analyze the effects of fixed effects on different litter size phenotypes. The analysis model was as follows:

[0008] Y ijk =u+s i +b j +e ijk

[0009] Among them, Yijk is the phenotypic value of the growth trait, u is the population mean, s i is the gender effect, b j is the batch effect, e ijk is the random residual effect.

[0010] Advantages of the present invention:

[0011] The method can be used to identify the litter size trait of Dazu black goats and to breed goats with high fertility traits. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for describing the embodiments. Obviously, the drawings described below are only six of the drawings of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0013] Figure 1 This is a diagram showing the effect of agarose detection of genomic DNA extraction according to an embodiment of the present invention;

[0014] Figure 2 This is a population structure analysis of an embodiment of the present invention; A. Principal component analysis diagram; B. Neighborhood tree;

[0015] Figure 3 is a normal distribution quantile-quantile diagram of an embodiment of the present invention;

[0016] Figure 4 This is the correlation analysis of different lambing phenotypes according to the embodiment of the present invention;

[0017] Figure 5 The GWAS results of the lambing number trait and SNPs variation sites in the embodiment of the present invention are as follows;

[0018] Figure 6 The GWAS results of the lambing number trait and Indels variation sites in the embodiments of the present invention are as follows; DETAILED DESCRIPTION

[0019] In order to deepen the understanding of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. The embodiments are only used to explain the present invention and do not limit the scope of protection of the present invention.

[0020] Example

[0021] 1. Experimental Animals and Phenotype Sources

[0022] This example uses Dazu black goats as the example object. By screening ewes with complete lambing records from Chongqing Dazu District Tengda Animal Husbandry Co., Ltd., 73 ewes in good health were selected.

[0023] 2. Main reagents and instruments

[0024] Detailed information on the main experimental consumables, reagents, and instruments involved in this example is shown in Tables 1, 2, and 3.

[0025] Table 1. Main test consumables and manufacturers

[0026]

[0027] Table 2 Main reagents and manufacturers

[0028]

[0029]

[0030] Table 3 Main instruments and manufacturers

[0031]

[0032] 3. Sample Collection and DNA Extraction

[0033] Blood was collected from 73 ewes that were screened based on lambing record information in the early stage. After the ewes were tied up, 5 mL of neck venous blood was collected with a syringe, transferred to a 2% EDTA anticoagulant tube, gently inverted to mix, and placed in a foam box containing an ice pack. After the experiment, it was immediately transferred to a -20°C refrigerator for storage. Whole blood DNA was extracted using a DNA extraction kit (Takara, USA), the sample DNA concentration was determined by NanoDropTM One (NanoDrop Technologies, USA), and the DNA quality and integrity were detected by agarose gel electrophoresis. The results of the whole blood genomic DNA test of the 73 goats tested in this example showed that the quality was qualified, and the OD values ​​(260 / 280) were all between 1.8 and 2.0. The integrity of the DNA was detected by 1.2% agarose gel electrophoresis, and the results showed that the bands of each sample were bright and single, and the positions were consistent, such as Figure 1 As shown, it shows that the DNA extraction effect is good and the integrity is excellent.

[0034] 4. Whole-Genome Sequencing

[0035] All DNA samples were used Universal DNA Library Preparation Kit ( Universal DNA Library Prep Kit v2.0, Genomic sequencing libraries were constructed using the BGI MGISEQ-T7-PE150 platform. The average sequencing data volume for each individual was approximately 28.3 Gb.

[0036] 5. Genomic Data Analysis

[0037] The whole genome variation data in this embodiment is processed as follows:

[0038] (1) Raw data quality control: Fastp (version 0.20.1) software was used to perform quality control analysis on the sequencing data, and reads with a length of less than 40 and low-quality bases were removed to obtain high-quality reads.

[0039] (2) Sequence alignment: High-quality reads were aligned to the goat reference genome (Capra_hircus.ARS1.2) using BWA-MEM2 (version 2.2.1) software and bam files were obtained.

[0040] (3) Removal of duplicates: Use Samtools software to remove PCR duplicate sequences and sort the deduplicated bam files and build an index.

[0041] (4) Splitting chromosomes: The “HaplotypeCaller” module of the GATK (version 4.5.0.0) software was used to split the bam file of each goat genome after deduplication into single chromosome files.

[0042] (5) Merging chromosomes: The “CombineGVCFs” module of GATK (version 4.5.0.0) software was used to merge the same chromosome data files of 73 goats.

[0043] (6) Detection of mutations: The “GenotypeGVCFs” module of GATK (version 4.5.0.0) software was used to detect the mutations of each chromosome.

[0044] (7) Identification of variants: SNP variants were identified using the “SelectVariants” module of GATK (version 4.5.0.0) software.

[0045] (8) Preliminary filtering: The variant files were preliminarily filtered using the “VariantFiltration” module of the GATK (version 4.5.0.0) software with the following parameters: “QD<2.0||MQ<40.0||FS>60.0||SOR>3.0||MQRankSum<-12.5||ReadPosRankSum<-8.0”.

[0046] (9) Merge VCF files: Merge the filtered chromosome variation files into a total vcf file.

[0047] (10) SNP site filtering quality control: All SNP site VCF files obtained in the previous step were filtered using PLINK v1.90 software, with the specific parameters set as --maf 0.05 --geno 0.05 --mind 0.05 --hwe 1e-5. The filtered files were filled with missing genotype sites using Beagle v5.4 software, and 13,103,093 high-quality SNP sites were retained for subsequent analysis. Whole-genome resequencing was performed on 73 Dazu black goat ewes, and a total of 2.07 Tb of raw data was obtained. After quality control filtering, 2.02 Tb of high-quality data (High quality reads, HQRs) were obtained. The GC content distribution of the individual genomes measured ranged from 41.48% (DH43) to 43.2% (DH8). The specific sequencing data of each sample individual is shown in Table 4. A total of 32,662,996 SNPs and indels were identified by aligning high-quality reads with the goat reference genome (Capra_hircus.ARS1.2). After quality control and filtering using PLINK software, 14,224,297 high-quality SNPs and indels were retained, including 13,103,093 SNPs and 1,118,126 indels.

[0048] Table 4. Sequencing status of 73 samples

[0049]

[0050]

[0051]

[0052] 6. Principal Component Analysis and Group Proximity Tree Construction

[0053] Principal component analysis (PCA) was performed using PLINK v1.90 software with the --pca and --cluster --matrix parameters. The PCA results were visualized using the "ggplot2" R package in RStudio. The "ggtree" R package was used in RStudio to construct a .nwk file based on the IBS matrix file. The .nwk file was then imported into the iTOL online platform (https: / / itol.embl.de / personal_page.cgi) to generate a proximity tree visualization.

[0054] The genetic structure of the Dazu black goat population was analyzed based on whole genome data, using neighbor-joining (NJ) tree construction and principal component analysis (PCA) methods. Figure 2 B shows the population structure results based on the NJ tree, in which 73 ewes are clustered into 6 kinship groups. PCA diagram ( Figure 2 A) In the results, the samples were distributed in different regions, indicating that there was a genetic structure stratification phenomenon in the population.

[0055] VII. Calculation and Analysis of Common Indicators of Genetic Diversity

[0056] Population genetic diversity parameters (expected heterozygosity, observed heterozygosity, minimum allele frequency, and consanguinity coefficient) were analyzed using PLINK v1.90 software. Specific parameters included the following: --hardy for expected and observed heterozygosity; --freq for allele frequency; and --het for consanguinity coefficient. Descriptive statistical analysis was performed using the awk command-line tool to calculate basic statistics for each parameter, such as maximum, minimum, mean, and standard deviation. The results are shown in Table 5. The H values ​​of the 73 Dazu Black goats were quantified. O The mean is 0.2733, H E The mean of the observed heterozygosity is 0.3104, and the observed heterozygosity is slightly lower than the expected heterozygosity; the population inbreeding coefficient F is The mean value was 0.1194, and the mean MAF value was 0.2251, suggesting that this population may be facing a genetic bottleneck or structural differentiation, and mild inbreeding may have had a potential impact on production performance (such as decreased reproductive rate and growth rate).

[0057] Table 5. Genetic diversity statistics of Dazu black goat population

[0058]

[0059] 8. Phenotypic Data Processing and Variance Analysis

[0060] The litter size phenotypes were arranged as follows: average litter size (AL), average multiparous litter size (AML), average litter size in first 3 parities (3AL), and average litter size in first 4 parities (4AL). The fixed effects mainly considered in this embodiment were the year of birth and season of birth of the ewes. The general linear model in SAS software was used to perform variance analysis by least squares regression to analyze the effects of fixed effects on different litter size phenotypes. The analysis model was as follows:

[0061] Y ijk =u+s i +b j +e ijk

[0062] Among them, Y ijk is the phenotypic value of the growth trait, u is the population mean, s i is the gender effect, b j is the batch effect, is the batch effect, e ijk is a random residual effect. On the one hand, this experiment continuously collected and statistically analyzed the phenotypic data of the number of lambs of 73 Dazu black goats. The results are shown in Table 6. Figure 3 The normal distribution QQ plot and Figure 4 As can be seen from the correlation heat map, the phenotypes of each lamb number trait are basically in line with the normal distribution, and there is a significant positive correlation between each phenotype. On the other hand, this experiment conducted a significance analysis of the two fixed effects of birth year and birth season on the different lamb number traits of 73 Dazu black goats. The analysis results found that the birth season did not have a significant effect on these lamb number phenotypes, but the birth year effect had a significant effect on the average number of lambs born in the first three litters and the average number of lambs born in the first four litters (Table 7). The above results show that the birth year effect of the ewe has a certain influence on the individual phenotypic data of the Dazu black goats in this embodiment.

[0063] Table 6. Descriptive statistics of litter size of Dazu black goats

[0064]

[0065] Table 7. Significance analysis of different effects of litter size traits in Dazu black goats (P values)

[0066]

[0067] 9. Genome-wide Association Analysis

[0068] This example uses a multiple loci linear mixed model (MLMM) to conduct genome-wide association analysis and screen for markers associated with growth trait phenotypic data. In order to correct the influence of sex effect and batch effect on the accuracy of the results, sex and batch were included in the analysis as covariates. Based on the MLMM model, a GWAS was performed on the litter size trait and SNPs variant sites of Dazu black goats. The analysis results showed that the QQ plot showed that most points were closely distributed near the diagonal line, with good fit and reliable analysis results ( Figure 5 ). According to the significant threshold line P < 7.63×10 -8, a total of 5 variant sites were identified to be significantly associated with the AL phenotype (CHR4_53479635, CHR11_63709261, CHR15_78514978, CHR15_6114208, CHR17_8788644), and annotated to 4 candidate genes, including CHN2, ALDH2, CSTPP1, and LRP4. GWAS analysis of AML identified a total of 7 significant association sites, of which 3 sites on chr4, chr11, and chr15 were consistent with the results of AL traits, and the other 4 were: CHR2_44956710, CHR5_34556753, CHR19_46676205, and CHR24_1072130, and annotated to 5 candidate genes: BMPR2, CRY1, GRIP1, MRC2, and TLK2. GWAS analysis detected eight significant association loci (CHR1_31001515, CHR3_1385534, CHR4_86259818, CHR7_6542179, CHR10_83282334, CHR14_68545425, CHR18_14967122, and CHR21_31261773), but only three candidate genes (OTOS, AP3B1, and SCAPER) were annotated. Regarding 4AL, three significant loci (CHR4_12693823, CHR5_69101314, and CHR8_72470160) were detected, and five candidate genes (TPK1, CRY1, DOCK5, KCTD9, and GNRH1) were annotated.

[0069] Similarly, a GWAS was conducted on the litter size trait and Indels variant sites of Dazu black goats based on the MLMM model ( Figure 6 ), according to the significant threshold line P < 8.94×10 -7 A total of three variants were identified that were significantly associated with the AL phenotype (CHR1_96770792, CHR9_31439153, and CHR18_1546384), and these variants were annotated to two candidate genes, including GPR160 and PHC3. A GWAS analysis of AML identified three variants (CHR1_104646627, CHR9_31439153, and CHR18_1546384), and annotated to one candidate gene, SHOX2. A GWAS analysis of AL detected two significant association loci (CHR3_30928016 and CHR5_46725836), but only annotated to two candidate genes (OTOS and GRIP1). Regarding 4AL, two significant loci (CHR10_91957410, CHR13_1141268) were detected, and one candidate gene, AP3B1, was annotated.

[0070] 10. Candidate gene identification and gene function enrichment analysis

[0071] This embodiment adopts the Bonferroni correction method to improve the accuracy of identification and reduce the influence of false positives. Specifically, two threshold lines are set: when the P value is less than (1 / N), the site is considered to be significantly associated with the phenotype, where N is the total number of sites. ANNOVAR software is used to annotate candidate genes in the 0.5Mb region upstream and downstream of the significant site. Due to the incomplete database of goat enrichment analysis, the goat gene symbol is converted into a mouse homologous gene in the ensemble, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Otology (GO) enrichment analysis is performed on the candidate genes using the "Gene-list Enrichment" module of the KOBAS online website (http: / / bioinfo.org / kobas / genelist / ), and the corrected P value (Corrected P-Value) <0.05 is set as the significant enrichment threshold. Functional enrichment results show that the above-mentioned candidate genes are enriched in 89 KEGG pathways and 99 GO entries. Among them, 13 KEGG pathways were significantly enriched (Table 8), including GnRHsecretion and Axon guidance, etc. In addition, 185 GO terms were significantly enriched, including nucleuslocalization, intracellular signal transduction, receptor complex, and signaltransduction (Table 9).

[0072] Table 8. KEGG pathways significantly enriched in GWAS candidate genes

[0073]

[0074]

[0075] Table 9. Top 20 GO terms significantly enriched in GWAS candidate genes

[0076]

[0077] The technical features of the above embodiments can be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0078] The above-described embodiments merely represent several implementation methods of the present application. While the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the present invention. It should be noted that a person skilled in the art could make various modifications and improvements without departing from the spirit of the present application, all of which fall within the scope of protection of the present application. Therefore, the scope of protection of the present patent application shall be determined by the appended claims.

Claims

1. A method for mining goat genetic diversity and multi-lamb gene based on whole genome analysis, characterized in that: It includes whole genome sequencing of goat blood DNA, whole genome variation data analysis, principal component analysis and population neighbor tree construction and genetic diversity analysis, phenotypic data processing and variance analysis and whole genome association analysis.

2. The method for mining goat genetic diversity and multi-lamb gene based on whole genome analysis according to claim 1, wherein: The whole genome variation data analysis includes raw data quality control, sequence alignment, duplication removal, chromosome splitting, chromosome merging, variation detection, variation identification, preliminary filtering and SNP site filtering quality control.

3. The method for mining goat genetic diversity and multiple lamb genes based on whole genome analysis according to claim 1, wherein: Population genetic diversity parameters include: The hardy parameter calculates expected heterozygosity and observed heterozygosity; The freq parameter calculates the allele frequency; The het parameter calculates the inbreeding coefficient; Descriptive statistical analysis uses the awk command line tool to calculate the basic statistics of each parameter, including maximum value, minimum value, mean value and standard deviation.

4. The method for mining goat genetic diversity and multi-lamb gene based on whole genome analysis according to claim 1, wherein: The phenotypic data processing and variance analysis, the lambing phenotypes sorted out include average lambing number, average lambing number per parity, average lambing number in the first three parities and average lambing number in the first four parities; The fixed effects were mainly the year and season of birth of the ewes. The general linear model in SAS software was used to perform variance analysis by least squares regression to analyze the effects of fixed effects on different litter size phenotypes. The analysis model was as follows: Y ijk =u+s i +b j +e ijk Among them, Y ijk is the phenotypic value of the growth trait, u is the population mean, s i is the gender effect, b j is the batch effect, e ijk is the random residual effect.

Images