Molecular marker for milk goat breeding and application

Abstract

The invention discloses a molecular marker for milk goat breeding and application, and belongs to the technical field of agricultural biology. The nucleotide sequence of the molecular marker for detecting milk goat breeding is as shown in SEQ ID NO. 1. The invention also provides application of the molecular marker for milk goat breeding in milk goat breeding. According to the technical scheme, eQTL positioning in a larger range is carried out, and key genes and molecular markers influencing milk production traits are identified in combination with a co-positioning result.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of agricultural biotechnology, and relates to a molecular marker for breeding of dairy goats and application thereof. BACKGROUND

[0002] Since the breeding of dairy goats started late and the breeding foundation is weak, advanced breeding techniques need to be used to accelerate the breeding progress of dairy goats, such as molecular marker assisted breeding. The breeding technique is initially bred through phenotypic traits; in the 19th century, breeders bred by designing crossbreeding experiments; in modern breeding, molecular tools such as genetic engineering and genome selection are widely used in breeding programs, which greatly improves the accuracy and efficiency of animal improvement; in the 21st century, with the exponential growth of genomic data and the progress of AI algorithms, AI-driven tools can be used to analyze massive data, accurately predict genotype-phenotype correlations, identify new gene combinations, and optimize breeding with unprecedented accuracy and efficiency. Modern biological breeding techniques include embryo engineering, genetic engineering and molecular engineering, among which the molecular marker technology is obtained by using a heritable, identifiable DNA sequence, which is generally generated by mutation of several nucleotides of a gene or a nucleotide sequence, and has the advantages of accuracy and polymorphism. The molecular markers currently widely used mainly include microsatellite markers, single nucleotide polymorphisms, insertion-deletion markers and quantitative trait loci.

[0003] In order to obtain the molecular marker required for molecular marker assisted breeding, eQTL and GWAS co-localization can be combined. The whole genome association study is a technical means for identifying genetic variation sites related to specific traits or diseases. It is mainly based on the phenomenon of linkage disequilibrium (LD) between adjacent single nucleotide polymorphisms (SNP), and uses linear regression and logistic regression statistical methods to detect variation sites with significant effect and genetic force. In the past research, GWAS identified thousands of risk sites closely related to human diseases, but most of the sites are distributed in the intergenic region (non-coding region), which makes it very complex to analyze the association between these sites and diseases. Therefore, eQTL analysis, as a technology to connect disease risk site variation and gene expression, emerged as the times required, and laid a foundation for analyzing the pathogenic mechanism of the risk sites identified by GWAS.

[0004] In addition, traditional eQTL studies are usually based on bulk RNA-seq, which uses the average expression of different cell types and states to mask the differences in expression levels. Therefore, it is urgent to identify these additional eQTLs missed in tissue analysis to expand the scope of eQTLs and make more explicit inferences about pathogenic genes and variations. In order to facilitate the identification of cell type-specific and environment-related eQTLs, statistical methods such as Decon-eQTL and CSeQTL have been developed, which can be used for both bulk samples obtained by deconvolution analysis and single-cell data. These methods can provide more detailed cell type resolution and capture dynamic effects of eQTLs.

[0005] Co-localization first preprocesses the phenotype and gene expression data of different samples, and genotype data is generally represented in 0 / 1 / 2 encoding mode. Second, use association analysis methods to obtain GWAS summary statistics and eQTL data. Finally, use Bayesian theorem and other methods to calculate the colocalization posterior probability (CLPP) of SNPs to obtain candidate genes affecting specific traits. That is, a SNP is assumed to be related to a complex trait and the expression level of a specific gene. In this case, the gene may be considered as a potential candidate gene for the trait.

[0006] In summary, in order to locate a wider range of eQTLs and identify key genes and molecular markers affecting milk production traits in combination with co-localization results, single-cell sequencing and transcriptome sequencing of goat mammary gland tissue were performed, and the transcriptome data was deconvoluted based on the single-cell results for single-cell eQTL analysis. The obtained eQTLs were co-localized with GWAS sites, and potential molecular markers for goat breeding were mined in combination with the co-localization results. SUMMARY

