GWAS and multi-omics data-based larch pine moth-resistant gene mining method

Abstract

The invention provides a larch pine caterpillar resistance gene mining method based on GWAS and multi-omics data, and belongs to the technical field of gene mining. The method comprises the following steps: S1, carrying out insect-resistant phenotype identification on larch individuals, and carrying out variation identification on the larch individuals based on a liquid-phase gene chip sequencing technology; s2, using a mixed linear model to perform correlation analysis by taking a plot position and phenotype observation time as an interaction environment, taking population density at different observation time as a covariable and taking a larch anti-deciduous pine caterpillar index as a phenotype, and identifying a candidate QTL region; and S3, carrying out transcriptome sequencing on resistance extreme individuals, carrying out variation recognition, comparing read segments to the candidate QTL region, and manually analyzing and annotating a gene structure in the QTL region by using IGV. According to the method disclosed by the invention, low-cost variation recognition, environment noise weakening whole genome association analysis and gene structure accurate annotation on complex genome species are realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of gene mining, and particularly relates to a larch anti-Dendrolimus punctatus Wang et Cramb gene mining method based on GWAS and multi-omics data. BACKGROUND

[0002] Larch (Larix gmelinii) is an important fast-growing timber species with economic and ecological protection value. The economic value of larch mainly lies in its fast-growing and high-yield characteristics, and its wood quality is excellent, which is an important raw material for the construction, furniture and papermaking industries, bringing considerable economic benefits to forestry production. In terms of ecological value, larch forests play a key role in water conservation, soil conservation, wind prevention and sand fixation. Dendrolimus punctatus Wang et Cramb is one of the main forest pests that harm larch and is a third-class harmful forest pest. Its outbreak can cause large-scale death of larch and cause serious economic and ecological losses. In addition, it can also cause tree branches to die, increasing the risk of forest fires. Therefore, breeding larch varieties resistant to Dendrolimus punctatus Wang et Cramb is one of the important goals of larch breeding research. Larix spp. Dendrolimus superans Accurate identification of quantitative trait loci (QTL) related to insect resistance is a core goal of larch insect resistance breeding. Traditional QTL mapping relies on parent-derived genetic populations, but due to the large genome of trees and narrow genetic background, the results have poor reproducibility and limited application in different populations. Genome-wide association analysis (GWAS) can break through the limitations of parent-derived populations, but it faces three technical bottlenecks in larch: (1) The larch genome is large, and whole-genome resequencing is costly and difficult to scale; (2) Trees have a long growth cycle and high environmental heterogeneity, and environmental noise masks the true genetic signal; (3) The larch reference genome is not fully annotated, and it is difficult to link the associated loci to functional genes. The above bottlenecks result in low reliability and difficulty in locating key genes for QTL identified by existing GWAS in larch.

[0003] In order to reduce the cost of practical GWAS in tree genetic breeding and improve the reliability and applicability of the results, multi-omics analysis has become an effective solution. Through liquid chip technology, a large number of variations in the larch genome are uniformly and randomly identified to provide a basis for GWAS analysis. By using a mixed linear model to eliminate false positives caused by genetic relationships and to assess the interaction effects between genes and environments, the model reliability is improved. By using transcriptome sequencing data to accurately depict the gene structure around the QTL identified by GWAS, conditions are provided for downstream applications.

[0004] In order to reduce the cost of practical GWAS in tree genetic breeding and improve the reliability and applicability of the results, multi-omics analysis has become an effective solution. Through liquid chip technology, a large number of variations in the larch genome are uniformly and randomly identified to provide a basis for GWAS analysis. By using a mixed linear model to eliminate false positives caused by genetic relationships and to assess the interaction effects between genes and environments, the model reliability is improved. By using transcriptome sequencing data to accurately depict the gene structure around the QTL identified by GWAS, conditions are provided for downstream applications. SUMMARY

