Queen bee breeding value evaluation and core population selection method based on multi-omics data

Abstract

The invention discloses a queen bee breeding value evaluation and core population selection method based on multi-omics data, and relates to the technical field of bee genetic breeding and multi-omics data application, and the method comprises the steps: obtaining candidate queen bee genome, intestinal microbiome and host transcriptome data; analyzing data and screening disease-resistant related core microbial markers, immune genes and bacterium-gene interaction network characteristics; by taking the characteristics as input, constructing a machine learning model to calculate a disease-resistant genome-microbiome combined breeding value; combining the breeding value, a core group genetic coefficient and an intestinal flora diversity index to construct a comprehensive selection index; according to the method, accurate evaluation of disease-resistant breeding values is achieved through multi-omics fusion and machine learning, genetic relationship and flora diversity optimization selection are combined, the disease resistance and genetic diversity of the core populations are improved, and a scientific basis is provided for queen bee breeding.

Description

Technical Field

[0001] This invention relates to the field of bee genetics and breeding and multi-omics data application technology, and more specifically, to a method for evaluating queen bee breeding value and selecting core populations based on multi-omics data. Background Technology

[0002] As the core of the bee colony's genetic information transmission, the queen bee's disease resistance directly determines the colony's survival and the economic benefits of apiculture. Diseases such as *Apis cerana* causing premature aging and decreased egg production in queen bees have led to annual losses of 15%-20% in traditional apiaries. Currently, in large-scale beekeeping, the evaluation of queen bee breeding value relies heavily on phenotypic observation, making it difficult to accurately capture recessive genetic traits such as disease resistance. This has resulted in a misconception in core population selection that prioritizes honey production over disease resistance.

[0003] Although genomic breeding technology has been applied to queen bee selection, the evaluation system based solely on SNP markers ignores the key regulatory role of the gut microbiome. Studies have confirmed that the abundance of lactic acid bacteria in the midgut of queen bees is significantly positively correlated with microsporidia resistance, while the expression of host immune genes is regulated by microbial metabolites. Single genome data cannot elucidate this gene-microbe synergistic disease resistance mechanism. Traditional transcriptome analysis is mostly limited to single-point detection after challenge and fails to combine dynamic changes in the microbiome to construct a regulatory network.

[0004] Existing selection methods suffer from three major bottlenecks: First, the calculation of breeding values ​​does not integrate microbiome information, resulting in a disease resistance prediction accuracy of less than 60%. Second, relying solely on pedigree records to control inbreeding fails to accurately quantify genomic relationships, leading to a decline in the genetic diversity of the core population. Third, the impact of gut microbiota diversity on queen health is ignored, with high-yielding lines often experiencing decreased resistance due to microbiota imbalance. Furthermore, artificial challenge experiments are costly and time-consuming, making them difficult to apply to large-scale candidate queen evaluation. Breakthroughs in multi-omics technologies offer a potential solution to these problems. Genomic SNP markers can locate disease resistance-related sites, 16S rRNA sequencing can analyze the dynamics of microbiota structure, and transcriptome data can capture differences in immune responses. The fusion of these three technologies with machine learning modeling holds promise for accurate prediction of disease resistance, providing a multi-dimensional scientific basis for core population selection.

[0005] Therefore, existing technologies suffer from problems such as lack of microbial association in disease resistance assessment, inaccurate phylogenetic quantification, and neglect of microbial diversity in selection. Summary of the Invention

[0006] To overcome the problems of existing technologies, such as lack of microbial community association in disease resistance assessment, inaccurate kinship quantification, and neglect of microbial community diversity in selection, this invention discloses a method for queen bee breeding value assessment and core population selection based on multi-omics data, which can effectively solve the above-mentioned technical problems.

[0007] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows: Methods for evaluating queen bee breeding value and selecting core populations based on multi-omics data include: Obtain multi-omics data of candidate queen bees, including genomic data, gut microbiome data, and host transcriptome data; The multi-omics data were analyzed to screen core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics related to queen bee disease resistance. Based on the core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics, a machine learning model is constructed, and the disease-resistant genome-microbiome joint breeding value of candidate queen bees is calculated through the machine learning model. A comprehensive selection index was constructed by combining the disease-resistant genome-microbiome joint breeding value, the kinship coefficient between the candidate queen bee and the existing core colony members, and the gut microbiome diversity index. Candidate queen bees are ranked according to the comprehensive selection index, and core populations are selected based on the ranking results to generate mating plans.

[0008] Preferably, the acquisition of multi-omics data of candidate queen bees includes: Before the artificial infection and challenge experiment with *Apis cerana*, midgut tissue contents were collected from candidate queen bees, and gut microbiome data were obtained by 16S rRNA sequencing. Midgut tissue was collected at specific time points after viral challenge, and host transcriptome data was obtained through transcriptome sequencing. The genomic data of candidate queen bees were obtained by whole-genome sequencing, and the genomic data included whole-genome SNP markers.

[0009] Preferably, the analysis of the multi-omics data to screen core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics related to queen bee disease resistance includes: The differences in gut microbiota structure between candidate queen bees in the disease-resistant group and the susceptible group before and after viral challenge were compared and analyzed. Multi-omics statistical methods were used to screen core probiotic genera that were significantly positively correlated with high disease resistance and harmful bacterial genera that were related to susceptibility, and the microbiome disease resistance index was calculated. Screening immune-related genes that were significantly differentially expressed in disease-resistant queen bees after challenge with the virus from transcriptome data; By using weighted gene co-expression network analysis, differentially expressed immune genes were associated with differentially abundant microbial taxa to identify key bacterial-gene regulatory pairs and construct microbial-gene interaction network features.