[0007] In order to overcome the defects in the prior art, the present application provides a molecular marker for goat breeding and its application, which performs single-cell sequencing and transcriptome sequencing of goat mammary gland tissue, deconvolutes the transcriptome data based on the single-cell results for single-cell eQTL analysis, and co-localizes the obtained eQTLs with GWAS sites. Potential molecular markers for goat breeding are mined in combination with the co-localization results.

[0008] The technical scheme is as follows:

[0009] First, embodiments of the present invention provide a molecular marker for dairy goat breeding, the nucleotide sequence of which is shown in SEQ ID NO.1. The molecular marker is TATGGGTGTAGTTTGTGGCTCC.

[0010] Secondly, this invention provides an application of the molecular markers for dairy goat breeding described in this invention in dairy goat breeding.

[0011] Third, the present invention provides a method for screening molecular markers in dairy goat breeding, comprising the following steps:

[0012] First, mammary gland tissue was collected from six dairy goats for single-cell sequencing. The obtained data underwent quality control to remove environmental contaminants and diploids. SCVI was used to de-batch process the samples, and Seurat was used for standardization and principal component analysis (PCA). Principal components (PCs) were selected based on elbow plots, and UMAP plots were used to visualize the clustering results. The "FindAllMarkers" function in Seurat was used to identify differentially expressed genes (DEGs) in one cluster relative to all other clusters, and annotation was performed based on classic marker genes for cell type.

[0013] Bulk RNA-seq data from mammary gland tissues of seven self-tested and collected dairy goats were quality controlled and filtered using FASTP. Then, a goat ARS1.2 reference genome index was constructed using STAR, and clean reads were aligned to the goat reference genome. GCF_001704415.2_ARS1.2_genomic.gtf from NCBI was used as the gene annotation file, and the normalized expression level (TPM) and counts expression level of each gene were estimated.

[0014] Then, single-cell expression matrices were extracted from breast tissue, and the `buildSignatureMatrixMAST` function from the DWLS package was used to obtain gene expression matrices corresponding to different genotypes in breast tissue. Next, the proportions of different cell types in different individuals within the bulk sample were calculated using the DWLS package. For cell-type-interaction (cseQTL), analysis was performed using OmiGA (v1.0.1). The `cis_interaction` mode was used, with the corresponding genotype, phenotype, covariates (age, sex, genotype pc, phenotype pc), and the predicted proportions of cell types with a median greater than 0.1 as input. The results were filtered using `pval_g1 < 0.05`, `pval_g2 < 0.2`, and `qval_g2 < 0.2`. For cseQTL, the TPM matrix was first obtained from the RNA-seq data of the breast bulk sample, and then, based on the obtained proportions of different cell types in the bulk sample, these two factors were used as input to predict the cell type-specific expression matrix. The cell type-specific expression matrix was then used as input and analyzed in OmiGa's cis mode. Significant cseQTLs were identified by using pval_g1_acat < 0.05 as a threshold.

[0015] Next, to identify shared associated variants between eQTL and GWAS sites in breast tissue, colocalization analysis was performed using the Bayesian statistical procedure Coloc (v5.1). All sites within 1000000 bp upstream and downstream of each lead SNP in the GWAS results that overlapped with ieQTL / cseQTL were extracted. For the genes whose expression levels were affected by these sites, the posterior probability PP4 was obtained from the "coloc.abf" function using default parameters. Variants corresponding to ieGene-trait combinations with PP4 > 0.5 were defined as significant colocalization sites.

[0016] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0017] Advantages over existing technologies: Traditional eQTL localization is typically based on bulk RNA sequencing (bulk RNA-seq) of large tissue samples. The association analysis data used represents the average gene expression levels across different cell types, masking differences in expression levels across cell types and states. Therefore, there is an urgent need to identify these additional eQTLs missed in tissue analysis to broaden the eQTL scope and make more definitive inferences about pathogenic genes and variants. To facilitate the identification of cell type-specific and environment-related eQTLs, statistical methods such as Decon-eQTL and CSeQTL have been developed. These methods can be used for both bulk samples obtained through deconvolution analysis and single-cell data, providing finer cell type resolution and capturing the dynamic effects of eQTLs.