[0005] ​The present application aims to provide a larch anti-moth gene mining method based on GWAS and multi-omics data, which combines liquid gene chip, mixed linear model and transcriptome-based gene structure analysis, breaks through the limitations of single-omics whole genome association analysis in forest tree genetic breeding, and realizes low-cost variation identification, environment noise reduction of whole genome association analysis and accurate gene structure annotation for complex genome species.

[0006] In order to achieve the above-mentioned application purposes, the present application provides the following technical solutions. The present application provides a larch anti-moth gene mining method based on GWAS and multi-omics data, comprising the following steps: S1, identifying the insect resistance phenotype of larch individuals, calculating the larch anti-larch moth index, and identifying the variation of larch individuals based on liquid gene chip sequencing technology; S2, using mixed linear model, taking plot location and phenotype observation time as interactive environment, taking insect density at different observation times as covariate, and taking larch anti-larch moth index as phenotype for association analysis to identify candidate QTL region; S3, performing transcriptome sequencing on the extreme individuals, identifying the variation, aligning the reads to the candidate QTL region, and using IGV to artificially analyze and annotate the gene structure in the QTL region.

[0007] Preferably, the sample number of the larch individuals in step S1 is 400-600, and at least comes from 250 different genetic backgrounds.

[0008] Preferably, the method for identifying the insect resistance phenotype in step S1 is: placing 3-7 larch moth larvae and tender larch branches in a net bag for 8-12 days, recording the initial and final body weight of the larch moth larvae, the loss of larch needles, the weight of feces produced during the bagging process, and calculating the larch anti-larch moth index. The calculation method of the larch anti-larch moth index is:

[0009] In the formula, PL represents the loss of larch needles, unit: root; WC represents the weight change of larch moth, unit: gram; FW represents the weight of feces produced by larch moth during the bagging process, unit: gram.

[0010] Preferably, step S1 further comprises a step of quality control of the sequencing data and the variation identification result, and the quality control standard of the variation identification result is: site missing rate < 5%, minor allele frequency > 5%, Hardy-Weinberg equilibrium test P value < 1e-6, and sample missing rate < 5%.

[0011] Preferably, the QTL region in step S2 is a QTL region with a significance P value < 1e-5 / SNP number and a linkage disequilibrium coefficient r² > 0.2 under a multi-environment model.

[0012] Preferably, the method used in the variation identification in step S3 is a k-mer method.

[0013] Preferably, the method for annotating the gene structure in the QTL region in step S3 is: based on the alignment information of the transcriptome sequencing data in the candidate QTL region, the transcription start and end sites and the alternative splicing information are identified manually, the specific CDS sequence is spliced, and sequence alignment is performed in the NCBI and Uniprot databases to explore the function of the gene.

[0014] The application provides an application of the method in breeding of larches resistant to pine caterpillars, and uses the candidate QTL region and the gene structure for SNP marker development in whole genome selection or CRISPR gene editing target point design.

[0015] Compared with the prior art, the application has the following beneficial effects: 1. The multi-omics analysis method for resistance to larch caterpillars based on whole genome association analysis of a larch breeding population provided by the application realizes low-cost variation identification based on gene chip sequencing technology, eliminates the influence of environmental factors on complex traits by using a mixed linear model, accurately captures candidate QTL, and identifies related gene structures in the QTL by using transcriptome sequencing data. Reliable target sites are provided for breeding of larch varieties resistant to larch caterpillars.

[0016] 2. The multi-omics analysis method for resistance to larch caterpillars based on whole genome association analysis of a larch breeding population provided by the application integrates the advantages of gene chips, mixed linear models and transcriptome sequencing, obtains candidate QTLs significantly related to insect resistance traits, and analyzes gene structures around the QTLs. Compared with the QTL results positioned in trees by using GWAS, the identified QTL region has linkage characteristics, and the gene structure in the region is clear, which provides convenience for downstream research. It is expected to further improve the efficiency of selection breeding in species with long life cycles, inconsistent growth environments and complex genomic backgrounds by using GWAS, and has a wide application prospect. BRIEF DESCRIPTION OF DRAWINGS