[0010] Preferably, the step of constructing a machine learning model based on the core microbial markers, immune-related genes, and microbial-gene interaction network characteristics, and calculating the disease-resistant genome-microbiome joint breeding value of candidate queen bees through the machine learning model, includes: Genomic SNP markers, microbiome disease resistance index before challenge, absolute abundance of key probiotics, and association strength of key bacteria-gene pairs were used as model input features. Survival time, spore load, and intestinal pathological changes of candidate queen bees after viral infection were used as disease resistance phenotypic values, and a machine learning model was constructed using gradient boosting decision tree or random forest algorithms. The model parameters were optimized through cross-validation. The optimized model was then applied to candidate queen bees that had not undergone challenge experiments to predict their potential disease resistance phenotype values, which were used as the joint breeding values ​​for disease resistance genome and microbiome.

[0011] Preferably, the formula for calculating the comprehensive selection index is: Where I represents the comprehensive selection index, The values ​​represent the combined breeding values ​​of the disease-resistant genome and microbiome. w1, w2, and w3 are the weights of each indicator and can be adjusted according to the breeding objectives.

[0012] Preferably, the comprehensive selection index is constructed by combining the disease-resistant genome-microbiome joint breeding value, the kinship coefficient between the candidate queen bee and existing core colony members, and the gut microbiome diversity index, including: Calculate the average kinship coefficient between the candidate queen bee and the existing core colony members to control the degree of inbreeding; The diversity index of the gut microbiome of candidate queen bees was calculated using the Shannon diversity index formula, and the diversity index is positively correlated with queen bee health and stress resistance.

[0013] Preferably, the step of ranking candidate queen bees according to the comprehensive selection index, selecting core populations based on the ranking results, and generating mating plans includes: All candidate queen bees were sorted in descending order based on a comprehensive selection index. The top-ranked candidate queen bees are selected to form the core population. Under the premise of ensuring high disease resistance, the mating scheme is generated by optimizing the algorithm so that the core population maintains a relatively distant kinship at the genetic level and high microbial diversity at the microbiome level.

[0014] Preferably, the multi-omics statistical method includes LEfSe analysis, which identifies microbial groups with significant abundance differences between the disease-resistant and susceptible groups, serving as core microbial biomarkers.

[0015] Preferably, the training process of the machine learning model includes: dividing the multi-omics data into a training set and a validation set, training the model using the training set, evaluating the model performance using the validation set, adjusting the model parameters based on the evaluation results, until the model prediction accuracy reaches a preset threshold.

[0016] Compared with existing technologies, the beneficial effects of this invention are as follows: Addressing the problem of insufficient microbial community association in disease resistance assessment, this invention innovatively integrates multi-omics data from the genome, gut microbiome, and host transcriptome. By comparing the differences in microbial community structure between disease-resistant and susceptible groups, core microbial biomarkers are screened. Microbial-gene interaction network features are constructed by combining differentially expressed immune genes from the transcriptome. These features, along with genomic SNP markers, are then used as input to a machine learning model to achieve accurate prediction of disease resistance phenotypes. This allows breeding value assessment to move beyond a single genome level and incorporate microbial community regulatory mechanisms, improving the comprehensiveness and accuracy of disease resistance assessment. Regarding the problem of inaccurate kinship quantification, this invention abandons the crude method of traditional pedigree recording. It accurately calculates the average kinship coefficient between candidate queen bees and existing core colony members using whole-genome data. This indicator is incorporated into the comprehensive selection index, and mating schemes are generated through an optimized algorithm. This strictly controls inbreeding at the genetic level, effectively preventing the decline of genetic diversity in the core population. Finally, addressing the problem of neglecting microbial community diversity in selection, this invention utilizes 16S... rRNA sequencing acquires gut microbiome data and calculates the Shannon diversity index, which is used as one of the core selection indicators. Combined with disease resistance breeding value and kinship coefficient, a comprehensive selection index is constructed. This ensures that the selection process not only guarantees high disease resistance but also screens out candidate queen bees with rich microbial diversity. Furthermore, by optimizing the mating program, the high diversity of the core population's microbiome is maintained, thereby improving the overall health and stress resistance of the queen bees. The above technical means work together to form a complete technical system of multi-omics fusion assessment, precise kinship quantification, and microbial diversity consideration. Ultimately, this achieves the beneficial effects of significantly improving the disease resistance of the core population, maintaining stable genetic diversity, and comprehensively enhancing health and stress resistance, providing a scientific, efficient, and sustainable solution for queen bee breeding. Attached Figure Description

[0017] To more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are merely exemplary. For those skilled in the art, other embodiments can be derived from the provided drawings without creative effort.

[0018] Figure 1 This is a diagram illustrating the steps of the method of the present invention. Detailed Implementation

[0019] The accompanying drawings are for illustrative purposes only and should not be construed as limiting the scope of this patent. To better illustrate this embodiment, some parts in the accompanying drawings may be omitted, enlarged, or reduced, and do not represent the actual product dimensions; It will be understood by those skilled in the art that certain well-known structures and their descriptions may be omitted in the accompanying drawings.

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

[0021] Example

[0022] This embodiment addresses the problems of reduced queen egg production (30%-50%) and weakened colony strength caused by Nosemaceranae disease in Oriental honeybee farming in the subtropical regions of southern my country. The aim is to cultivate a core queen bee population with high disease resistance and genetic diversity. A large-scale breeding practice was conducted, relying on a breeding base within the honeybee industry technology system. Two hundred one-day-old virgin Oriental honeybee queens from five demonstration apiaries (four standard colonies per apiary) were selected as candidate populations and randomly divided into: a model training group (120 queens, used to build a machine learning model and optimize parameters), a breeding verification group (60 queens, used to verify the accuracy of breeding value evaluation and selection effect), and a blank control group (20 queens, without challenge treatment, used as basic data reference).

