A high-throughput sequencing-based screening method for crab aggression genes

Abstract

A method for screening crab aggressive genes based on high-throughput sequencing, comprising evaluating the aggressiveness of crabs through mirror experiments and an aggressive evaluation model, and selecting crabs with different aggressive phenotypes; using transcriptome analysis to identify different genes of the thoracic ganglion of crabs with different aggressive phenotypes, with no less than 3 biological repeats in each group; performing KEGG enrichment analysis on the different genes to screen the first n different signal pathways to obtain a subset of aggressive different genes; through gene set enrichment analysis, screening genes in the subset of different genes that contribute greatly to the enrichment score as a target gene set; using expression correlation analysis to obtain a key gene set for the regulation of aggressive phenotypes; and taking the receptor protein in the key gene set as a core regulatory protein. The present application mines core genes related to aggressiveness through hierarchical screening, reveals and verifies the target protein-ligand signal pathway for the regulation of aggressiveness, and provides theoretical and practical support for the application of molecular breeding technology in crustaceans.

Description

TECHNICAL FIELD

[0001] The present application relates to a crab aggressive gene screening method based on high-throughput sequencing, and belongs to the field of animal behavior and molecular biology. BACKGROUND

[0002] Aggression is an important personality trait that is ubiquitous in the animal kingdom and has profound effects on resource acquisition, growth, life history, and fitness. Individuals with different levels of aggression often exhibit behaviors that match their competitive abilities, such as strong aggressive individuals exhibiting competitive and resource monopolizing behaviors, and weak aggressive individuals avoiding direct competition for resources and costly social interactions.

[0003] Aggression, as a complex quantitative trait, is influenced by many genes. For example, there are 12 gene modules that are significantly associated with the duration of aggression in zebrafish (Danio rerio); cytochrome P450 gene family and genes related to electron transport, voltage-gated potassium channels, catabolism, G protein-coupled receptor signaling, and basic cellular and metabolic processes are all associated with aggressive behavior. A large number of gene expression patterns can be seen as a portrait of the neural genomic state that controls aggression, with different aggression phenotypes corresponding to different gene expression patterns in the central nervous system. Understanding these mechanisms not only helps to reveal the behavioral characteristics of crustaceans in depth, but also provides a theoretical basis for crustacean breeding and resource management.

[0004] Under traditional pond culture conditions, crabs often face problems such as limited space, food, and repeated changes in group composition, and they are more likely to engage in intraspecific cannibalism than in natural environments, leading to limb mutilation and even death. Mutilated individuals have difficulty feeding, so their growth rate is often affected, leading to a decrease in the quality of aquaculture and a reduction in production. Therefore, it is necessary to breed new varieties of crabs based on their aggressive traits. The present application combines behavioral and molecular biology research to establish a crab aggressive gene screening method, which identifies key regulatory genes through a hierarchical screening mechanism, providing new ideas for future crab germplasm innovation and new variety breeding. SUMMARY

[0005] The purpose of the present application is to apply behavioral and molecular biology techniques to provide a method for screening aggressive genes in crabs, providing molecular technical support for the breeding of excellent crab varieties.

[0006] To achieve the above-mentioned application purposes, the technical solutions adopted by the present application are as follows:

[0007] A method for screening aggressive genes in crabs based on high-throughput sequencing, comprising:

[0008] Step 1, evaluate the aggressiveness of crabs by mirror experiment and aggressive assessment model, select crabs with different aggressiveness phenotypes;

[0009] The features further include the following steps:

[0010] Step 2, use transcriptomic analysis to identify differentially expressed genes in the thoracic ganglion of crabs with different aggressiveness phenotypes, with no less than 3 biological replicates per group;

[0011] Step 3, perform KEGG enrichment analysis on the differentially expressed genes to screen the top n differentially expressed signaling pathways to obtain a subset of aggressiveness differentially expressed genes;

[0012] 1. Perform KEGG pathway enrichment analysis on the differentially expressed genes obtained in Step 2 using R language;

[0013] 2. Perform Fisher's exact test on the above analysis results, and when P < 0.05, it indicates that there is significant enrichment of the KEGG pathway function, and the number of significantly enriched KEGG pathways is denoted as Q1;

[0014] 3. Sort the significantly enriched KEGG pathways in ascending order of P value, and the number of genes in the final differentially expressed gene subset is denoted as Q2:

[0015] Case 1: Q1 ≥ 20

[0016] If n = 20 and Q2 ≤ 500, then select the genes enriched in the top 20 differentially expressed signaling pathways as the differentially expressed gene subset;

[0017] If n = 20 and Q2 > 500, then calculate the genes enriched in the top n-1 differentially expressed signaling pathways until Q2 ≤ 500; Case 2: Q1 < 20

[0018] If n = Q1 and Q2 ≤ 500, then select the genes enriched in the top Q1 differentially expressed signaling pathways as the differentially expressed gene subset;

[0019] If n = Q1 and Q2 > 500, then calculate the genes enriched in the top n-1 differentially expressed signaling pathways until Q2 ≤ 500.