[0017] In order to make the technical solutions in the embodiments of the present application or the prior art clearer, the accompanying drawings needed in the embodiments or prior art description will be briefly introduced below. Obviously, the accompanying drawings in the following description only aim at the embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative effort on the basis of the provided drawings.

[0018] Figure 1 Panorama of the results of the whole genome association analysis in Example 1; Figure 2 Panorama of the results of the gene structure analysis in Example 1. DETAILED DESCRIPTION

[0019] The present application provides a larch anti-Dendrolimus punctatus Walker gene mining method based on GWAS and multi-omics data, comprising the following steps: S1, identifying the insect resistance phenotype of larch individuals, calculating the larch anti-Dendrolimus punctatus Walker index, and identifying the variation of larch individuals based on liquid gene chip sequencing technology; S2, using a mixed linear model, taking the plot location and the phenotype observation time as the interactive environment, taking the insect density at different observation times as the covariate, and taking the larch anti-Dendrolimus punctatus Walker index as the phenotype for association analysis to identify candidate QTL regions; S3, performing transcriptome sequencing on the resistance extreme individuals, identifying the variation, aligning the reads to the candidate QTL region, and using IGV to artificially analyze and annotate the gene structure in the QTL region.

[0020] In the present application, the purpose of step S1 is to obtain as many sample genotypes and phenotype data as possible at low cost in the forest population, the purpose of step S2 is to eliminate the false positive results of association analysis caused by genetic relationship and population structure and identify candidate QTL regions, and the purpose of step S3 is to focus on the candidate QTL region using transcriptome sequencing data, to improve the resolution of variation identification on the basis of analyzing the gene structure, and to further verify the correspondence between different genotypes and phenotypes in the resistance extreme individuals.

[0021] In the present application, the number of samples of the larch individuals in step S1 is 400-600, preferably 450-550, and further preferably 511; at least from 250 different genetic backgrounds, preferably at least from 280 different genetic backgrounds, and further preferably at least from 293 different genetic backgrounds.

[0022] In this invention, the method for identifying the insect resistance phenotype in step S1 is as follows: 3-7 larch caterpillars and young larch branches are placed in a net bag and bagged for 8-12 days. The initial and final weight of the caterpillars, the amount of pine needle loss, and the weight of insect excrement produced during the bagging process are recorded. The larch resistance index against larch caterpillars is calculated. The method for calculating the larch resistance index against larch caterpillars is as follows:

[0023] In the formula, PL This indicates the amount of larch needle loss, expressed in units of needles. WC This indicates the change in body weight of the larch caterpillar, in grams. FW This indicates the weight of the excrement produced during the bagging process of larch caterpillars, in grams.

[0024] In this invention, when identifying the insect-resistant phenotype, it is preferable to place 4-6 larch caterpillars and young larch branches in a net bag and cover it for 9-11 days, and more preferably to place 5 larch caterpillars and young larch branches in a net bag and cover it for 10 days.

[0025] In this invention, during the identification of the insect-resistant phenotype, the spatial location (forest plot) of the larch individual, the phenotype observation time, and the insect population density after the pine caterpillar infestation are also recorded.

[0026] In this invention, step S1 further includes a step of quality control of sequencing data and variant identification results. The quality control criteria for variant identification results are: site deletion rate < 5%, minor allele frequency > 5%, Hardy-Weinberg equilibrium test P value < 1e-6, and sample deletion rate < 5%.

[0027] In this invention, the QTL region mentioned in step S2 is a QTL region with a significance P-value < 1e-5 / SNP number (based on Bonferroni correction) and a linkage disequilibrium coefficient r² > 0.2 under the multi-environment model.

[0028] In this invention, the method used for mutation identification in step S3 is the k-mer method.