[0023] Before the experiment, a complete technology platform was set up: ① A biosafety challenge laboratory (P2 level), equipped with a constant temperature and humidity system (temperature 34±0.2℃, humidity 65±3%), a sterile operating table, and an independent ventilation system to avoid cross-contamination; ② A multi-omics sequencing analysis platform, including an Illumina NovaSeq 6000 sequencer (for genome and transcriptome sequencing), a PacBio SequelII sequencer (for 16S rRNA full-length sequencing), a Nanodrop 2000 nucleic acid detector, an Agilent 2100 bioanalyzer, and other equipment; ③ A high-performance computing workstation, equipped with an Intel Xeon Gold 6330 processor (28 cores), an NVIDIA A100 GPU (80GB VRAM), and 512GB of RAM, pre-installed with professional analysis software such as QiIME 2, WGCNA, Scikit-learn, PLINK, and GATK; ④ An intelligent breeding management system, integrating modules for sample tracing, data storage, breeding value calculation, kinship analysis, and mating scheme generation.

[0024] The microsporidia strain used in the experiment was isolated and purified from the midgut of a diseased queen bee. After three generations of single-spore cloning culture, the spore concentration was calibrated to 1.5 × 10⁻⁶ using a hemocytometer. 6The samples were stored at 4°C for later use. The core experimental consumables included: DNA extraction kit (Qiagen DNeasy Blood & Tissue Kit), RNA extraction kit (ThermoFisher TRIzol Reagent), 16S rRNA amplification primers (V3-V4 region: 338F 5'-ACTCCTACGGGAGGCAGCAG-3', 806R5'-GGACTACHVGGGTWTCTAAT-3'), and real-time PCR reagent (TaKaRa SYBR Premix Ex TaqII). All consumables were sterilized.

[0025] Please see Figure 1 Methods for evaluating queen bee breeding value and selecting core populations based on multi-omics data include: Obtain multi-omics data of candidate queen bees, including genomic data, gut microbiome data, and host transcriptome data; The multi-omics data were analyzed to screen core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics related to queen bee disease resistance. Based on the core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics, a machine learning model is constructed, and the disease-resistant genome-microbiome joint breeding value of candidate queen bees is calculated through the machine learning model. A comprehensive selection index was constructed by combining the disease-resistant genome-microbiome joint breeding value, the kinship coefficient between the candidate queen bee and the existing core colony members, and the gut microbiome diversity index. Candidate queen bees are ranked according to the comprehensive selection index, and core populations are selected based on the ranking results to generate mating plans.

[0026] The acquisition of multi-omics data of candidate queen bees includes: Before the artificial infection and challenge experiment with *Apis cerana*, midgut tissue contents were collected from candidate queen bees, and gut microbiome data were obtained by 16S rRNA sequencing. Midgut tissue was collected at specific time points after viral challenge, and host transcriptome data was obtained through transcriptome sequencing. The genomic data of candidate queen bees were obtained by whole-genome sequencing, and the genomic data included whole-genome SNP markers.

[0027] The analysis of the multi-omics data, screening for core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics related to queen bee disease resistance, includes: The differences in gut microbiota structure between candidate queen bees in the disease-resistant group and the susceptible group before and after viral challenge were compared and analyzed. Multi-omics statistical methods were used to screen core probiotic genera that were significantly positively correlated with high disease resistance and harmful bacterial genera that were related to susceptibility, and the microbiome disease resistance index was calculated. Screening immune-related genes that were significantly differentially expressed in disease-resistant queen bees after challenge with the virus from transcriptome data; By using weighted gene co-expression network analysis, differentially expressed immune genes were associated with differentially abundant microbial taxa to identify key bacterial-gene regulatory pairs and construct microbial-gene interaction network features.

[0028] The process involves constructing a machine learning model based on the core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics. This model is used to calculate the disease-resistant genome-microbiome combined breeding value of candidate queen bees, including: Genomic SNP markers, microbiome disease resistance index before challenge, absolute abundance of key probiotics, and association strength of key bacteria-gene pairs were used as model input features. Survival time, spore load, and intestinal pathological changes of candidate queen bees after viral infection were used as disease resistance phenotypic values, and a machine learning model was constructed using gradient boosting decision tree or random forest algorithms. The model parameters were optimized through cross-validation. The optimized model was then applied to candidate queen bees that had not undergone challenge experiments to predict their potential disease resistance phenotype values, which were used as the joint breeding values ​​for disease resistance genome and microbiome.

[0029] The formula for calculating the comprehensive selection index is as follows: Where I represents the comprehensive selection index, The values ​​represent the combined breeding values ​​of the disease-resistant genome and microbiome. w1, w2, and w3 are the weights of each indicator and can be adjusted according to the breeding objectives.

[0030] The comprehensive selection index is constructed by combining the disease-resistant genome-microbiome joint breeding value, the kinship coefficient between the candidate queen bee and existing core colony members, and the gut microbiome diversity index, including: Calculate the average kinship coefficient between the candidate queen bee and the existing core colony members to control the degree of inbreeding; The diversity index of the gut microbiome of candidate queen bees was calculated using the Shannon diversity index formula, and the diversity index is positively correlated with queen bee health and stress resistance.

[0031] The step of ranking candidate queen bees according to the comprehensive selection index, selecting core populations based on the ranking results, and generating mating plans includes: All candidate queen bees were sorted in descending order based on a comprehensive selection index. The top-ranked candidate queen bees are selected to form the core population. Under the premise of ensuring high disease resistance, the mating scheme is generated by optimizing the algorithm so that the core population maintains a relatively distant kinship at the genetic level and high microbial diversity at the microbiome level.

[0032] The multi-omics statistical method includes LEfSe analysis, which identifies microbial groups with significant abundance differences between the disease-resistant and susceptible groups, serving as core microbial biomarkers.

[0033] The training process of the machine learning model includes: dividing multi-omics data into training set and validation set, training the model using the training set, evaluating the model performance using the validation set, adjusting the model parameters based on the evaluation results, until the model prediction accuracy reaches a preset threshold.