[0020] Step 4, screen genes in the differentially expressed gene subset that contribute significantly to the enrichment score using gene set enrichment analysis (GSEA) to obtain a target gene set:

[0021] 1. If Q2 ≤ 50, then all genes in the differentially expressed gene subset are the target gene set, and proceed directly to Step 5;

[0022] If Q2 > 50, then perform GSEA analysis on the differentially expressed gene subset using R language;

[0023] 2. The screening criteria for GSEA are |NES|>1 and FDR<0.25, where NES is the standardized enrichment score and FDR is the p-value after multiple hypothesis testing correction; the larger the |NES|, the smaller the FDR value, and the more significant the enrichment.

[0024] 3. For genes in the differential gene subset, the Signal2Noise sorting algorithm is used based on the NES score, and the size of the prior gene set (set according to the number of genes in the differential gene subset) is 15 to Q2;

[0025] 4. Select the gene combination that contributes the most to the enrichment score as the target gene set.

[0026] Step 5: Use expression level correlation analysis to obtain the key gene set regulating the aggressive phenotype;

[0027] Step 6: Focus on receptor proteins with a concentration of key genes as core regulatory proteins.

[0028] In step 1, the attack assessment is performed in the following manner:

[0029] 1. Isolate and temporarily keep the crab in an aquarium for two weeks;

[0030] 2. Select healthy crabs with intact appendages and in the intermolt period for mirror imaging experiments, and repeat the images 3 times, with an interval of no less than 24 hours between each image;

[0031] 3. Record and quantify the aggressive behavior of crabs using behavior analysis software, including the cumulative duration of stillness and the relative distance traveled (the ratio of the distance traveled to the width of the carapace);

[0032] 4. Based on the crab aggression assessment model, calculate the aggression score to obtain crabs with different aggression phenotypes;

[0033] In step 1, the mirror experiment is conducted in the following manner:

[0034] Select an observation cylinder with an opaque white inner wall. Fix the camera directly above the observation cylinder. Take two mirrors with a width equal to the inner diameter of the observation cylinder, and place them vertically in the middle of the observation cylinder with their backs touching, so that the two sides can be mirrored simultaneously. Figure 2 3) Insert a partition on each side of the mirror in the observation tank, and place two crabs on one side of the partition respectively, separating them from the mirror. After acclimatizing for 10 minutes, gently remove the partition, exposing the crabs to the mirror. Film for 20 minutes, then stop. After each group of films, change the seawater in the observation tank and clean the tank. Filming should be conducted under completely dark conditions, and the seawater temperature and salinity should be kept consistent with the temporary holding conditions during the experiment.

[0035] In step 1, the aggressive behavior of the crabs was recorded and quantified using behavioral analysis software as follows:

[0036] The detection sampling rate is set to 8 frames / second, and the detection method is either dynamic silhouette or grayscale gradient. Lost frame correction and smooth trajectory settings are enabled. The trajectory smoothing configuration is as follows:

[0037] a. Smoothing is performed based on 10 samples before and after each sampling point;

[0038] b. If the moving distance is less than 3cm, the sampling point will still be set to the previous position;

[0039] c. If the maximum distance moved is greater than 20cm, the sampling point will be set as missing.

[0040] In step 1, the aggression assessment model is:

[0041] Y = 0.023X1 - 0.001X2 - 0.002

[0042] Where Y is the aggression score, X1 is the relative movement distance, and X2 is the cumulative duration of stationary movement.

[0043] In step 1, after conducting the attack assessment,

[0044] After anesthetizing the crab with cryotherapy for 5 minutes, the thoracic ganglion was quickly dissected on ice and placed in a 2ml cryovial. It was then immediately placed in liquid nitrogen for brief flash freezing and subsequently stored in a -80℃ freezer for later testing.

[0045] Step 2 includes collecting raw sequencing data and transcriptomics analysis:

[0046] The raw sequencing data collected includes:

[0047] 1. Randomly select at least 2 crabs of each aggressive phenotype from step 1 for pooling, and set at least 3 biological replicates for each group;

[0048] 2. Total RNA was extracted using the TRIzol method. RNA integrity was determined by 1% agarose gel electrophoresis. The concentration and purity of total RNA were detected using a Nanodrop 2000 nucleic acid and protein analyzer to ensure that the A260 / A280 ratio was between 1.8 and 2.0.

[0049] 3. After extracting 1 μg of total RNA from the sample, base pairing was performed with the poly A at the 3' end of the mRNA using magnetic beads with Oligo. After purifying the mRNA, fragmentation buffer was added to randomly break the mRNA into small fragments of about 300 bp.

[0050] 4. Using mRNA as a template, one-stranded cDNA was synthesized by reverse polymerization, followed by the synthesis of two-stranded cDNA. End-repair mix was added to complete the ends, and then an A base was added to the 3' end. The product after adapter ligation was then purified and fragment sorted. The sorted product was used for PCR amplification, and the final purified library was obtained.

[0051] 5. The constructed library was sequenced using Illumina HiSeq Xten / NovaSeq6000 to obtain raw sequencing data.

[0052] Transcriptomics analysis includes:

[0053] 1. Use SeqPrep and Sickle software to filter the raw sequencing data to obtain high-quality sequencing data (clean reads); 2. Use HiSat2 and TopHat2 software to align the quality-controlled clean reads with the reference genome to obtain mapping data, and perform quality assessment on the alignment results;

[0054] 3. Gene expression levels were calculated by the number of clean reads located in the genomic region. RSEM software was used to perform quantitative analysis of gene expression levels. The quantitative results were expressed in FPKM (female fraction of kilobases).

[0055] 4. Differential gene expression levels were analyzed using DESeq2 software to identify differentially expressed genes, and then genes with significant differential expression were screened. The screening criteria were: FDR < 0.05 & |log2FC| ≥ 2, where FDR is the p-value after multiple hypothesis testing.

[0056] In step 5, obtaining the key gene set regulating the aggressive phenotype using expression level correlation analysis includes the following steps:

[0057] 1. Use R language to perform expression correlation analysis on the target gene set, generate a correlation matrix, where each element in the matrix represents the correlation between two genes, and construct a co-expression matrix;

[0058] 2. Spearman correlation analysis was performed on the co-expression matrix, and Benjamini-Hochberg multiple tests were used to correct the analysis results. The significance level was P<0.05, and the absolute threshold for correlation was 0.8.

[0059] 3. Calculate the connectivity of each gene in the target gene set and denote its absolute value as D; connectivity refers to the sum of the correlations between a gene and other genes. Key genes in the network are identified by the connectivity of each gene.

[0060] 4. Calculate the mean D value for all genes using the following formula:

[0061] DM = (D1 + D2 + ... + Dq) / q

[0062] Where DM is the mean of D, D1,…,Dq are gene connectivity, and q is the number of genes in the target gene set that have significant correlations.

[0063] 5. Use R packages (such as WGCNA, igraph) to visualize the gene co-expression matrix;

[0064] 6. Using D≥DM as the screening criterion, obtain the key gene set regulating aggressive phenotype.

[0065] The application of the high-throughput sequencing-based crab aggressive gene screening method in crab germplasm innovation and new variety breeding.

[0066] The present invention has the following advantages:

[0067] 1. A novel method for screening aggressive genes in crabs based on high-throughput sequencing was proposed for the first time, employing hierarchical screening (…). Figure 1 This study aims to identify core genes related to aggression, locate the core receptor proteins regulating aggression in the central nervous system of crabs, and reveal and verify the target protein-ligand signaling pathway of aggression regulation, providing theoretical and practical support for the application of molecular breeding technology in crustaceans.

[0068] 2. The correlation between aggressive phenotype and core genes was established, and molecular markers for identifying aggressive phenotypes in crabs were identified and verified, providing new ideas for future crab germplasm innovation and new variety breeding. Attached Figure Description

[0069] Figure 1 A flowchart for screening gene hierarchies of aggressive core regulatory genes.

[0070] Figure 2 This is a diagram illustrating the collection of data for aggressive behavior.

[0071] Figure 3 Capture a top-down view for aggressive behavior.

[0072] Figure 4 Transcriptomic results of the thoracic ganglia of the swimming crab with different aggressive phenotypes.

[0073] In this diagram, A is a differential gene volcano plot, where red represents significantly upregulated genes, blue represents significantly downregulated genes, and gray represents no significant difference; B is a PCA analysis diagram between samples, where GS represents the highly aggressive group and GW represents the weakly aggressive group; and C is the top 20 pathways in the KEGG enrichment analysis.

[0074] Figure 5 The GSEA results are shown for differential gene subsets of the swimming crab with different aggressive phenotypes.

[0075] Among them, AD represents the top 4 enriched pathways of GSEA, GS represents the highly aggressive group, and GW represents the weakly aggressive group.

[0076] Figure 6 A diagram outlining the neural regulation pathway of the NMDAR-CaMKⅡ pathway in the swimming crab *Portunus trituberculatus*.

[0077] Figure 7 To illustrate the differences in aggression scores of the three-spined swimming crab under different experimental treatments.

[0078] Figure 8 This image shows the localization and fluorescence intensity of NMDAR in the thoracic ganglia of the swimming crab *Portunus trituberculatus*.

[0079] Among them, NMDAR, DAPI, and FITC are different fluorescent dyes:

[0080] NMDAR (N-methyl-D-aspartic acid receptor) fluorescent dye,