[0018] This method combines the sample cell proportions of the bulk obtained by deconvolution with the cell type-specific gene expression matrix predicted by bMIND to obtain ieqtl (cell type interaction eqtl) and cseqtl (cell type stratification eqtl). The eqtl obtained based on these methods can discover more refined biological signals. Attached Figure Description

[0019] Figure 1 This invention provides a technical roadmap for the development of molecular markers in dairy goat breeding.

[0020] Figure 2 The effect of different genotypes at this site in cseQTL on the expression level of gene RAET1E in luminal progenitor cells;

[0021] Figure 3 The effect of different genotypes at this site in cseQTL on the expression level of the LOC102177708 gene in luminal progenitor cells;

[0022] Figure 4 To determine the location and p-value distribution of the gwas site within the colocalization region, and the location and p-value distribution of the cseqtl site that influences the expression level of the LOC102177708 gene;

[0023] Figure 5 To determine the location and p-value distribution of the gwas site within the colocalization region, and the location and p-value distribution of the cseqtl site that influences RAET1E gene expression. Detailed Implementation

[0024] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0025] Unless otherwise specified, the experimental and detection methods in the following embodiments are conventional methods; the reagents and materials mentioned are commercially available unless otherwise specified; and the index data are measured using conventional methods unless otherwise specified.

[0026] First, embodiments of the present invention provide a molecular marker for dairy goat breeding, the nucleotide sequence of which is shown in SEQ ID NO.1. The molecular marker is TATGGGTGTAGTTTGTGGCTCC.

[0027] Secondly, this invention provides an application of the molecular markers for dairy goat breeding described in this invention in dairy goat breeding.

[0028] Third, this invention provides a method for screening molecular markers in dairy goat breeding as described in this invention.

[0029] First, mammary gland tissue was collected from six dairy goats for single-cell sequencing. The obtained data underwent quality control to remove environmental contamination and diploids. Batch processing was performed using SCVI (Gayoso A et al. 2022), followed by normalization and principal component analysis (PCA) using the Seurat package. Principal components (PCs) were selected based on elbow plots, and the clustering results were visualized using UMAP plots. The "FindAllMarkers" function in Seurat was used to identify differentially expressed genes (DEGs) in one cluster relative to all other clusters, and annotation was performed based on classic marker genes for cell type.

[0030] Bulk RNA-seq data from mammary gland tissues of seven self-tested and collected dairy goats were quality controlled and filtered using FASTP (Chen S2023). Then, a goat ARS1.2 reference genome index was constructed using STAR (Chen S 2023), and clean reads were aligned to the goat reference genome. The gene annotation file GCF_001704415.2_ARS1.2_genomic.gtf from NCBI was used, and the normalized expression level (TPM) and counts expression level of each gene were estimated.

[0031] Then, single-cell expression matrices were extracted from breast tissue, and the `buildSignatureMatrixMAST` function from the DWLS package was used to obtain gene expression matrices corresponding to different genotypes in breast tissue. Next, the proportions of different cell types in different individuals within the bulk sample were calculated using the DWLS package. For cell-type-interaction (cseQTL), analysis was performed using OmiGA (v1.0.1). The `cis_interaction` mode was used, with the corresponding genotype, phenotype, covariates (age, sex, genotype pc, phenotype pc), and the predicted proportions of cell types with a median greater than 0.1 as input. The results were filtered using `pval_g1 < 0.05`, `pval_g2 < 0.2`, and `qval_g2 < 0.2`. For cseQTL, the TPM matrix was first obtained from the RNA-seq data of the breast bulk sample, and then, based on the obtained proportions of different cell types in the bulk sample, these two factors were used as input to predict the cell type-specific expression matrix. The cell type-specific expression matrix was then used as input and analyzed in OmiGa's cis mode. Significant cseQTLs were identified by using pval_g1_acat < 0.05 as a threshold.

[0032] Next, to identify shared associated variants between eQTL and GWAS sites in breast tissue, colocalization analysis was performed using the Bayesian statistical procedure Coloc (v5.1). All sites within 1000000 bp upstream and downstream of each lead SNP in the GWAS results that overlapped with ieQTL / cseQTL were extracted. For the genes whose expression levels were affected by these sites, the posterior probability PP4 was obtained from the "coloc.abf" function using default parameters. Variants corresponding to ieGene-trait combinations with PP4 > 0.5 were defined as significant colocalization sites.