[0029] In this invention, the method for annotating the gene structure within the QTL region in step S3 is as follows: based on the alignment information of transcriptome sequencing data in the candidate QTL region, the transcription start and stop sites and alternative splicing information are manually identified, and after splicing the specific CDS sequence, the sequence is aligned in the NCBI and Uniprot databases to explore the function of the gene.

[0030] This invention provides an application of the method described in larch breeding for resistance to pine caterpillars, using candidate QTL regions and gene structures for SNP marker development in genome-wide selection, or for CRISPR gene editing target design.

[0031] The technical solutions provided by the present invention will be described in detail below with reference to the embodiments, but they should not be construed as limiting the scope of protection of the present invention.

[0032] Example 1

[0033] A method for mining larch resistance genes against pine caterpillars based on GWAS and multi-omics data, the steps of which are as follows: 1. Observation of insect resistance phenotype Using 511 larch individuals belonging to 293 genetic backgrounds as the research subjects, 80 cm long young larch branches were taken, retaining lateral branches with a distance of less than 30 cm between the base and the terminal bud, and removing other lateral branches. Five larch caterpillars were taken, and their weight m1 was recorded. They were then placed in a mesh bag with the lateral branches facing inwards. The bag was tied tightly to isolate the larvae and lateral branches from the outside environment. The lateral branches covered with mesh bags were inserted into a bucket of water, with the water level no higher than 20 cm above the mesh bag. The water was changed every 2 days. After 10 days, the bag was removed, and the weight m2 of the pine caterpillars, the amount of pine needle loss n, and the weight w of the frass were recorded. At the same time, the spatial location (forest plot), the time of phenotypic observation, and the population density after pine caterpillar infestation were also recorded.

[0034] 2. Gene chip mutation identification

[0035] DNA was extracted using the CTAB method (concentration >300ng / μL), and variants were identified at 60,000 loci in larch using liquid-phase microarray sequencing technology. The vcf format files were quality controlled and converted using PLINK2. Quality control included deletion rate control (less than 5%), allele frequency (greater than 5%), HW balance test (less than 1e-6), and high deletion samples (less than 5%).

[0036] 3. Genome-wide association analysis

[0037] Using 3VmrMLM with cell location and phenotypic observation time as the interaction environment, insect population density at different observation times as the covariate, and larch resistance to larch caterpillar index as the phenotypic cricket, association analysis was performed to identify candidate QTL regions. The QTL regions are those with a significance P-value < 1e-5 / SNP number (based on Bonferroni correction) and a linkage disequilibrium coefficient r² > 0.2 under the multi-environment model.

[0038] The method for calculating the larch resistance index against larch caterpillars is shown in Equation 1.

[0039]

[0040] Formula 1

[0041] In the formula, PL Loss of larch needles (roots); WC : 10-day weight change (grams) of larch caterpillar, WC = m2-m1; FW The weight (in grams) of frass produced by the larch caterpillar over 10 days.

[0042] 4. Gene structure analysis

[0043] Transcriptome sequencing was performed on nine resistant individuals (three biological replicates). For quality-controlled transcriptome reads, a k-mer-based method was used to identify variants in candidate QTL regions (k=27). Splash2 was used to analyze the variation of target kmers in the anchor-target model. A generalized linear model was used to statistically analyze whether there were significant differences in the frequency of target kmers corresponding to anchor kmers between individuals in the resistant and susceptible groups. Regions containing significantly different anchor kmers were selected for gene structure analysis. Transcriptome reads aligned to candidate QTL regions were extracted and visualized using IGV. Candidate related gene structures (coverage >20×) were manually assembled based on the transcription start and stop sites and alternative splicing information displayed by the reads. Gene CDS sequences were extracted and sequence aligned in the NCBI and Uniprot databases to analyze gene function.

[0044] 5. Results