[0081] DAPI (4',6-Diamidino-2-Phenylindole, diphenylaminoacridine) fluorescent dye,

[0082] Merge display: Combines images of two fluorescence channels.

[0083] Figure 9 Verification results for the NMDAR-CaMKⅡ signal path.

[0084] Where A represents the difference in average fluorescence intensity of NMDAR protein between the two aggressive phenotypes under different experimental treatments; BL represents the difference in relative expression levels of 11 key genes; the significance of differences between the strong aggressive groups under different experimental treatments is represented by the corresponding lowercase letters, and the significance of differences between the weak aggressive groups is represented by the corresponding uppercase letters. Detailed Implementation

[0085] The present invention will now be described in detail below using the swimming crab (Portunus trituberculatus) as an example, in conjunction with the accompanying drawings. The described experiments are only some embodiments of the present invention, not all embodiments. All other embodiments obtained in the art based on the embodiments of the present invention fall within the scope of protection of the present invention.

[0086] Step 1: Assess the aggression of the swimming crab *Portunus trituberculatus* using mirror experiments and an aggression assessment model, selecting crabs with both strong and weak phenotypes:

[0087] 1. Isolate and temporarily keep the male three-spined swimming crab (Portunus trituberculatus) in an aquarium (170×80×20cm) for two weeks;

[0088] 2. Select healthy swimming crabs with intact appendages and in the intermolting period (carapace width 91.32±4.27mm, n=146) for mirror imaging experiments, repeat the images 3 times, with an interval of 24h between each time;

[0089] 3. Mirroring Experiment: Select an observation tank with an opaque white inner wall. Fix the camera directly above the tank. Take two mirrors with a width equal to the inner diameter of the tank and place them vertically in the center of the tank with their backs touching, so that the two sides can be mirrored simultaneously. Insert a partition into each side of the mirror in the tank, and place two crabs on one side of the partition, separating the crabs from the mirrors. Figure 2 ,3);

[0090] 4. After fasting for 24 hours, place the swimming crabs in the observation tank. After acclimatizing for 10 minutes, gently remove the partition to expose the crabs to the mirror. Film for 20 minutes, then save the video for later analysis.

[0091] 5. After each group of photos is taken, replace the seawater in the observation tank and clean the tank. All photos should be taken in a completely dark environment, and the seawater temperature and salinity should be maintained consistent with the temporary holding conditions throughout the experiment.

[0092] 6. Record and quantify the aggressive behavior of swimming crabs using behavior analysis software, including the cumulative duration of stillness and the relative movement distance (the ratio of movement distance to carapace width). Set the detection sampling rate to 8 frames / second, and the detection method to dynamic silhouette or grayscale gradient. Enable lost frame correction and trajectory smoothing settings; the trajectory smoothing configuration is as follows:

[0093] a. Smoothing is performed based on 10 samples before and after each sampling point;

[0094] b. If the moving distance is less than 3cm, the sampling point will still be set to the previous position;

[0095] c. If the maximum distance moved is greater than 20cm, the sampling point will be set as missing.

[0096] 7. Based on the crab aggression assessment model, aggression scores are calculated to obtain swimming crabs with strong and weak aggression phenotypes. The aggression assessment model is as follows:

[0097] Y = 0.023X1 - 0.001X2 - 0.002

[0098] Where Y is the aggression score, X1 is the relative movement distance, and X2 is the cumulative duration of stationary movement.

[0099] 8. After the mirror image experiment was completed, the swimming crabs were returned to the temporary holding aquarium to recover for 24 hours. Then, the swimming crabs were cryo-anesthetized for 5 minutes, and the thoracic ganglia were quickly dissected on ice and placed in 2ml cryovials for short-term flash freezing.

[0100] It was then stored in a -80℃ freezer for later testing.

[0101] Step 2: Using transcriptomics analysis, detect and identify differentially expressed genes in the thoracic ganglia of the swimming crab with strong and weak aggression phenotypes:

[0102] 1. Randomly select 2 swimming crabs of each aggressive phenotype from step 1 for pooling, with 3 biological replicates for each group; 2. The thoracic ganglia of crabs are located in the cephalothorax and are connected to the cerebral ganglia. They send nerves to the appendages and participate in coordinating movement and local reflex activities. Therefore, we selected the thoracic ganglia as the detection site.

[0103] 3. RNA was extracted using the TRIzol method, and RNA integrity was determined by 1% agarose gel electrophoresis. The concentration and purity of total RNA were detected using a Nanodrop2000 nucleic acid and protein analyzer to ensure that the A260 / A280 ratio was between 1.8 and 2.0.

[0104] 4. After extracting 1 μg of total RNA, use magnetic beads with Oligo(dT) to pair the bases with the poly A at the 3' end of the mRNA. After purifying the mRNA, add fragmentation buffer to randomly break the mRNA into small fragments of about 300 bp.

[0105] 5. Using mRNA as a template, reverse synthesize one-stranded cDNA, then synthesize two-stranded cDNA, add end repair reagent to make it blunt ends, and then add a base A to the 3' end;

[0106] 6. Next, the product after ligation of the adapter is purified and fragment sorted. The sorted product is used for PCR amplification, and the final library is obtained after purification.

[0107] 7. The constructed library was sequenced using Illumina HiSeq Xten / NovaSeq6000 to obtain raw sequencing data;

[0108] 8. Use SeqPrep and Sickle software to filter the raw sequencing data to obtain high-quality sequencing data (clean reads). Use HiSat2 and TopHat2 software to align the quality-controlled clean reads with the reference genome (GCF_017591435.1) to obtain mapping data. Perform quality assessment on the alignment results.

[0109] 9. Gene expression levels were calculated by the number of clean reads located in genomic regions. RSEM software was used for quantitative analysis of gene expression levels, with the results expressed as FPKM (female fractional kbps).