[0033] go through Figure 1 The described route leads to a C-to-T mutation at position 73952482 on chromosome 9. No related papers were found. The following results support the view that this site is valuable:

[0034]

[0035]

[0036] Both Table 1 and Table 2 show that... Figure 1 The different genotypes at the sites shown are significantly different in luminal progenitor cells for... Figure 2 , Figure 3The expression levels of the genes shown (LOC102177708, RAET1E) have different effects. Table 1 presents the GWAS results for the milk fat percentage trait in 974 dairy goat populations, along with co-localization probability information obtained from CSeqTL co-localization with luminal progenitor cells. Table 2 shows the p-value, minor allele frequency, effect size, and standard deviation of this locus in the GWAS results for the milk fat percentage trait in the 974 dairy goat populations. Figure 4 , Figure 5 The LD (linkage disequilibrium strength) of all sites within a 1mbp range upstream and downstream of this site in gwas and cseqtl, and the corresponding p value.

Claims

1. A molecular marker for breeding dairy goats, characterized in that, The nucleotide sequence of the molecular marker is shown as SEQ ID NO.

1.

2. The application of the molecular marker for breeding dairy goats according to claim 1 in breeding dairy goats.

3. The method of screening for a molecular marker of dairy goat breeding according to claim 1, wherein the marker is a marker of the nucleotide sequence of SEQ ID NO: 1 or 2. Comprising the following steps: Step 1, collecting the mammary gland tissues of 6 dairy goats for single cell sequencing; Step 2, using FASTP to perform quality control and filtering on the Bulk RNA-seq data of the mammary gland tissues of 7 self-tested and collected dairy goats; Step 3, extracting the single cell expression matrix of the mammary gland tissues, and using the function buildSignatureMatrixMAST in DWLS package to obtain the gene expression matrix corresponding to different genotypes in the mammary gland tissues; Step 4, performing co-localization analysis using the Bayesian statistical program of Coloc v5.

1.

4. The method of claim 3, wherein, In step 1, the obtained data is subjected to quality control to remove environmental pollution and diploids, SCVI is used to process the samples for batch removal, Seurat package is used for standardization, principal component analysis, and a certain number of principal components are selected according to the elbow diagram, UMAP diagram is used for visualization of the clustering results, and the "FindAllMarkers" function in Seurat is used to identify the differentially expressed genes of a cluster relative to all other clusters, and then the classical marker genes of cell types are used for annotation.

5. The method of claim 3, wherein, In step 2, STAR is used to build an index of the goat ARS1.2 reference genome, and the clean reads are aligned to the goat reference genome; GCF_001704415.2_ARS1.2_genomic.gtf in NCBI is used as a gene annotation file, and then the normalized expression level and Counts expression level of each gene are estimated.

6. The method of claim 3, wherein, In step 3, DWLS package is used to calculate the proportion of different cell types in different individuals of the Bulk sample, and for cell-type-interaction eQTL, OmiGA software is used for analysis, and the cis_interaction mode therein is used, the corresponding genotype, phenotype, covariate and predicted cell type proportion greater than 0.1 are input, and the obtained results are filtered by pval_g1 < 0.05, pval_g2 < 0.2, qval_g2 < 0.2, and for cseQTL, first, the TPM matrix is obtained from the RNA-seq data of the Bulk sample of the mammary gland, and then the proportion of different cell types in the Bulk sample is obtained, which are used as input to predict the cell type-specific expression matrix, and then the obtained cell type-specific expression matrix is used as input to analyze the cis mode of OmiGa, and the obtained results are filtered by pval_g1_acat < 0.05 as a threshold to obtain significant cseQTL.

7. The method of claim 6, wherein, The covariates include age, gender, genotype pc, and phenotype pc.

8. The method of claim 3, wherein, In step 4, all the sites overlapping with ieQTL / cseQTL within 500000bp upstream and downstream of each lead SNP in gwas results and with LD > 0.1 are extracted, and the posterior probability PP4 is obtained from the coloc.abf function for the genes affected by the expression of these sites using default parameters. The variation sites corresponding to the ieGene-trait combinations with PP4 > 0.5 are defined as significant colocalization sites.

Images