[0045] A panoramic view of the genome-wide association analysis results is as follows: Figure 1 The image shows a scatter plot of a Manhattan plot; points with high vertical axes represent genomic regions containing variant events related to the corresponding phenotype (insect resistance), while stacked scatter points indicate that multiple phenotype-associated, interconnected loci, i.e., QTL regions, have been identified in this region. Figure 1 As can be seen from the above association analysis, a QTL hotspot consisting of 10 mutually linked SNPs was identified in larch with a complex genomic background. Some of the association analysis results are shown in Table 1.

[0046] Gene structure analysis results as follows Figure 2 As shown, the alignment information of transcriptome reads in QTL hotspot regions reveals the transcription start sites, exon and intron locations, and transcription termination sites in these regions, thus elucidating the gene structure.

[0047] As can be seen from the above, this invention can accurately and reliably identify QTLs associated with complex traits and analyze the related gene structures in breeding populations with long life histories, inconsistent growth environments, and complex genomic backgrounds, providing a reference for marker-assisted selection of larch for resistance to larch caterpillars. Furthermore, the method of this invention can be extended to GWAS analysis of complex traits in different forest trees, showing broad application prospects.

[0048] Table 1. Partial Association Analysis Results

[0049] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for mining larch resistance to pine caterpillar genes based on GWAS and multi-omics data, characterized in that, Includes the following steps: S1. Identify insect-resistant phenotypes in larch individuals, calculate the larch caterpillar resistance index, and identify variations in larch individuals based on liquid-phase gene chip sequencing technology. S2. Using a mixed linear model with plot location and phenotypic observation time as the interaction environment, insect population density at different observation times as the covariate, and larch resistance to larch caterpillar index as the phenotypic cricket index, association analysis was performed to identify candidate QTL regions. S3. Transcriptome sequencing was performed on antagonistic extremist individuals to identify variants. Reads were aligned to the candidate QTL regions, and the gene structures within the QTL regions were manually analyzed and annotated using IGV.

2. The method as described in claim 1, characterized in that, The number of larch individuals in step S1 is 400 to 600, and they come from at least 250 different genetic backgrounds.

3. The method as described in claim 1, characterized in that, The method for identifying the insect resistance phenotype in step S1 is as follows: Place 3-7 larch caterpillars and young larch branches in a net bag and bag it for 8-12 days. Record the initial and final weight of the caterpillars, the amount of pine needle loss, and the weight of insect excrement produced during the bagging process. Calculate the larch caterpillar resistance index. The method for calculating the larch resistance index against larch caterpillars is as follows: In the formula, PL This indicates the amount of larch needle loss, expressed in units of needles. WC This indicates the change in body weight of the larch caterpillar, in grams. FW This indicates the weight of the excrement produced during the bagging process of larch caterpillars, in grams.

4. The method as described in claim 1, characterized in that, Step S1 also includes a quality control step for sequencing data and variant identification results. The quality control criteria for variant identification results are: site deletion rate < 5%, minor allele frequency > 5%, Hardy-Weinberg equilibrium test P value < 1e-6, and sample deletion rate < 5%.

5. The method as described in claim 1, characterized in that, The QTL region mentioned in step S2 is a QTL region with a significance P-value < 1e-5 / SNP number and a linkage disequilibrium coefficient r² > 0.2 under the multi-environment model.

6. The method as described in claim 1, characterized in that, The method used for mutation identification in step S3 is the k-mer method.

7. The method as described in claim 1, characterized in that, The method for annotating gene structures within QTL regions in step S3 is as follows: based on the alignment information of transcriptome sequencing data in candidate QTL regions, transcription start and stop sites and alternative splicing information are manually identified, specific CDS sequences are spliced ​​together, and sequence alignment is performed in the NCBI and Uniprot databases to explore gene function.

8. The application of the method according to any one of claims 1 to 7 in the breeding of larch resistant to pine caterpillars, characterized in that, Candidate QTL regions and gene structures can be used for SNP marker development in genome-wide selection or CRISPR gene editing target design.

Images