[0110] Unit;

[0111] 10. Differential gene expression levels were analyzed using DESeq2 software to identify differentially expressed genes, and then genes with significant differential expression were screened. The screening criteria were: FDR < 0.05 & |log2FC| ≥ 2, where FDR is the p-value after multiple hypothesis testing.

[0112] 11. Results showed that transcriptome sequencing of six samples from the strong and weak aggression groups of *Portunus trituberculatus* yielded 39.43 Gb of clean reads, with each sample having over 6.07 Gb of clean reads and a Q30 base percentage exceeding 95.89%. A total of 3850 differentially expressed genes (DEGs) were identified between the strong and weak aggression groups, of which 2066 were significantly upregulated and 1784 were significantly downregulated (P < 0.05). Figure 4 A, B).

[0113] Step 3: Perform KEGG enrichment analysis on the differentially expressed genes, and screen the top 20 differentially expressed signaling pathways to obtain a subset of aggressive differentially expressed genes:

[0114] 1. The differentially expressed genes obtained in step 2 were subjected to KEGG pathway enrichment analysis using R language. Fisher's exact test was used for the test. When P < 0.05, it was found that there was significant enrichment of the KEGG pathway function.

[0115] 2. The results showed that the top 20 significantly enriched signaling pathways were mainly concentrated in energy metabolism, such as oxidative phosphorylation, PI3K-Akt, AMPK, and FoxO signaling pathways; and neurodegenerative diseases, such as Alzheimer's disease, Parkinson's disease, Huntington's disease, and neurodegenerative pathways in various other diseases. Figure 4 C);

[0116] 3. Among the top 20 significantly enriched signaling pathways, 257 differentially expressed genes were selected that are associated with energy metabolism and neurodegenerative diseases, which constitutes the subset of aggressive differentially expressed genes.

[0117] Step 4: Using Gene Set Enrichment Analysis (GSEA), select genes from the gene subset that contribute significantly to the enrichment score as the target gene set.

[0118] 1. Gene set enrichment analysis (GSEA) was performed on differentially expressed gene subsets using R language. The prior gene set was set to 15-257, and the sorting algorithm was Signal2Noise.

[0119] 2. Based on the screening criteria of |NES|>1 and FDR<0.25, the key signaling pathways enriched in the differentially expressed gene subsets were identified as oxidative phosphorylation, retrograde endocannabinoid signaling, axonal regeneration, and cAMP signaling pathways. Figure 5 Among them, 50 genes contributed the most to the enrichment score, which are the target gene set (Table 1).

[0120] Table 1. Target gene set

[0121]

[0122] Note: ES indicates the degree of enrichment of gene set members at both ends of the sort list; RM (Rank at MAX) is the position of the gene set's ES value in the sort list; LE (Leading-edge) is the gene member that contributes the most to the enrichment score.

[0123] Step 5: Use expression level correlation analysis to obtain the key gene set regulating the aggressive phenotype:

[0124] 1. Expression correlation analysis of the target gene set was performed using R language. The analysis method was Spearman correlation analysis, with Benjamini-Hochberg multiple correction. The significance level was P<0.05, and the absolute correlation threshold was 0.8. 2. The connectivity of each gene in the target gene set was calculated, and its absolute value was denoted as D. Connectivity refers to the sum of the correlations between a gene and other genes. Key genes in the network are identified through gene connectivity.

[0125] 3. Calculate the mean D value for all genes using the following formula:

[0126] DM = (D1 + D2 + ... + Dq) / q

[0127] Where DM is the mean of D, Dq is the gene connectivity, and q is the number of genes in the target gene set that have significant correlations.

[0128] 4. Use the igraph package to visualize the gene co-expression network and further analyze these genes with high connectivity;

[0129] 5. Using D≥DM as the screening criterion, obtain the key gene set regulating the aggressive phenotype;

[0130] 6. The results showed that 13 key genes were associated with differences in aggressive phenotypes, including NR2B, CaMKII, CREB, SERCA,

[0131] 5-HTR1, DAR1, etc., which are key gene sets (Table 2).

[0132] Table 2. Key gene set

[0133]

[0134] Step 6: Focus on receptor proteins with concentrated key genes as core regulatory proteins:

[0135] The results showed that the key gene cluster consisted of receptor proteins NMDAR, 5HTR1, and DAR1. Using these as core regulatory proteins, the key gene cluster formed an aggressively regulated pathway centered on NMDAR-CaMKII, cascading through cAMP, oxidative phosphorylation, PI3K-Akt, AMPK, and FoxO signaling pathways. Figure 6 ).

[0136] Step 7: Verify the modulatory effect of NMDAR on aggressive behavior through in vivo injection experiments:

[0137] 1. After the temporary rearing period, randomly select no fewer than 60 three-spined crabs in the intermolt stage with intact appendages. After starvation for 24 hours, determine their initial aggression through a mirror image experiment. Repeat the mirror image experiment 3 times for each crab, and the photography requirements are the same as in step 1;