[0034] Gut microbiome data acquisition and sequencing Sample collection: After numbering 200 candidate queen bees, they were fasted for 24 hours before the challenge experiment. The midgut tissue was separated using aseptic dissection techniques. The midgut tissue was gently rinsed three times with sterile PBS buffer (pH 7.4) to remove surface impurities. Approximately 50 mg of midgut contents was scraped off using a sterile scalpel and immediately placed in a 2 mL centrifuge tube containing 1 mL of sterile PBS buffer. Three sterile glass beads with a diameter of 1 mm were added, and the bacterial aggregates were broken by vortexing for 1 min. The tube was then centrifuged at 5000 r / min for 5 min to collect the precipitate. 200 μL of CTAB extraction buffer was added to each tube, and the tubes were frozen at -80℃ for later extraction.

[0035] DNA Extraction and Quality Assay: Total DNA from intestinal microorganisms was extracted using the CTAB-chloroform method. The specific steps were as follows: After thawing frozen samples at room temperature, 20 μL of proteinase K (20 mg / mL) was added, and the mixture was incubated at 65°C for 30 min; an equal volume of chloroform-isoamyl alcohol (24:1) was added, and the mixture was gently inverted and centrifuged at 12000 r / min for 10 min; the supernatant was collected, and 0.8 volumes of isopropanol were added, and the mixture was incubated at -20°C for 30 min, and centrifuged at 12000 r / min for 15 min; the precipitate was washed twice with 75% ethanol, dried, and dissolved in 50 μL of TE buffer. The DNA integrity was verified by 1% agarose gel electrophoresis (clear bands without tails), the purity was detected by Nanodrop 2000 (A260 / A280 = 1.8-2.0, A260 / A230 ≥ 1.8), and the DNA concentration was quantified by Qubit 4.0 (≥ 50 ng / μL).

[0036] PCR amplification and sequencing: PCR amplification of the 16S rRNA V3-V4 region was performed using extracted DNA as a template. The reaction system (50 μL) consisted of: 25 μL of 2×Taq Plus MasterMix, 2 μL each of forward and reverse primers (10 μmol / L), 4 μL of DNA template, and 17 μL of sterile water. The reaction program was: 95℃ pre-denaturation for 5 min; 95℃ denaturation for 30 s, 55℃ annealing for 30 s, 72℃ extension for 45 s, for a total of 35 cycles; and a final extension at 72℃ for 10 min. After 1% agarose gel electrophoresis, impurities were removed using magnetic bead purification. An Illumina sequencing library (insert fragment length 400-500 bp) was constructed, and paired-end sequencing (2×300 bp) was performed using an Illumina NovaSeq 6000 sequencer. Each sample had a sequencing depth ≥60,000 reads, ensuring a bacterial community coverage ≥99%.

[0037] Host transcriptome data acquisition and sequencing Artificial challenge treatment: 120 queen bees in the model training group were artificially challenged with venom. Each queen bee was fed 50 μL of a sucrose solution containing spores of *Apis cerana* microsporidia* (spore concentration 1.5 × 10⁻⁶) via a micro-feeder. 6 (Number of bees / mL), the blank control group was fed an equal amount of sterile sucrose solution. After challenge, the queen bees were placed in an independent queen rearing frame and fed in a standard nursery colony. The rearing environment was the same as that of the original apiary.

[0038] Sample collection and RNA extraction: On day 7 post-challenge (the critical period for microsporidia infection of midgut epithelial cells), 80 queen bees were randomly selected from the model training group and 20 queen bees from the blank control group. Approximately 100 mg of midgut tissue was obtained using aseptic dissection and immediately placed in a pre-cooled centrifuge tube containing 2 mL of TRIzol reagent in liquid nitrogen. After adding sterile grinding beads, the tissue was homogenized using a tissue homogenizer (60 Hz, 2 min). After standing at room temperature for 5 min, total RNA was extracted according to the kit instructions. The specific steps were as follows: 400 μL of chloroform was added, vigorous shaking was performed for 15 s, and the tissue was stood at room temperature for 3 min; centrifuged at 12000 rpm for 15 min, and the supernatant was collected and an equal volume of isopropanol was added. The tissue was then stood at -20℃ for 20 min; centrifuged at 12000 rpm for 10 min, and the precipitate was washed twice with 75% ethanol, dried, and dissolved in 50 μL of enzyme-free water.

[0039] RNA quality testing and sequencing: RNA integrity (RIN value ≥ 8.5) was tested using an Agilent 2100 bioanalyzer, and purity (A260 / A280 = 1.9-2.1) was tested using a Nanodrop 2000. mRNA was enriched using Oligo(dT) magnetic beads, fragmented into 200-300bp fragments with fragmentation buffer, reverse transcribed to synthesize cDNA, added adapters, and then amplified by PCR to construct a transcriptome sequencing library. Paired-end sequencing (2 × 150bp) was performed using an Illumina NovaSeq 6000 sequencer, with a sequencing yield of ≥ 12Gb per sample, ensuring low-expression gene coverage of ≥ 95%.

[0040] Genome data acquisition and sequencing Sample collection and DNA extraction: Approximately 20 mg of footpad tissue was collected from 200 candidate queen bees (to avoid affecting the health of the queen bees). The tissue was placed in a centrifuge tube containing 200 μL of tissue lysis buffer, and 20 μL of proteinase K (20 mg / mL) was added. The tube was incubated overnight at 56°C. Whole-genome DNA was extracted using the Qiagen DNeasy Blood & Tissue Kit, following these steps: 200 μL of buffer AL was added, vortexed, and incubated at 70°C for 10 min; 200 μL of anhydrous ethanol was added, and vortexed; the mixture was transferred to a centrifuge column and centrifuged at 8000 rpm for 1 min; 500 μL of buffer AW1 was added, and centrifuged at 8000 rpm for 1 min; 500 μL of buffer AW2 was added, and centrifuged at 14000 rpm for 3 min; after the centrifuge column was dried, 50 μL of elution buffer was added, and the tube was incubated at room temperature for 5 min, then centrifuged at 8000 rpm for 1 min to collect the DNA.

