Method for screening bee heat stress character related genes based on time sequence transcriptome sequencing and weighted gene co-expression network analysis
Through temporal transcriptome sequencing and weighted gene co-expression network analysis, genes related to honeybee heat stress traits were screened, which solved the problem of lack of honeybee capping heat stress response mechanism, improved the understanding of honeybee heat stress response, and guided bee colony health management and sustainable development of the beekeeping industry.
Patent Information
- Application Number
- CN202510799228.X
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-06-16
- Publication Date
- 2025-09-26
AI Technical Summary
Existing research lacks a systematic analysis of the heat stress response mechanism of honeybees during the capping stage, especially the molecular regulatory network, which affects the physiological development and ecological adaptability of honeybees, resulting in heat stress significantly inhibiting the growth, development and behavioral performance of honeybees, threatening the survival of the bee colony.
Time-sequential transcriptome sequencing and weighted gene co-expression network analysis methods were used to screen genes related to honey bee heat stress traits. By sequencing the transcriptomes of capped honey bees after heat stress treatment at different times, differentially expressed genes were screened, gene co-expression networks were constructed, key gene sets were identified, and screening was carried out in combination with biological characteristics and module characteristics.
It has improved our understanding of honeybees' response to and adaptation to heat stress, analyzed the molecular mechanism of developmental abnormalities caused by heat stress, provided theoretical support for honeybees' environmental adaptability, and guided bee colony health management and the sustainable development of the beekeeping industry.
Smart Images
Figure BDA0005450731500000071 
Figure BDA0005450731500000091 
Figure BDA0005450731500000101
Abstract
Description
Technical Field
[0001] The present invention belongs to the field of biotechnology, and in particular relates to a method for screening genes related to honeybee heat stress traits based on temporal transcriptome sequencing and weighted gene co-expression network analysis. Background Art
[0002] As global climate change intensifies, honeybee populations worldwide are experiencing a significant decline, with native bee species in tropical regions and temperate honeybees facing severe challenges in adapting. This crisis stems in part from the fact that the evolution of honeybee heat tolerance mechanisms lags behind the rate of global warming, rendering physiological defenses less able to withstand environmental stress. Heat stress poses multiple threats to honeybee physiology and ecological adaptability. It not only significantly inhibits individual growth and development and reduces foraging efficiency, but also induces systemic oxidative damage by inducing excessive accumulation of reactive oxygen species (ROS), seriously threatening the survival of the colony. Studies have shown that heat stress significantly impacts honeybee growth and development, leading to prolonged larval development, increased mortality during the capping stage, shortened lifespan, and impacts on foraging and learning. However, existing research has primarily focused on the effects of heat stress on post-hive behavior and the effects of long-term heat stress on honeybee development. Studies on the response mechanisms of Italian worker honeybees to heat stress during the capping stage remain lacking.
[0003] The optimal developmental temperature for capping in Apis mellifera (Apis mellifera) is 29-38°C, while that for Apis cerana (Apis cerana) is 30-38°C. Both species develop normally when cultured at 33-35°C, but bees exposed to temperatures below 29°C or above 39°C exhibit developmental deformities or even die. Two- to three-day-old capping bees are most sensitive to low temperatures; exposure to 20°C for 12 hours results in 50% mortality, with mortality increasing with increasing exposure time. Within the 29-38°C range, the developmental period shortens with increasing temperature when cultured below 35°C, while it increases slightly above 35°C. There are no significant differences in the birth weight of capping worker bees at different developmental temperatures, but the developmental period is shortest for drones at 35°C. High temperature (>35°C) stress results in abnormal pigmentation in the appendages of emerging bees (e.g., pale yellow tarsus of the mid- and hind legs). However, external morphology (snout length, wing width, and dorsal plate length) of bees emerging at 35°C is significantly superior to that of bees exposed at other temperatures, indicating that this temperature is crucial for optimizing economic traits. Under high temperature stress (40°C), all 2-day-old capped Italian honeybees died after only 32 hours of treatment, while 5-day-old individuals had the strongest tolerance. Although honeybees have the ability to adapt to temperature, their ability to regulate body temperature is limited. Continuous or extreme high temperatures will interfere with their development, behavior, physiological functions and reproductive capacity, leading to a decline in colony vitality and weakened immune function. In severe cases, it can lead to Colony Collapse Disorder (CCD). The capping stage is a critical period in the development of honeybees, and individuals at this developmental stage are particularly sensitive to heat stress. However, existing studies have mostly focused on phenotypic observations of the effects of heat stress on capping developmental phenotypes (such as survival rate and emergence rate) and adult bee behavior, while there is still a lack of systematic analysis of its intrinsic response mechanisms (especially molecular regulatory networks). Summary of the Invention
[0004] In order to solve the above technical problems, the present invention provides a method for screening genes related to honeybee heat stress traits based on temporal transcriptome sequencing and weighted gene co-expression network analysis.
[0005] To achieve the above object, the present invention adopts the following technical solutions:
[0006] A method for screening genes related to honey bee heat stress traits based on temporal transcriptome sequencing and weighted gene co-expression network analysis, comprising the following steps:
[0007] S1. After capping the bees, heat stress treatment was performed for different periods of time. After returning to normal temperature and culturing them outside the hive, RNA was extracted from the bees in each group. At the same time, the bees treated normally served as a control group.
[0008] S2. Collect biological characteristics and RNA from honey bee samples from the heat stress treatment group and the control group, sequence them, obtain transcriptome sequencing data, and analyze them against the reference genome;
[0009] S3, differentially expressed genes were screened based on the screening criteria of |log2fold change|≥1 and p-value<0.05;
[0010] S4, reverse transcribe the RNA obtained in step S2 to obtain cDNA, and select several differentially expressed genes for qRT-PCR verification to ensure that the sequencing results are true and reliable;
[0011] S5. Using the WGCNA software package in R language to read the gene expression matrix data obtained by transcriptome sequencing in step S2, perform weighted calculation based on the pick Soft Threshold function to determine the optimal soft threshold;
[0012] S6. Perform sample clustering based on Euclidean distance, calculate the expression similarity coefficient between genes based on the TOM Similarity module, and construct a co-expression network module;
[0013] S7. Display the gene co-expression network model through heat maps, combine the sample biological characteristics and the expression profile characteristics of the module characteristic genes, calculate the correlation coefficient and P value between them, and screen modules related to biological characteristics;
[0014] S8. Core gene screening: To screen for key genes with high correlation with weight and deformity traits, modules with correlation coefficients |r| > 0.5 and P < 0.05 were defined as significantly correlated modules. MM values (module membership) were calculated by calculating the correlation coefficients between genes and module eigengenes, and GS values (gene significance) were calculated by calculating the correlation coefficients between gene expression values and traits. Key genes in the modules were screened based on |MM| > 0.8 and |GS| > 0.2.
[0015] S9. Perform Venn diagram intersection analysis on the S8 module genes and the differentially expressed genes in the S3 newborn bee heat stress treatment group to screen Hub genes.
[0016] In the screening method as described above, optionally, in step S1, the number of days after capping is 2 days.
[0017] In the screening method as described above, optionally, in step S1, the temperature of the heat stress treatment is 42°C; and the temperature of the normal treatment is 35°C.
[0018] In the screening method as described above, optionally, in step S2, the biological characteristics include birth weight, body length, body width, forewing length, forewing width, snout length, hind tibia length, hind femur length and hind tarsus length.
[0019] In the screening method described above, optionally, in step S2, the analysis includes quality control and filtering of the sequencing data, wherein the quality control includes counting the mapping ratio of each sample to the mellifera reference genome (Amel_HAv3.1), and samples with a ratio greater than 80% can be used for subsequent quantitative analysis; the sequencing quality of the obtained raw data files is assessed using FASTQ analysis software, and the relevant quality parameters of the samples, such as the total number of bases with a base recognition accuracy of more than 99.9% (Q30, bp), the percentage of ambiguous bases (N, %), the percentage of bases with a base recognition accuracy of more than 99% (Q20, %), and the percentage of bases with a base recognition accuracy of more than 99.9% (Q30, %), are counted;
[0020] Based on the statistical results, the raw data were filtered to remove adapter sequences, empty reads, low-quality reads, reads that were too long or too short, reads with a copy number of 1, and reads with an average quality score below Q20. High-quality sequences (Clean reads) were obtained through the above filtering.
[0021] In the screening method as described above, optionally, the reference genome refers to Amel_HAv3.1.
[0022] In the screening method described above, optionally, the sequencing results of the heat stress treatment and the control group are compared to select the up-regulated genes and down-regulated genes with large differences.
[0023] The screening method described above, optionally, in step S4, uses qRT-PCR to verify the expression of differentially expressed genes, designs primers for differentially expressed genes, amplifies cDNA, and uses 2 -ΔΔCt The analysis method calculates the relative expression of genes, requiring that the expression trends of key differentially expressed genes in qRT-PCR and sequencing results are consistent, and if both are up-regulated or down-regulated, the results are accurate.
[0024] In the screening method as described above, optionally, in step S8, the deformity determination criterion is the morphological abnormality of the snout, wings or feet being incomplete, bent or folded.
[0025] In the screening method as described above, optionally, in step S8, the correlation coefficient |r| of the significantly correlated modules is greater than 0.5, and P is less than 0.05.
[0026] The beneficial effects of the present invention are:
[0027] The present invention provides a method for screening genes related to honeybee heat stress traits based on temporal transcriptome sequencing and weighted gene co-expression network analysis. By analyzing the traits and genes of honeybees emerging from the hive after heat treatment, temporal transcriptome sequencing is performed, and the sequencing results are associated with the weight and deformity traits of newborn bees, and a gene co-expression network is constructed. Using WGCNA to predict gene interaction networks from quality-checked transcriptome sequencing data, we identified a core set of candidate genes (Pgrp-s2, Tsf1, Apid1, gukh, Vkr, puc, Pdk1, CYP9S1, CYP6AS4, LOC406142, LOC410626, LOC408650, LOC406144, LOC100576758, LOC724126, LOC726803, Rgl, LOC100578437, LOC410614, LOC724644, LOC552797, LOC410370, and LOC724717) that may contribute to differences between the two traits. This method improves our understanding of honey bees' responses and adaptations to capping heat stress and provides a direct theoretical basis for understanding the molecular mechanisms underlying developmental abnormalities (increased malformation rates and reduced body weight) caused by heat stress. This method not only provides theoretical support for analyzing the molecular mechanism of honeybee environmental adaptability, but also provides technical guidance for bee colony health management (such as breeding of heat-resistant bee species and regulation of beehive microenvironment), which is of great significance to the sustainable development of the beekeeping industry in the context of climate change. BRIEF DESCRIPTION OF THE DRAWINGS
[0028] Figure 1 Schematic diagram of the capping process for heat stress treatment.
[0029] Figure 2 The effect of heat stress on the hatching rate and development period of capped rice.
[0030] Figure 3 This is the result graph of the effect of heat stress on the lid deformity rate.
[0031] Figure 4 is the birth weight of newborn bees in different treatment groups.
[0032] Figure 5 The effect of heat stress on the morphological development of the cap.
[0033] Figure 6 The bar graph shows the differentially expressed genes in newborn bees.
[0034] Figure 7 A Venn diagram of differentially expressed genes in newborn bees.
[0035] Figure 8 The results are for the verification of differentially expressed genes.
[0036] Figure 9Soft threshold determination for weighted gene co-expression networks.
[0037] Figure 10 is the gene co-expression module and the number of genes it contains.
[0038] Figure 11 Correlation analysis between modules and traits.
[0039] Figure 12 Identification of Hub genes associated with body weight and malformation phenotypes for the module.
[0040] Figure 13 The figure shows the Venn diagram of genes in the module and differentially expressed genes.
[0041] Note: HS in the figure means heat stress; CW means control group, FW means 4-hour heat stress group, EW means 8-hour heat stress group, and TW means 12-hour heat stress group. DETAILED DESCRIPTION
[0042] Using Italian worker honey bees as their capping model, the study first analyzed the dynamic effects of heat stress on capping phenotypes and physiological parameters. Furthermore, using time-sequential transcriptome sequencing, the study revealed key genes and signaling pathways involved in heat stress responses. This research not only provides theoretical support for understanding the molecular mechanisms of honey bee environmental adaptability but also offers technical guidance for bee colony health management (e.g., breeding heat-tolerant bee strains and regulating the beehive microenvironment), and is of great significance for the sustainable development of the beekeeping industry in the face of climate change.
[0043] The following examples are used to further illustrate the present invention, but should not be construed as limiting the present invention. Without departing from the spirit and substance of the present invention, modifications or substitutions made to the present invention all fall within the scope of the present invention.
[0044] Unless otherwise specified, the technical means used in the examples are conventional means well known to those skilled in the art. Unless otherwise specified, the reagents used in the examples are of analytical grade or above.
[0045] Example 1
[0046] 1. Sample processing:
[0047] Six groups of Italian honey bees with equal strength, good health and no disease were selected. The honeycombs were placed in the bottom box to let the worker bees clean out the queen's egg-laying area, and then the queen-limited egg-laying frame was placed to allow the queen to lay eggs. The capping situation was observed around 8.5 days, and the caps that were capped within 4 hours were recorded as 0 days old. They were cut into 5 parts, packed in gauze bags, and hung vertically on the partition of the incubator. They were then placed in an incubator at 35°C (control group) and 75%±5%RH for culture. When the caps were 2 days old, they were subjected to heat stress treatment at 42°C (heat treatment experimental group) for 4h (FW), 8h (EW), 12h (TW), and 16h (SW). The humidity remained unchanged. After the treatment, they were placed in a 35°C incubator and continued to be cultured until they were released from the hive. The control group was cultured at 35°C until they were released from the hive. Samples of the bees that had left the hive (newborn bees) were collected and their birth weight and development were recorded. The experimental flow chart is as follows. Figure 1 shown.
[0048] Finally, RNA was extracted using the traditional method (TRIZOL), with the following specific steps: (1) 0.1 g of tissue sample (whole bee) was taken, 1 mL of RNA extraction reagent was added, magnetic beads were placed, and the homogenizer was used for 300 sec at 50 Hz. The homogenate was repeatedly blown with a pipette; the homogenate was transferred to a new centrifuge tube, the tissue homogenate was allowed to stand at room temperature for 5 min, and then centrifuged at 12000 g / min at 4°C for 5 min. The supernatant was transferred to a new centrifuge tube;
[0049] (2) Add 200 μL of chloroform, mix thoroughly, and let stand at room temperature for 5 min; centrifuge at 12,000 g, 4°C for 15 min; aspirate the supernatant into a new centrifuge tube, add 500 μL of isopropanol, shake and mix, let stand for 10 min, and then centrifuge at 12,000 g / min, 4°C for 10 min;
[0050] (3) Discard the supernatant and wash the RNA precipitate and the tube wall with 1 mL of 80% pre-cooled ethanol. Centrifuge at 7500g and 4°C for 5 min and discard the supernatant. Open the tube cap and dry the precipitate at room temperature (approximately 5 min). Add 50-100 μL of DEPC water to dissolve the RNA according to the amount of RNA precipitate. Detect the RNA quality by 1.5% agarose gel electrophoresis and the nucleic acid concentration by NanoDrop 2000 spectrophotometer. Qualified samples should be stored at -80°C for future use.
[0051] 1.1 Calculation of occupancy rate
[0052] The hatch rate = the number of normal hatches / the number of sealed lids during heat stress treatment.
[0053] 1.2 Developmental period
[0054] The developmental period is the time it takes for Italian worker bees to emerge from the hive from capping. After the first bee emerges, bees are observed every hour and any deformities are noted.
[0055] The results of the emergence rate and developmental period showed that with the extension of heat stress time, the emergence rate of worker bees showed a downward trend, while the developmental period showed a prolonged trend. Compared with the control group, the emergence rate of capped bees in the heat stress treatment groups of 8h and 12h was significantly reduced, which was reduced by 72.54% and 95.11% respectively compared with the control group (P<0.001). The emergence rate of capped bees in the heat stress treatment group of 16h was 0 (P<0.001). Figure 2 A). Compared with the control group, the developmental period of the heat stress treatment groups of 8h and 12h was significantly prolonged, by 0.47d and 0.97d respectively. Figure 2 B, P < 0.001). The above results indicate that heat stress affects the development of capping by reducing the rate of capping emergence and prolonging the development period of Italian worker bees, and that 16 h is the limit of the time for honey bees to resist heat stress.
[0056] The results of deformity rate showed that no deformity was observed in the newborn bees in the control group and the group treated with heat stress for 4 hours. However, the newborn bees in the group treated with heat stress for 8 hours showed obvious morphological abnormalities, mainly manifested as wing deformities (such as Figure 3 A). Among them, the deformity rate of newborn bees in the group treated with heat stress for 8 hours was 25.56%, and the deformity rate of newborn bees in the group treated with heat stress for 12 hours reached 50% ( Figure 3 B) These results indicate that heat stress can affect the morphological development of newborn bees and may affect subsequent hive activities and behaviors.
[0057] 1.3 Determination of newborn bee birth weight
[0058] While observing the bees emerging from the hive and recording their developmental stages, the newborn weight of the bees emerging from the hive was weighed using an electronic balance and the newborn weight data were recorded. The results are shown in Tables 1 and Figure 4 .
[0059] Table 1 Comparison of birth weight of newborn bees in different groups
[0060]
[0061] The results of the newborn bee birth weight showed that as the heat stress time prolonged, the newborn bee birth weight showed a time-dependent decreasing trend. Among them, the newborn bee birth weight of the heat stress (42℃) treatment group for 8h and 12h was significantly lower than that of the control group (P<0.05, P<0.001), respectively, and was reduced by 18.95mg and 30.05mg ( Figure 4 The above results show that heat stress treatment can reduce the birth weight of bees after the caps are developed, and the effect is the greatest when the treatment is 12 hours.
[0062] 1.4 Determination of morphological development
[0063] At least six bees were collected from each treatment group, placed in 75% alcohol, labeled, and stored in a -20°C freezer. Whole, uniformly sized bees were placed in a Petri dish containing 75% alcohol. The bees were placed under a microscope, and forceps and ophthalmic scissors were used to gently separate the tissues to be measured. Excess tissue was gently peeled off using forceps. To prevent the exfoliated tissue from drying out, the dissected tissues were placed on a glass slide coated with a thin layer of petroleum jelly. A coverslip was then placed over the slide and pressed firmly, squeezing out any excess petroleum jelly. After wiping clean, the two slides were secured with transparent tape to create specimens, labeled, and stored. The specimens were observed under a stereomicroscope according to the measurement criteria proposed by Ruttner in 1978. After adjusting the viewing resolution and setting the scale parameters, images were acquired, classified, labeled, and saved to the corresponding folder to obtain preliminary morphological image data. Data statistics and measurements were performed on the images of the Italian honey bee specimens using Image J software. In this study, eight morphological indices were measured, including body length, body width, forewing length, forewing width, snout length, hind tibia length, hind femur length, and hind tarsus length, to evaluate the effects of heat stress on the capping morphology of Italian worker bees. Figure 5 As shown in Figure 2, as the heat stress time prolonged, most morphological indices showed a gradual downward trend, and the extent of the decline was positively correlated with the treatment time. Figure 5 AH (body length, body width, forewing length, forewing width, snout length, hind tibia length, hind femur length and hind tarsus length) of the newborn bees in the 4-h heat stress treatment group had no significant effect (P>0.05); the body length, body width, snout length, hind tibia length, hind femur length and hind tarsus length of the newborn bees in the 8-h heat stress treatment group were significantly reduced (P<0.01), and the various indicators were reduced by 1.32 mm (body length), 0.57 mm (body width), 0.67 mm (body width) and 0.71 mm (body width) respectively compared with the control group. The results showed that heat stress significantly reduced the morphological development of Italian worker bees after capping and emergence, with the body length, body width, forewing length, forewing width, foreskin length, foreskin length, foreskin width, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length, foreskin length
[0064] 2. Transcriptome sequencing: Samples of honey bees from the 35°C (control group) and 42°C (heat stress experimental group) treatment groups were sent to Nanjing Paisonno Gene Technology Co., Ltd. for RNA-seq. Specifically, the sequencing data were first quality controlled and filtered. Quality control included counting the mapping ratio of each sample to the Italian honey bee reference genome. Samples with a ratio greater than 80% could be used for subsequent quantitative analysis. The sequencing quality of the obtained raw data files was assessed using FASTQ analysis software. The relevant quality parameters of the samples, including the total number of bases with a base recognition accuracy of more than 99.9%, the percentage of ambiguous bases, the percentage of bases with a base recognition accuracy of more than 99%, and the percentage of bases with a base recognition accuracy of more than 99.9%, were statistically analyzed.
[0065] Based on the statistical results, the original offline data were filtered: S21, remove adapter sequences; S22, remove empty reads; S23, remove low-quality reads; S24, remove reads that are too long or too short; S25, remove reads with a copy number of 1; S26, remove reads with an average quality score lower than Q20; after the above screening, high-quality sequences Cleanreads were finally obtained.
[0066] Based on the screening criteria of |log2fold change|≥1 and p-value<0.05, differentially expressed genes were screened. In the comparison between the heat stress group and the control group, a total of 194 up-regulated genes and 96 down-regulated genes were detected (e.g. Figure 6 ), a total of 225 differentially expressed genes were obtained in the heat stress treatment group (e.g. Figure 7 ), where Fold change refers to the ratio of the gene FPKM value of the treatment group to that of the control group.
[0067] 3. Transcriptome Sequencing Data Verification: Eight genes were selected from the screened genes with high fold-difference to validate the sequencing results. Sequences of the relevant target genes were retrieved from the National Center for Bioinformation (NCBI) website. Primers were designed using Premier 5.0 software and synthesized by the Wuhan branch of Beijing Qingke Biotechnology Co., Ltd. Sequence information is shown in Table 2.
[0068] Table 2 Primer information
[0069]
[0070]
[0071] The RNA samples that passed the quality inspection at 35°C (control group) and 42°C (heat treatment experimental group) were reverse transcribed into cDNA using the reverse transcription kit (Cat. No.: AG11728) of Hunan Aikerui Bioengineering Co., Ltd. and stored at -20°C for future use. Using the primers in Table 2, the cDNA obtained by reverse transcription of the RNA of the newborn bee samples in the control group and the heat treatment experimental group was used as a template for RT-qPCR amplification. The PCR reaction system was 10 μL: 1 μL cDNA, 5 μL 2×SYBR Green Pro Taq HS Premix, 0.2 μL each of upstream and downstream primers, and ddH2O was added to make up the total volume to 10 μL. After the system was prepared, the real-time fluorescence quantitative PCR instrument of the American Bio-Rad Company was started to detect. The program was set as follows: 95°C pre-denaturation for 30s; 95°C denaturation for 15s, 60°C annealing for 30s, 40 cycles; the melting curve used the default settings of the instrument; and the results were exported. After the amplification was completed, the corresponding Ct value was exported from the program of the PCR instrument, and 2 -ΔΔCt The relative expression level of the target gene in each sample was calculated by the analytical method. The calculation formula is: ΔCt = Ct (target gene) - Ct (reference gene), ΔΔCt = ΔCt (treatment group) - ΔCt (control group). Based on the quantitative results of the gene, the Log2fold change value was calculated and finally compared with the Log2fold change value of the related gene calculated based on the sequencing results. The quantitative results were consistent with the sequencing results, proving that the sequencing results were reliable. Figure 8 shown.
[0072] 2.2 Number of differentially expressed genes in different groups: The number of differentially expressed genes in each comparison combination was displayed in a bar graph. Compared with the control group, the number of up-regulated genes (meaning the expression level of a gene in the treatment group was significantly higher than that in the control group) was 18, 77, and 99, respectively, and the number of down-regulated genes (meaning the expression level of a gene in the treatment group was significantly lower than that in the control group) was 23, 39, and 34, respectively. Figure 6 shown.
[0073] 2.3 Transcriptome sequencing Venn diagram: The differentially expressed genes in each group were plotted in a Venn diagram to obtain the number of differentially expressed genes after heat stress treatment (excluding the common differentially expressed genes). The co-expression Venn diagram shows the number of uniquely expressed genes in each group, and the overlapping area shows the number of co-expressed genes in each treatment group. There are 225 differentially expressed genes between the heat stress group and the control group, such as Figure 7 shown.
[0074] 2.4WGCNA Analysis
[0075] Weighted gene co-expression network construction (weighted gene co-expression network analysis, WGCNA) has been widely used in the downstream analysis of transcriptome sequencing. This method can identify genes with the same expression trend in the same group of samples, and associate these genes with biological phenotypes in a modular form, so as to achieve the purpose of identifying major regulatory genes, screening functional candidate genes, or assisting in the annotation of gene functions. The present invention compares the sequencing data of 16 newborn bee samples to the Italian honey bee reference genome (Amel_HAv3.1) in the transcriptome data quality control link, associates the sequencing results with the newborn bee weight and deformity traits, and constructs a gene co-expression network. WGCNA is used to predict the gene interaction network from the transcriptome sequencing data that has passed the quality inspection, and finds the core candidate gene set that may cause the difference between the two traits. The specific operations are as follows:
[0076] (1) The WGCNA package (Langfelder and Horvath 2008) in R language was used to read the gene expression matrix data obtained by transcriptome sequencing. The weighted calculation was performed based on the pick Soft Threshold function, and the optimal soft threshold was determined to be 7. The results are shown in Figure 2. Figure 9 shown.
[0077] (2) Samples were clustered based on Euclidean distance, expression similarity coefficients between genes were calculated based on the TOM Similarity module, and a co-expression network module was constructed. The results are shown in Figure 2. Figure 10 shown.
[0078] (3) The gene co-expression network model is displayed through a heat map. The biological characteristics of the sample and the expression profile characteristics of the module characteristic genes are combined to calculate the correlation coefficient and P value between them, and the modules related to the biological characteristics are screened. The results are as follows Figure 11 shown.
[0079] (4) Core gene screening: To screen key genes with high correlation with weight and deformity traits, modules with correlation coefficient |r|>0.5 and P<0.05 were defined as significantly correlated modules. The MM value (module membership) was obtained by calculating the correlation coefficient between genes and module characteristic genes, and the GS value (gene significance) was obtained by calculating the correlation coefficient between gene expression values and traits. Key genes in the module were screened according to |MM|>0.8 and |GS|>0.2. The results are as follows: Figure 12 shown.
[0080] (5) Venn diagram intersection analysis was performed with the differentially expressed genes of the newborn bees to screen Hub genes. The results are as follows Figure 13 shown.
[0081] The results showed that the magenta module was significantly positively correlated with the weight of newborn bees and significantly negatively correlated with deformities, while the red module was significantly negatively correlated with weight and significantly positively correlated with deformities. The magenta module contained 180 genes, and the red module contained 321 genes.
[0082] Venn diagram intersection analysis of red and magenta module genes and 225 differentially expressed genes based on WGCNA screening ( Figure 13 ), a total of 23 core genes were identified: Pgrp-s2, Tsf1, Apid1, gukh, Vkr, puc, Pdk1, CYP9S1, CYP6AS4, LOC406142, LOC410626, LOC408650, LOC406144, LOC100576758, LOC724126, LOC726803, Rgl, LOC100578437, LOC410614, LOC724644, LOC552797, LOC410370, and LOC724717.
[0083] The key core genes for heat stress traits identified in this screening may be directly involved in the developmental regulation of honey bee capping under heat stress, providing a direct theoretical basis for understanding the molecular mechanisms underlying heat stress-induced developmental abnormalities (increased deformity rates and weight loss). The expression patterns of these core genes may serve as molecular markers for assessing honey bee heat tolerance. They can also provide technical guidance for bee colony health management (such as breeding heat-tolerant bees and regulating the beehive microenvironment), which is of great significance for the sustainable development of the beekeeping industry in the context of climate change.
Claims
1. A method for screening genes related to honey bee heat stress traits based on temporal transcriptome sequencing and weighted gene co-expression network analysis, characterized in that: It includes the following steps: S1. After capping the bees, heat stress treatment was performed for different periods of time. After returning to normal temperature and culturing them outside the hive, RNA was extracted from the bees in each group. At the same time, the bees treated normally served as a control group. S2. Collect biological characteristics and RNA from honey bee samples from the heat stress treatment group and the control group, sequence them, obtain transcriptome sequencing data, and analyze them against the reference genome; S3, differentially expressed genes were screened based on the screening criteria of |log2fold change|≥1 and p-value<0.05; S4, reverse transcribe the RNA obtained in step S2 to obtain cDNA, and select several differentially expressed genes for qRT-PCR verification to ensure that the sequencing results are true and reliable; S5. Using the WGCNA software package in R language to read the gene expression matrix data obtained by transcriptome sequencing in step S2, perform weighted calculation based on the pick Soft Threshold function, and determine the optimal soft threshold to be 7; S6. Perform sample clustering based on Euclidean distance, calculate the expression similarity coefficient between genes based on the TOM Similarity module, and construct a co-expression network module; S7. Display the gene co-expression network model through heat maps, combine the sample biological characteristics and the expression profile characteristics of the module characteristic genes, calculate the correlation coefficient and P value between them, and screen modules related to biological characteristics; S8. Core gene screening: To screen key genes with high correlation with weight and deformity traits, modules with correlation coefficient |r|>0.5 and P<0.05 were defined as significantly correlated modules; The MM value was obtained by calculating the correlation coefficient between the gene and the module characteristic gene, and the GS value was obtained by calculating the correlation coefficient between the gene expression value and the trait; the key genes in the module were screened according to |MM|>0.8 and |GS|>0.2; S9. Perform Venn diagram intersection analysis on the S8 module genes and the differentially expressed genes in the S3 newborn bee heat stress treatment group to screen Hub genes.
2. The screening method according to claim 1, wherein In step S1, the number of days after capping is 2 days.
3. The screening method according to claim 1, wherein In step S1 , the temperature for heat stress treatment is 42° C.; the temperature for normal treatment is 35° C.
4. The screening method according to claim 1, wherein In step S2, the biological characteristics include birth weight, body length, body width, forewing length, forewing width, snout length, hind tibia length, hind femur length and hind tarsus length.
5. The screening method according to claim 1, wherein In step S2, the analysis includes quality control and filtering of the sequencing data, wherein the quality control includes counting the mapping ratio of each sample to the reference genome of the Apis mellifera ligustica. Samples with a ratio greater than 80% can be used for subsequent quantitative analysis; the sequencing quality of the obtained raw data files is assessed using FASTQ analysis software, and the relevant quality parameters of the samples, such as the total number of bases with a base recognition accuracy of more than 99.9%, the percentage of ambiguous bases, the percentage of bases with a base recognition accuracy of more than 99%, and the percentage of bases with a base recognition accuracy of more than 99.9%, are counted; Based on the statistical results, the raw data were filtered to remove adapter sequences, empty reads, low-quality reads, reads that were too long or too short, reads with a copy number of 1, and reads with an average quality score below Q20. After the above filtering, high-quality clean reads were finally obtained.
6. The screening method according to claim 1, wherein The reference genome refers to Amel_HAv3.
1.
7. The screening method according to claim 1, wherein After comparing the sequencing results of heat stress treatment with those of the control group, up-regulated and down-regulated genes with significant differences were selected.
8. The screening method according to claim 1, wherein In step S4, qRT-PCR was used to verify the expression of differentially expressed genes, primers for differentially expressed genes were designed, and cDNA was amplified using 2 -ΔΔCt The analysis method calculates the relative expression of genes, requiring that the expression trends of key differentially expressed genes in qRT-PCR and sequencing results are consistent, and if both are up-regulated or down-regulated, the results are accurate.
9. The screening method according to claim 1, wherein In step S9, the deformity determination standard is the abnormal morphology of the snout, wings or feet, such as incompleteness, bending or folding.
10. The screening method according to claim 1, wherein In step S9 , the correlation coefficient |r| of the significantly correlated modules is greater than 0.5, and P is less than 0.05.