[0138] 2. After at least 24 hours of recovery, the drugs were injected. The experimental group was injected with NMDAR activator (NMDA) and inhibitor (MK801), while the control group was injected with PBS. Each treatment group consisted of 20 swimming crabs, with 10 strong and 10 weak crabs.

[0139] 3. Use a pointed microsyringe to inject 100 μL into the base of the swimming foot. The injection volume is 0.1 mM and the injection depth is 4 mm. Ensure that the solution is injected into the pericardial cavity and completely released.

[0140] 4. The aggressiveness of the swimming crabs was measured and evaluated again 15 minutes after injection;

[0141] 5. Results showed that the trends in aggression in *Portunus trituberculatus* with different aggressive phenotypes differed after injection of NMDA activators and inhibitors. Compared with the control group injected with PBS, NMDA injection significantly reduced aggression in the highly aggressive group (P<0.05), while increasing aggression in the weakly aggressive group; similarly, MK801 injection reduced aggression in the highly aggressive group and increased aggression in the weakly aggressive group. Figure 7 );

[0142] 6. Immediately after the behavioral test, the swimming crabs were cryo-anesthetized, dissected on ice, and the thoracic ganglia were removed and placed in 2ml cryovials. One part was briefly frozen in liquid nitrogen and then stored at -80℃ for real-time quantitative PCR analysis; the other part was fixed with 4% paraformaldehyde for tissue immunofluorescence assay.

[0143] Step 8: The regulatory role of the NMDAR-CaMKII pathway on the aggressiveness of the swimming crab *Portunus trituberculatus* was verified by real-time quantitative PCR and tissue immunofluorescence.

[0144] Real-time quantitative PCR includes the following steps:

[0145] 1. Total RNA was extracted from the thoracic ganglion samples of swimming crabs using the Trizol method and analyzed by real-time quantitative PCR.

[0146] 2. RNA integrity was determined by 1% agarose gel electrophoresis, and the concentration and purity of total RNA were detected by a Nanodrop 2000 nucleic acid and protein analyzer, ensuring that the A260 / A280 ratio was between 1.8 and 2.0;

[0147] 3. Adopt III RT SuperMix kit was used to reverse transcribe and synthesize cDNA templates for testing;

[0148] 4. The relative expression levels of key genes were determined using the Taq Pro Universal SYBR qPCR Master Mix kit;

[0149] 5. Primers were designed using primer 5, and the amplification efficiency was verified to be between 1.8 and 2.2. β-actin was used as an internal reference gene. 6. Amplification conditions were: pre-denaturation at 95℃ for 30 s; cycling reaction at 95℃ for 10 s, 57℃ for 30 s, for 40 cycles; melting curve at 95℃.

[0150] 15s, 60℃, 60s, 95℃, 15s. Perform three technical replicates for each parallel sample. Quantitative fluorescence results are presented in 2... -△△Ct The method is used for analysis.

[0151] Tissue immunofluorescence assay includes the following steps:

[0152] 1. Fresh thoracic ganglion tissue samples were fixed with 4% paraformaldehyde for 24 hours, dehydrated with graded ethanol, cleared with xylene, embedded in paraffin, sectioned, and sectioned to a thickness of about 5 μm. The tissue was then mounted on a glass slide to prevent detachment and dried at 37°C overnight.

[0153] 2. Paraffin sections are first treated with environmentally friendly dewaxing solutions I, II, and III for 10 minutes each, and then rehydrated in a gradient of anhydrous ethanol.

[0154] 3. Next, the sections were microwaved in a citric acid retrieval solution at pH 6.0 for antigen retrieval using a microwave-safe program of 8 minutes on medium heat, 8 minutes off, and 7 minutes on medium-low heat, and then allowed to cool naturally.

[0155] 4. After washing the slides in PBS (pH 7.4), block them with 3% BSA for 30 minutes, then add the primary antibody diluted 1:50 and incubate overnight at 4°C.

[0156] 5. The next day, after washing the sections, add secondary antibody diluted 1:400 and incubate at room temperature for 50 minutes in the dark.

[0157] 6. DAPI staining solution is used for nuclear counterstaining. Incubate at room temperature in the dark for 10 minutes, then treat with autofluorescence quencher B solution for 5 minutes, and rinse with running water for 10 minutes to quench autofluorescence.

[0158] 7. Mount the slides using anti-fluorescence quenching mounting medium and observe and photograph them under a fluorescence microscope through the DAPI (blue) and 488 (green) channels. Figure 8 );

[0159] 8. The relative quantification of protein fluorescence intensity was performed using ImageJ software.

[0160] The above experimental results show that:

[0161] After NMDA injection, compared with the PBS group, the W group showed upregulation of the NR2B gene, a significant decrease in the average fluorescence intensity of NMDAR protein, and upregulation of the relative expression levels of aggression-related genes CaMKII, SERCA, PKA, COX2, ND1, CYTB, and FOXO, while downregulation of CREB and FOSLN. In the S group, the S group showed downregulation of the NR2B gene and a significant increase in the average fluorescence intensity of NMDAR protein. Figure 9 A) The relative expression levels of CaMKII, SERCA, PKA, CREB, FOXO, and FOSLN were downregulated, while the relative expression levels of genes related to the oxidative phosphorylation pathway, COX2, ATP6, ND1, and CYTB, were upregulated (P<0.05). Figure 9 ).