[0041] Genome sequencing and SNP detection: After passing quality testing, the DNA was fragmented into approximately 350bp fragments using a Covaris M220 ultrasonic disruptor to construct a whole-genome resequencing library. Paired-end sequencing (2×150bp) was performed using an Illumina NovaSeq 6000 sequencer with a sequencing depth ≥30×, ensuring genome coverage ≥98%. The raw sequencing data underwent quality control (filtering reads with Q20 <80%, removing adapter contamination and low-complexity sequences) and was aligned to the Eastern Honeybee reference genome (Amel_HAv3.1) using the BWA-MEM algorithm. Repetitive sequences were removed using the Picard tool, and SNP calling was performed using GATK software (using the HaplotypeCaller module). Detected SNP sites were filtered (MAF ≥0.05, deletion rate ≤0.1, Hardy-Weinberg equilibrium P ≥1e-6) to obtain a high-quality whole-genome SNP marker set for each queen bee (approximately 2.8 million SNP sites were retained per queen bee on average).

[0042] Screening of core microbial biomarkers Disease resistance phenotype identification: On day 21 after challenge with the virus, disease resistance was assessed in all queen bees of the training group: ① Survival status was recorded; ② Spore load was measured using a hemocytometer to count the number of microsporidia spores in the midgut tissue; ③ Intestinal pathological observation was performed, including preparing paraffin sections of the midgut tissue, staining with hematoxylin and eosin (HE), and observing the integrity of intestinal epithelial cells. Based on the assessment results, the training group was divided into a disease-resistant group (45 bees, surviving with a spore load <1×10⁻⁶). 4 (Spores / mg midgut tissue, with intact intestinal epithelial cell structure), susceptible group (35 mice, dead or with spore load >5×10) 4 (The number of bees was 1 / mg of midgut tissue, and the intestinal epithelial cells were severely damaged). The remaining 40 queen bees were excluded due to ambiguous phenotypes.

[0043] Microbial community structure analysis: 16S rRNA sequencing data were analyzed using QIIME 2 software. The specific steps were as follows: quality control and noise reduction were performed using the DADA2 plugin to obtain amplicon sequence variations (ASVs); species annotation was performed based on the Greengenes database; Alpha diversity (Shannon index, Simpson index) and Beta diversity (Bray-Curtis distance) were calculated. The results showed that the Shannon index of the disease-resistant group (4.82±0.35) was significantly higher than that of the susceptible group (2.36±0.28), indicating that the gut microbiota diversity of the disease-resistant queen bee was higher.

[0044] Differential Microbial Screening: LEfSe analysis (LDA threshold = 4.0, P < 0.01) was used to screen for differentially expressed microbial groups between the resistant and susceptible groups. Results showed that the relative abundance of *Lactobacillus*, *Bifidobacterium*, and *Apibacter* was significantly higher in the resistant group than in the susceptible group (18.2% ± 2.3% vs 3.5% ± 0.8%, 12.5% ​​± 1.7% vs 2.1% ± 0.5%, 8.7% ± 1.2% vs 1.3% ± 0.3%, respectively). In the susceptible group, the relative abundance of *Enterococcus*, *Staphylococcus*, and *Klebsiella* was significantly increased (15.6% ± 2.1% vs 4.2% ± 0.9%, 10.3% ± 1.5% vs 1.3% ± 0.3%). (2.8%±0.6%, 7.8%±1.1% vs 1.5%±0.4%), the above 6 genera were identified as core microbial markers related to the resistance of queen bees to microsporidia.

[0045] Microbiome disease resistance index calculation: To quantify the disease resistance potential of gut microbiota, a microbiome disease resistance index was constructed. The calculation formula is: Disease resistance index = (relative abundance of Lactobacillus + relative abundance of Bifidobacterium + relative abundance of Bee bacteria) ÷ (relative abundance of Enterococcus + relative abundance of Staphylococcus + relative abundance of Klebsiella) × 100. The calculation results showed that the average disease resistance index of the disease-resistant group was 35.2±4.8, and that of the susceptible group was 3.9±1.2. The difference between the two groups was extremely significant (P<0.001).

[0046] Screening of immune-related genes Differentially expressed gene analysis: After quality control, the transcriptome sequencing data were aligned to the reference genome of Apis cerana using HISAT2 software, and gene expression levels (FPKM values) were calculated using HTSeq software. Differential expression analysis between the disease-resistant and susceptible groups was performed using DESeq2 software. The screening criteria were |log2FC|>2.0 and P<0.01. A total of 1426 differentially expressed genes were obtained, of which 873 genes were upregulated and 553 genes were downregulated in the disease-resistant group.

[0047] Screening of immune-related genes: Combining the Apis cerana immune gene database (integrating NCBI, Swiss-Prot database and published literature on bee immune-related genes), functional annotation was performed on differentially expressed genes, and a total of 102 differentially expressed immune-related genes were screened. Among them, defensin family genes (Defensin-1, Defensin-2), lysozyme gene, Toll signaling pathway genes (Toll, MyD88, Tube), IMD signaling pathway genes (IMD, Relish), and prophenol oxidase gene (ProPO) were significantly upregulated in the disease-resistant group, with expression levels 4.5-10.2 times that of the susceptible group; while some inflammatory response suppressor genes (such as Cactus) were significantly upregulated in the susceptible group.

[0048] Construction of microbial-gene interaction network Co-expression network construction: Weighted Gene Co-expression Network Analysis (WGCNA) was used to construct an interaction network between immune genes and core microbial biomarkers. The FPKM values ​​of 102 immune-related genes and the relative abundance of 6 core microbial biomarkers were used as input data. A soft threshold β=12 was set (satisfying the scale-free network characteristic R²>0.85). The co-expression modules were divided by dynamic tree segmentation, and 6 key modules were obtained. Among them, the blue module showed a very strong positive correlation with the abundance of Lactobacillus (r=0.82, P<0.001), and this module contained 28 immune-related genes.