[0162] After injection of MK801, compared with the PBS group, the NR2B gene was upregulated in the W group, and the average fluorescence intensity of NMDAR protein was significantly reduced. Figure 9 A) All 11 key genes related to aggression were downregulated; in group S, NR2B gene was downregulated, the average fluorescence intensity of NMDAR protein was increased, CaMKII, SERCA, CREB, and FOSLN were significantly downregulated, PKA, COX2, ATP6, and CYTB were downregulated, and ND1 was upregulated (P<0.05). Figure 9 ).

[0163] The verification experiment confirmed that it is feasible to use receptor proteins with key gene concentrations as the core regulatory proteins in step 6.

[0164] This invention provides a validated method for screening crab aggression genes based on high-throughput sequencing. Through hierarchical screening, a set of key genes related to aggression regulation was constructed, key signaling pathways regulating crab aggression were identified, the correlation between aggression phenotypes and core genes was revealed, and molecular markers of crab aggression phenotypes were identified and validated. This method utilizes molecular biology techniques to achieve precise breeding at the gene level, providing theoretical and practical support for the application of molecular breeding technology in crustaceans and for future crab germplasm innovation and new variety breeding.

Claims

1. A method for screening crab aggressive genes based on high-throughput sequencing, comprising the following steps: Step 1: Assess crab aggression through mirror experiments and aggression assessment models, and select crabs with different aggression phenotypes; Its features also include the following steps: Step 2: Using transcriptomics analysis, identify differentially expressed genes in the thoracic ganglia of crabs with different aggressive phenotypes, with at least 3 biological replicates for each group; Step 3: Perform KEGG enrichment analysis on the differentially expressed genes and screen the top n differentially expressed signaling pathways to obtain a subset of aggressive differentially expressed genes. Step 3.1: Use R language to perform KEGG pathway enrichment analysis on the differentially expressed genes obtained in Step 2; Step 3.2: Apply Fisher's exact test to the above analysis results. When P When the value is less than 0.05, it indicates that the KEGG pathway function is significantly enriched. The number of significantly enriched KEGG pathways is denoted as Q1. Step 3.3, sort the significantly enriched KEGG pathways according to... P The values ​​are sorted from smallest to largest, and the number of genes in the final differential gene subset is denoted as Q2: Case 1: Q1 ≥ 20 If Q2 ≤ 500 when n = 20, then the genes enriched in the top 20 differential signaling pathways are selected as the differential gene subset. If Q2 > 500 when n = 20, then calculate the genes enriched in the first n-1 differential signaling pathways until Q2 ≤ 500; Case 2: Q1 < 20 If when n = Q1, Q2 ≤ 500, then the genes enriched in the first Q1 differential signaling pathways are selected as the differential gene subset. If Q2 > 500 when n = Q1, then calculate the genes enriched in the first n-1 differential signaling pathways until Q2 ≤ 500; Step 4: Using Gene Set Enrichment Analysis (GSEA), select genes that contribute significantly to the enrichment score from the differentially expressed gene subset as the target gene set. Step 4.1: If Q2 ≤ 50, then all genes in the differential gene subset are taken as the target gene set, and proceed directly to step 5. If Q2 > 50, then perform GSEA analysis on the differentially expressed gene subset using R language; Step 4.2, the screening criteria for GSEA are |NES| > 1 and FDR < 0.25, where NES is the standardized enrichment score and FDR is the score after multiple hypothesis testing correction. P The larger the |NES| value, the smaller the FDR value, and the more significant the enrichment. Step 4.3: For genes in the differential gene subset, use the Signal2Noise sorting algorithm based on the NES score. The size of the prior gene set is 15 ~ Q2. Step 4.4: Select the gene member combination that contributes the most to the enrichment score as the target gene set; Step 5: Use expression level correlation analysis to obtain the key gene set regulating the aggressive phenotype; Step 6: Focus on receptor proteins with a concentration of key genes as core regulatory proteins.

2. The method for screening crab aggressive genes based on high-throughput sequencing as described in claim 1, characterized in that: In step 1, the attack assessment is performed in the following manner: Step 1.1: Isolate and temporarily keep the crabs in an aquarium for two weeks; Step 1.2: Select healthy crabs with intact appendages and in the intermolting period for mirror imaging experiments. Repeat the imaging 3 times, with an interval of no less than 24 hours between each time. Step 1.3: Record and quantify the crab’s aggressive behavior using behavior analysis software, including the cumulative duration of stillness and the relative distance traveled—the ratio of the distance traveled to the width of the carapace. Step 1.4: Based on the crab aggression assessment model, calculate the aggression score to obtain crabs with different aggression phenotypes.

3. The method for screening crab aggressive genes based on high-throughput sequencing as described in claim 2, characterized in that: In step 1, the mirror experiment is conducted in the following manner: An observation tank with an opaque white inner wall was selected. The camera was fixed directly above the tank. Two mirrors with a width equal to the inner diameter of the tank were placed vertically in the middle of the tank with their backs touching, so that the two sides could be mirrored simultaneously. A partition was inserted into each side of the mirror in the tank. Two crabs were placed in the partitions on one side of the partitions, separating them from the mirrors. After acclimatizing for 10 minutes, the partitions were gently removed, exposing the crabs to the mirrors. Filming was conducted for 20 minutes. After each set of filming, the seawater in the observation tank was replaced and the tank was cleaned. Filming was conducted under completely dark conditions, and the seawater temperature and salinity were kept consistent with the temporary holding conditions during the experiment.

4. The method for screening crab aggressive genes based on high-throughput sequencing as described in claim 2, characterized in that: In step 1, the aggressive behavior of the crabs was recorded and quantified using behavioral analysis software as follows: The detection sampling rate is set to 8 frames / second, and the detection method is either dynamic silhouette or grayscale gradient. Lost frame correction and smooth trajectory settings are enabled. The trajectory smoothing configuration is as follows: a. Smoothing is performed based on 10 samples before and after each sampling point; b. If the moving distance is less than 3 cm, the sampling point will still be set to the previous position; c. If the maximum distance moved is greater than 20 cm, the sampling point will be set as missing.

5. The method for screening crab aggressive genes based on high-throughput sequencing as described in claim 2, characterized in that: In step 1, the aggression assessment model is: Y = 0.023X1 - 0.001X2 - 0.002 Where Y is the aggression score, X1 is the relative movement distance, and X2 is the cumulative duration of stationary movement.

6. The method for screening crab aggressive genes based on high-throughput sequencing as described in claim 1, characterized in that: Step 2 includes collecting raw sequencing data and transcriptomics analysis: The raw sequencing data collected includes: Step 2.1: Randomly select at least 2 crabs of each aggressive phenotype from Step 1 for pooling, with at least 3 biological replicates in each group; Step 2.2: Total RNA was extracted using the TRIzol method. RNA integrity was determined by 1% agarose gel electrophoresis. The concentration and purity of total RNA were detected using a Nanodrop 2000 nucleic acid and protein analyzer to ensure that the A260 / A280 ratio was between 1.8 and 2.

0. Step 2.3: After extracting 1 μg of total RNA from the sample, use Oligo magnetic beads to pair the bases with the poly A at the 3' end of the mRNA. After purifying the mRNA, add fragmentation buffer to randomly break the mRNA into small fragments of about 300 bp. Step 2.4: Using mRNA as a template, reverse synthesize one-stranded cDNA, then synthesize two-stranded cDNA, add end repair reagent EndRepair Mix to complete the blunt ends, then add a base A to the 3' end, then purify and sort the product after ligating the adapter, use the sorted product for PCR amplification, and purify to obtain the final library; Step 2.5: The constructed library was sequenced using Illumina HiSeq Xten / NovaSeq6000 to obtain raw sequencing data; Transcriptomics analysis includes: Step 2.6: Use SeqPrep and Sickle software to filter the raw sequencing data to obtain high-quality sequencing data (clean reads); Step 2.7: Use HiSat2 and TopHat2 software to align the quality-controlled clean reads with the reference genome to obtain mapping data, and evaluate the quality of the alignment results. Step 2.8: The gene expression level is calculated by the number of clean reads located in the genomic region (read counts). The gene expression level is quantitatively analyzed using RSEM software. The quantitative expression results are expressed in FPKM (female fraction of kilobases). Step 2.9: Differential gene expression levels were analyzed using DESeq2 software to identify differentially expressed genes, and then genes with significant differential expression were screened. The screening criteria were: FDR < 0.05 & |log2FC| ≥ 2, where FDR is the result after multiple hypothesis testing correction. P value.

7. The method for screening crab aggressive genes based on high-throughput sequencing as described in claim 1, characterized in that: In step 5, obtaining the key gene set regulating the aggressive phenotype using expression level correlation analysis includes the following steps: Step 5.1: Perform expression correlation analysis on the target gene set using R language to generate a correlation matrix, where each element in the matrix represents the correlation between two genes, and construct a co-expression matrix; Step 5.2: Spearman correlation analysis was performed on the co-expression matrix, and the results were corrected using Benjamini-Hochberg multiple tests at a significance level of [missing value]. P <0.05, the absolute threshold for correlation is 0.8; Step 5.3: Calculate the connectivity of each gene in the target gene set and denote its absolute value as D; connectivity refers to the sum of the correlations between a gene and other genes. The connectivity of each gene is used to identify key genes in the network. Step 5.4, calculate the mean D value for all genes using the following formula: DM = (D1+D2+…+Dq) / q Where DM is the mean of D, D1,…,Dq are gene connectivity, and q is the number of genes in the target gene set that have significant correlations. Step 5.5: Visualize the gene co-expression matrix using the R package; Step 5.6: Using D ≥ DM as the screening criterion, obtain the key gene set regulating the aggressive phenotype.

8. The application of any one of the high-throughput sequencing-based crab aggressive gene screening methods described in claims 1-7 in crab germplasm innovation and new variety breeding.

Images