[0049] Key Bacterium-Gene Pair Identification: Pearson correlation analysis was performed on 28 immune genes and 6 core microbial biomarkers in the blue module. Bacterium-gene pairs with a correlation strength r>0.75 and P<0.001 were screened, resulting in 45 key regulatory pairs. Among them, Lactobacillus and Defensin-1 had the strongest correlation (r=0.89, P<0.001), followed by Bifidobacterium and Lysozyme (r=0.85, P<0.001). Combined with GO functional annotation and KEGG pathway analysis, these key bacterial-gene pairs are mainly involved in biological processes such as innate immune response, antimicrobial peptide synthesis, and intestinal mucosal barrier repair.

[0050] Interaction network visualization: A microbial-gene interaction network was constructed using Cytoscape 3.9.1 software. Nodes represent core microbial biomarkers or immune-related genes, edges represent significant associations between the two, and the thickness of the edges indicates the strength of the association. Network analysis showed that Lactobacillus and Bifidobacterium are the core nodes in the network, regulating the expression of multiple immune genes and constituting the core regulatory network for queen bee resistance to microsporidia.

[0051] Model input and output feature determination Input feature selection: Based on the results of multi-omics analysis, the following features were selected as model inputs: ① Genomic features: 180 SNP loci that were significantly associated with disease resistance were screened from the whole genome SNP marker set (association analysis was performed using PLINK software, P<0.005). ② Microbiome characteristics: Microbiome disease resistance index, absolute abundance of Lactobacillus spp. (detected by real-time quantitative PCR, primer sequences F: 5'-AGAGTTTGATCCTGGCTCAG-3', R: 5'-CGGCTACCTTGTTACGACTT-3'), and absolute abundance of Bifidobacterium spp. before challenge; ③ Interaction network features: The association strength values ​​between Lactobacillus and Defensin-1, and between Bifidobacterium and Lysozyme were used to determine 184 input features.

[0052] Output characteristics (disease resistance phenotypic values) were determined: The weighted scoring method was used to quantify the queen bee's disease resistance into phenotypic values. The scoring indicators and weights were: survival time after challenge (40%, the longer the survival time, the higher the score), spore load (35%, the lower the load, the higher the score), and intestinal pathology score (25%, 1-5 points, the more intact the intestinal structure, the higher the score). The total score ranged from 0 to 100 points. A score ≥80 points indicated high disease resistance, 60-79 points indicated medium disease resistance, and <60 points indicated low disease resistance.

[0053] Model training and optimization Data preprocessing: The 184 input feature data of 80 queen bees in the model training group were standardized (using Z-score standardization, so that the mean of each feature is 0 and the standard deviation is 1), outliers were removed (using the 3σ criterion), and the data were divided into a training set (64 bees) and a validation set (16 bees) in an 8:2 ratio.

[0054] Model Construction: Gradient Boosting Decision Tree (GBDT) and Random Forest algorithms were used to construct machine learning models, respectively, based on the Scikit-learn library in Python. The initial parameters of the GBDT model were: learning rate 0.1, maximum tree depth 6, number of estimators 100, and minimum number of sample splits 20. The initial parameters of the Random Forest model were: number of decision trees 200, maximum tree depth 10, minimum number of sample splits 15, and random feature selection ratio 0.7.

[0055] Parameter optimization and model evaluation: Five-fold cross-validation was used to optimize the model parameters. The prediction accuracy and mean squared error (MSE) of the validation set were used as evaluation metrics. The optimal parameters of the GBDT model after optimization were: learning rate 0.08, maximum tree depth 8, number of estimators 180, and minimum number of sample splits 18. At this time, the prediction accuracy of the validation set reached 93.75%, and the MSE was 2.68. The optimal parameters of the random forest model were: number of decision trees 250, maximum tree depth 12, minimum number of sample splits 12, and random feature selection ratio 0.6. The prediction accuracy of the validation set was 87.5%, and the MSE was 4.32. After comprehensive evaluation, the GBDT model was selected as the final breeding value prediction model.

[0056] Calculation of the disease resistance genome-microbiome joint breeding value: 184 input feature data of 60 candidate queen bees in the breeding validation group that had not undergone challenge treatment were standardized and then imported into the optimized GBDT model to predict their disease resistance phenotype values. This predicted value is the disease resistance genome-microbiome joint breeding value. The results showed that the joint breeding values ​​of queen bees in the breeding validation group ranged from 62.3 to 91.7 points. Among them, 22 queen bees had breeding values ​​≥80 points (high disease resistance candidate individuals), 28 queen bees had breeding values ​​of 60-79 points (medium disease resistance), and 10 queen bees had breeding values ​​<60 points (low disease resistance). To verify the accuracy of the prediction, challenge experiments were conducted on 22 high disease resistance candidate queen bees in the breeding validation group. The results showed that 20 queen bees had actual disease resistance phenotype values ​​≥80 points, and the prediction accuracy rate reached 90.9%, indicating that the model has high reliability.

[0057] Construction of comprehensive selection index Kinship coefficient calculation: Based on whole-genome SNP data, the kinship coefficients of 60 candidate queen bees in the breeding validation group and 15 queen bees in the existing core colony of the base were calculated using PLINK software. The average value was taken as the average kinship coefficient of the candidate queen bees to control the degree of inbreeding. The calculation results showed that the average kinship coefficient ranged from 0.07 to 0.34, of which 15 candidate queen bees had an average kinship coefficient ≤0.15 (low inbreeding risk), 30 had an average kinship coefficient of 0.16-0.25 (medium inbreeding risk), and 15 had an average kinship coefficient >0.25 (high inbreeding risk).

[0058] Gut microbiome diversity index calculation: The Shannon diversity index calculation formula (H=-ΣPlnP, where P is the relative abundance of the i-th genus) was used to calculate the gut microbiome diversity index of 60 candidate queen bees in the breeding validation group. The results showed that the diversity index ranged from 2.45 to 4.98 and was significantly positively correlated with the queen bee's disease resistance (r=0.72, P<0.001), indicating that queen bees with a high diversity index have stronger resistance.

[0059] Comprehensive selection index calculation: Combining breeding objectives (prioritizing disease resistance while considering genetic diversity and gut microbiota health), the weights of each indicator are set as follows: disease resistance genome-microbiome joint breeding value (w1=0.6), average kinship coefficient (w2=0.2, the reciprocal is taken to convert to a positive indicator, i.e., the lower the kinship coefficient, the higher the converted value), and gut microbiome diversity index (w3=0.2). The comprehensive selection index calculation formula is: I=0.6×V'+0.2×(1 / R)'+0.2×H', where V' is the standardized value of joint breeding (0-100), R is the average kinship coefficient, (1 / R)' is the standardized value after conversion (0-100), and H' is the standardized value of diversity index (0-100).

[0060] Core Population Selection: Based on the comprehensive selection index, the 60 candidate queen bees in the breeding verification group were sorted in descending order, and the top 18 (I≥83 points) were selected to form a new core population. The basic characteristics of this core population are: ① Disease resistance: The average joint breeding value is 85.6 points, all of which are highly disease-resistant individuals; ② Genetic diversity: The average kinship coefficient is 0.13, which is 0.08 lower than the original core population, and the risk of inbreeding is significantly reduced; ③ Microbiome health: The average gut microbiome diversity index is 4.62, which is 0.85 higher than the ordinary population, and the microbiome structure is more stable. The core population includes 16 original highly disease-resistant candidate individuals and 2 individuals with medium breeding values ​​but low kinship coefficients and high microbiome diversity (numbered V18 and V45), achieving a balance between disease resistance, genetic diversity and microbiome health.

[0061] Mating scheme generation: To maintain the superior characteristics of the core population, the non-dominated sorting genetic algorithm (NSGA-II) is used to generate the optimal mating scheme. The objective function is to maximize the average disease resistance breeding value of the offspring, minimize the average kinship coefficient of the offspring, and maximize the predicted value of gut microbiota diversity of the offspring. Mating constraints are: ① Each queen bee is mated ≤ 3 times per year; ② Mating is prohibited when the kinship coefficient between the female and male queen bees is ≥ 0.2; ③ The female-to-male ratio of the core population is controlled at 3:1.

[0062] Ten optimal mating combinations and implementation details were generated: ①♀V03 (breeding value 91.7, kinship coefficient 0.11, diversity index 4.98) ×♂M05 (breeding value 88.5, kinship coefficient 0.09, diversity index 4.82), with an expected offspring breeding value of 90.1, kinship coefficient 0.10, and diversity index 4.90. Mating was conducted on the 6th day after the queen emerged, with an ambient temperature of 30℃ and humidity of 70%. ②♀V08 (breeding value 90.2, kinship coefficient 0.12, diversity index 4.85) ×♂M03 (breeding value 87.3, kinship coefficient 0.08, diversity index 4.76), with an expected offspring breeding value of 88.8, kinship coefficient 0.10, and diversity index 4.81. Three days prior to mating, the offspring were fed a nutritional diet containing lactic acid bacteria (concentration 1×10⁻⁶). 8③♀V12 (breeding value 89.5 points, kinship coefficient 0.10, diversity index 4.92) × ♂M07 (breeding value 86.8 points, kinship coefficient 0.10, diversity index 4.68), expected offspring breeding value 88.2 points, kinship coefficient 0.10, diversity index 4.80; ④♀V15 (breeding value 88.7 points, kinship coefficient 0.13, diversity index 4.78) × ♂M01 (breeding value 85.6 points, kinship coefficient 0.07, diversity index 4.72), expected offspring breeding value 87.2 points, kinship coefficient 0.10, diversity index 4.80. Number 4.75; ⑤♀V18 (breeding value 82.3 points, kinship coefficient 0.07, diversity index 4.95) × ♂M05 (breeding value 88.5 points, kinship coefficient 0.09, diversity index 4.82), expected offspring breeding value 85.4 points, kinship coefficient 0.08, diversity index 4.89; ⑥♀V22 (breeding value 87.6 points, kinship coefficient 0.11, diversity index 4.81) × ♂M03 (breeding value 87.3 points, kinship coefficient 0.08, diversity index 4.76), expected offspring breeding value 87.5 points, kinship coefficient 0.09, diversity index 4. .79; ⑦♀V25 (breeding value 86.9, kinship coefficient 0.12, diversity index 4.75) × ♂M07 (breeding value 86.8, kinship coefficient 0.10, diversity index 4.68), expected offspring breeding value 86.9, kinship coefficient 0.11, diversity index 4.72; ⑧♀V30 (breeding value 85.8, kinship coefficient 0.10, diversity index 4.83) × ♂M01 (breeding value 85.6, kinship coefficient 0.07, diversity index 4.72), expected offspring breeding value 85.7, kinship coefficient 0.09, diversity index 4.7 8; ⑨♀V45 (breeding value 81.5 points, kinship coefficient 0.08, diversity index 4.90) × ♂M05 (breeding value 88.5 points, kinship coefficient 0.09, diversity index 4.82), expected offspring breeding value 85.0 points, kinship coefficient 0.09, diversity index 4.86; ⑩♀V52 (breeding value 88.3 points, kinship coefficient 0.13, diversity index 4.79) × ♂M03 (breeding value 87.3 points, kinship coefficient 0.08, diversity index 4.76), expected offspring breeding value 87.8 points, kinship coefficient 0.11, diversity index 4.78.

[0063] The mating program also specifies the technical specifications for offspring rearing: ① During the larval stage, feed high-quality royal jelly (protein content ≥18%) three times a day; ② On the third day after the queen emerges, feed her a sucrose solution containing core probiotic preparations (lactic acid bacteria + bifidobacteria, total concentration 5×10⁻⁶). 7 ③ Regularly monitor the gut microbiota structure and disease resistance of the offspring queen bees, and update the core population by 10%-15% for each generation.

[0064] A one-year follow-up monitoring was conducted on the F1 generation queen bees bred from the newly selected core population according to the mating program. The results showed that: ① Disease resistance: The average survival time of the F1 generation queen bees after viral challenge was 52% longer than that of the ordinary population, the spore load was reduced by 75%, the degree of intestinal pathological damage was reduced by 68%, and the proportion of highly disease-resistant individuals increased from 15% in the ordinary population to 82%; ② Genetic diversity: After three generations of continuous breeding, the average kinship coefficient of the core population remained below 0.15, and no genetic decline was observed. The egg production capacity of the population increased by 23% compared with the original core population; ③ Microbial stability: The Shannon diversity index of the intestinal microbiome of the F1 generation queen bees was 4.75 on average, which was 0.92 higher than that of the ordinary population. The average abundance of lactic acid bacteria and bifidobacteria remained above 30%, and the microbial structure was less affected by environmental fluctuations. In addition, after the core population bred using this method was promoted and applied in five demonstration apiaries, the incidence of microsporidia in bee colonies decreased from 35% to 8%, and the average honey production increased by 18%, achieving significant economic and ecological benefits.

[0065] The same or similar labels correspond to the same or similar parts; The terms used to describe positional relationships in the accompanying drawings are for illustrative purposes only and should not be construed as limiting this patent. Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the implementation of the present invention. For those skilled in the art, other variations or modifications can be made based on the above description. It is neither necessary nor possible to exhaustively list all implementation methods here. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the claims of the present invention.

Claims

1. A method for evaluating queen bee breeding value and selecting core populations based on multi-omics data, characterized in that, include: Obtain multi-omics data of candidate queen bees, including genomic data, gut microbiome data, and host transcriptome data; The multi-omics data were analyzed to screen core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics related to queen bee disease resistance. Based on the core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics, a machine learning model is constructed, and the disease-resistant genome-microbiome joint breeding value of candidate queen bees is calculated through the machine learning model. A comprehensive selection index was constructed by combining the disease-resistant genome-microbiome joint breeding value, the kinship coefficient between the candidate queen bee and the existing core colony members, and the gut microbiome diversity index. Candidate queen bees are ranked according to the comprehensive selection index, and core populations are selected based on the ranking results to generate mating plans.

2. The method according to claim 1, characterized in that, The acquisition of multi-omics data of candidate queen bees includes: Before the artificial infection and challenge experiment with *Apis cerana*, midgut tissue contents of candidate queen bees were collected, and gut microbiome data were obtained by 16S rRNA sequencing. Midgut tissue was collected at specific time points after viral challenge, and host transcriptome data was obtained through transcriptome sequencing. The genomic data of candidate queen bees were obtained by whole-genome sequencing, and the genomic data included whole-genome SNP markers.

3. The method according to claim 1, characterized in that, The analysis of the multi-omics data, screening for core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics related to queen bee disease resistance, includes: The differences in gut microbiota structure between candidate queen bees in the disease-resistant group and the susceptible group before and after viral challenge were compared and analyzed. Multi-omics statistical methods were used to screen core probiotic genera that were significantly positively correlated with high disease resistance and harmful bacterial genera that were related to susceptibility, and the microbiome disease resistance index was calculated. Screening immune-related genes that were significantly differentially expressed in disease-resistant queen bees after challenge with the virus from transcriptome data; By using weighted gene co-expression network analysis, differentially expressed immune genes were associated with differentially abundant microbial taxa to identify key bacterial-gene regulatory pairs and construct microbial-gene interaction network features.

4. The method according to claim 1, characterized in that, The process involves constructing a machine learning model based on the core microbial biomarkers, immune-related genes, and microbial-gene interaction network characteristics. This model is used to calculate the disease-resistant genome-microbiome combined breeding value of candidate queen bees, including: Genomic SNP markers, microbiome disease resistance index before challenge, absolute abundance of key probiotics, and association strength of key bacteria-gene pairs were used as model input features. Survival time, spore load, and intestinal pathological changes of candidate queen bees after viral infection were used as disease resistance phenotypic values, and a machine learning model was constructed using gradient boosting decision tree or random forest algorithms. The model parameters were optimized through cross-validation. The optimized model was then applied to candidate queen bees that had not undergone challenge experiments to predict their potential disease resistance phenotype values, which were used as the joint breeding values ​​for disease resistance genome and microbiome.

5. The method according to claim 1, characterized in that, The formula for calculating the comprehensive selection index is as follows: Where I represents the comprehensive selection index, The values ​​represent the combined breeding values ​​of the disease-resistant genome and microbiome. w1, w2, and w3 are the weights of each indicator and can be adjusted according to the breeding objectives.

6. The method according to claim 1, characterized in that, The comprehensive selection index is constructed by combining the disease-resistant genome-microbiome joint breeding value, the kinship coefficient between the candidate queen bee and existing core colony members, and the gut microbiome diversity index, including: Calculate the average kinship coefficient between the candidate queen bee and the existing core colony members to control the degree of inbreeding; The diversity index of the gut microbiome of candidate queen bees was calculated using the Shannon diversity index formula, and the diversity index is positively correlated with queen bee health and stress resistance.

7. The method according to claim 1, characterized in that, The step of ranking candidate queen bees according to the comprehensive selection index, selecting core populations based on the ranking results, and generating mating plans includes: All candidate queen bees were sorted in descending order based on a comprehensive selection index. The top-ranked candidate queen bees are selected to form the core population. Under the premise of ensuring high disease resistance, the mating scheme is generated by optimizing the algorithm so that the core population maintains a relatively distant kinship at the genetic level and high microbial diversity at the microbiome level.

8. The method according to claim 3, characterized in that, The multi-omics statistical method includes LEfSe analysis, which identifies microbial groups with significant abundance differences between the disease-resistant and susceptible groups, serving as core microbial biomarkers.

9. The method according to claim 4, characterized in that, The training process of the machine learning model includes: dividing multi-omics data into training set and validation set, training the model using the training set, evaluating the model performance using the validation set, adjusting the model parameters based on the evaluation results, until the model prediction accuracy reaches a preset threshold.

Images