Identification of candidate genes associating with soybean cyst nematode in soybean (Glycine max L.) using BSA-seq

Plant materials and phenotype evaluation

A F₂ population containing 600 plants was constructed with a cross between Chidou4 (susceptible parental, FI <10%) and H-10 (resistant parental, FI >60%). Chidou4 and H-10 were cultured in Chifeng Academy of Agricultural and Animal Husbandry Sciences, and Shenyang Agricultural University, respectively. For nine differential cultivars, Peking, Pickett, PI90763, PI548316, PI89772, PI209332, PI437654, PI88788, and Lee68, their seeds were provided by Northeast Agricultural University, Key Laboratory of Soybean Biology, Ministry of Education.

H-10, Chidou4, and 600 differential F₂ plant were evaluated for their SCN (race 3) resistance within the climate room at Chifeng Academy of Agricultural and Animal Husbandry Sciences. SCN (race 3)-infected soil was added into plastic cups (Ø6 cm × h12 cm), separately. Lee68, Jiyu86, and Liaodou15 were seeded into the SCN-infested soil, cysts were harvested in roots at 30 days later with the 710–250 µm sieve tower from its 250 µm sieve. After washing, cysts were broken with the rubber stopper to collect eggs onto the sieve stack including 250, 75 and 25 µm sieves. Then, we rinsed the mixed sample in 25 µm sieve in the 50 mL plastic conical tube. After adding 50% sucrose solution, eggs were subjected to stirring and 3-min centrifugation at 2,000 rpm, while those in the supernatant or middle layer were harvested onto the 25 µm sieve.

Transplantation of H-10 and Chidou4 was completed, and three replicates, whereas that of 600 differential F₂ plants was conducted with no replicate being set. One plant was grown within the plastic cup (Ø6 cm × h12 cm) in one replicate. At five days following transplantation, seedlings were inoculated using 2,000 eggs/cup. Plant growth conditions were 28–24(L:D), 70–80% relative humidity, 16h:8 h (L:D) photoperiod, and daily watering. At 30 days post-inoculation, the nematode cysts in roots of every replicate were harvested to calculate cyst number using the image analysis software (Wei et al., 2022).

In this study, female index (FI) was determined and the phenotype data used to analyze QTL, below: FI (%) = (mean cyst number on every line/mean cyst number on Lee68) × 100. Lines were classified into resistant (FI <10), moderately resistant (10 ≤ FI <30), moderately susceptible (30 ≤ FI <60), or susceptible (FI ≥ 60) for classifying SCN response.

DNA preparation

When conducting disease resistance identification on plants, individual plants were labeled, leaf numbers were taken and stored in liquid nitrogen. Based on the disease resistance results, the corresponding top young leaves of the plants were taken. The genomic DNA of parent materials H-10, Chidou4, and offspring leaves was extracted using the improved cetyltrimethylammonium bromide (CTAB) method (Zeng et al., 2021b). The resultant solution was then diluted to the same concentration. Then the DNA of 30 extreme disease resistant families were mixed with 30 extreme disease susceptible families equally. Two parent pools and two extreme phenotype mixed pools were constructed for whole genome resequencing.

Library construction and BSA sequencing

Sample library construction and sequencing were completed by Gene Pioneer Biotechnology company ji (Gene Pioneer Biotechnology Co. Ltd., Nanjing, China). The qualified genomic DNA samples were randomly broken into specific sized fragments using Covaris ultrasound, and the fragments were slected using the Agency AMPure XP Medium Kit to concentrate the sample bands around 200–300 bp. The Qubit dsDNA HS AssayKit (500 assays) was used to detect the purified DNA sample size. End repair on the interrupted DNA fragment was performed, and the A base at the 3′ end as well as the specific sequencing primers to the DNA fragment were connected. The connecting products were purified, the appropriate sized DNA fragments were selected, and the DNA fragments on a cBot instrument were amplified to prepare the library. Quality testing was performed on the constructed library, and the qualified library on the machine with an average sequencing depth of 30× for parents and average sequencing depth of offspring 60× were sequenced. The reference genome was Wm82.a2.v1 version of soybean genome (Song et al., 2016). SNP detection was achieved using GATK (McKenna et al., 2010) tool, and mutation annotation and prediction of mutation impact were conducted with SnpEff software (Cingolani et al., 2012).

Sequencing data analysis

The raw reads obtained from sequencing were filtered, and BWA software was used to compare the filtered sequence numbers with soybean reference genome (https://www.ncbi.nlm.nih.gov/) sequences (Li & Durbin, 2010). The results were sorted, the quality was filtered, and duplicates were labeled for subsequent mutation type detection. GATK software was used to search for SNP and Index sites between reference genome and sample. SNP is used for variable filtering parameters in GATK (settings: - filterExpression “QD < 4.0 ∥FS > 60.0 ∥ MQ < 40.0”, - G_filter “GQ < 20”, - cluster Windows Size 4). InDel filters through variable filtering parameters (setting: - Filter expression “QD < 4.0 ∥FS > 200.0 ∥ Read PosRankSum < − 20.0 || Inbreeding coefficient < − 0.8”). ANNOVAR (Wang, Li & Hakonarson, 2010) is an efficient software tool for annotating SNPs or InDel based on reference genome GFF3 files. Homozygous SNPs/InDel between two parents should be extracted from the VCF file of SNP/InDel. The reading depth information of homozygous SNPs/InDel above the offspring library can be obtained to calculate the SNP/InDel index.

The calculation of ΔSNP-index was as follows: SNP index (aa) = Maa/(Paa + Maa), SNP index. In the formula, M and P represent parent plants; Maa (Paa) represents the depth of the aa group derived from M (P); And Mab (Pab) refers to the depth of the ab group derived from M (P). We use the genotype of one parent as a reference and use statistical readings from that parent’s genotype or other genotypes in the offspring pool. Then, we calculated the ratio of different read counts to the total number, i.e., the SNP/InDel index of the base site. We filtered out points in both pools where the SNP/InDel index is less than 0.3. Sliding window method was used to present the SNP/InDel index of the entire genome. The average value of all SNP/InDel indexes in each window is the SNP/InDel index of that window. Usually, we use a 1 Mb window size and a 10 kb step size as default settings. The difference in SNP/InDel indices between two pools is calculated as the delta SNP/InDels index. Breakdancer (McKenna et al., 2010) software was used to detect SV sites between reference genome and sample. Control FREEC (Boeva et al., 2012) software was used to search for CNV regions between reference genome and sample. In order to obtain more accurate mutation information, further filtering was required for each mutation site, and the mutation site annotation software Anovar was used to annotate the location region of the mutation site on the genome.

Candidate gene localization and analysis

Based on parents and disease-resistant/susceptible single plant mixed pools, the SNP allele frequency of the mixed pool as well as the difference in allele frequency of the same site between two mixed pools were calculated. The positioning interval was obtained. The SNP index of both offspring in the extreme mixing pool, which is less than 0.30 or greater than 0.75 was filter out. An SNP that was not detected in one of the pools was subtracted from the SNP indices of the two extreme mixing pools to obtain a Delta SNP index. Delta SNP index within sliding window was calculated using a sliding window method, and the mean of diverse SNP sites within the sliding window was obtained. The sliding window size was 1.00 Mb in steps of 0.05 Mb. Genes within candidate regions were annotated with the NCBI database (https://www.ncbi.nlm.nih.gov) and clustering analysis was performed using Mercator software (https://www.plabipd.de/mercator_main.html) to classify and calculate the metabolic pathways involved in genes.

qRT-PCR for candidate genes

The disease-resistant variety H-10 and the susceptible variety Chidou4 was planted in a flowerpot containing SCN Race 3 diseased soil (about 50 cysts/100g soil). After 14 days of emergence, equal samples were taken from the same position at the root of the plant, with a sampling amount of 5 g. This was repeated three times and the samples were quickly placed in liquid nitrogen for freezing. The screened genes were selected, primers were designed, RNA was extracted, and reverse transcription was performed. With ACTIN7 as the internal reference gene, relative quantitative analysis of the data was performed using the 2^−ΔΔCt method. The two-step fluorescence quantitative PCR amplification procedure was conducted as follows: for the amplification curve, 5 min under 95 °C, with one cycle; 10 s under 95 °C, 30 s under 60 °C, 40 cycles. For the melting curve, 15 s under 95 °C, 60 s under 60 °C, and 15 s under 95 °C, and the signal was continuously detected.

Phenotypic variation

The plants were cultured as described in the section of materials and methods. There was a significant differences in FI between the parents (Figs. 1A–1B). The average female index of male parent H-10 was 0.61%, making it a highly resistant material; the female index of Chidou4 is 71.28%, indicating high sensitivity. In the F₂ population of the hybrid offspring of H-10 and Chidou4 variation range was 0–313.13%, with an average variation rate of 68.02%, between the resistant and susceptible parents (Fig. 1C). The statistical analysis of F₂ population resistance level showed a normal distribution. These juvenile nematodes infect the soybean plant’s vascular tissue. Female nematodes can be observed at the plant’s roots and some of the symptoms of this disease include severe growth retardation, stunting, and a yellowish appearance (Fig. 1D).

Sequencing quality evaluation

Two DNA hybridization libraries were constructed from the F₂ population established from disease resistant parents (H-10) and susceptible parents (Chidou4), as well as 30 disease resistant and susceptible materials. There is a significant difference between the resistance pool (female index <10%) and the susceptibility pool (female index >60%). The entire genome of four samples from a mixed library of parents and offspring was re sequenced using high-throughput sequencing technology, resulting in a total of 112.47 GB of raw sequencing data. After joint filtering, low-quality filtering, and N-filtering using SOAPnuke software, high-quality sequence number data of 110.88Gb was finally obtained. The GC content ranges from 33.95% to 35.32%, indicating high quality sequencing data. The average mass value of the vast majority of base sequences is above 30 (Q30 ≥ 92.05%) (Table 1). This suggests an enough data volume for the four samples, reliable sequencing data, qualified sequencing quality, and normal GC distribution, so data can be analyzed in the next step. Following quality control with the soybean reference genome, BWA software (https://bio-bwa.sourceforge.net/) was used to compare filtered sequence numbers from the four samples and Qualimap software was used to compare and statistically analyze the results, The results showed that the genome size was 978,495,272 bp, the effective genome size was 954,832,134 bp (the reference sequence does not contain N). Meanwhile, it was found that the GC content of the reference genome was over 95%, with a comparison rate between 99.28% and 99.41% for the four samples, and with an average sequencing depth of 16.97% (Table 2). The high data comparison rate is conducive to the next step of SNP detection. After sequencing sequence numbers were compared against reference genome, the coverage and sequencing depth distribution of different chromosome regions was calculated on the reference genome. 1× coverage percentage exceeds 95.51%, 5× coverage percentage exceeding 93.02%, 10× coverage percentage exceeds 81.43% (Table 2). By comparing the results with Qualimap software, it was found that the similarity between each sample and the reference gene met the requirements of resequencing analysis. The sequencing data had good coverage depth and coverage, and the reference genome was evenly covered with good randomness, which was conducive to SNP filtering and screening.

SNP detection and annotation

The analysis results showed that a total of 1,152,230 SNPs were obtained from the four samples, with the least SNPs in the resistant offspring library and the most SNPs in H-10 (resistant parents). SNP loci are distributed on 20 chromosomes of soybean. A total of 128,911 non synonymous SNPs were obtained between the disease resistant material H-10 and the susceptible material Chidou4, and a total of 101,061 non synonymous single nucleotides were obtained from the mixed pool. The statistical data of SNP location and coding information indicate that the number and proportion of SNPs of the same type (between parents and offspring) in the same type of material are usually equal (Table 3). SNP loci are distributed in 15 regions, with the majority located in introns, gene intervals, upstream and downstream of genes, while the rest are distributed in variable splicing sites, synonymous mutations, non synonymous mutations, 3′-untranslated regions, 5′-untranslated regions, and unknown regions. The SNP ratio of homozygous mutations is significantly higher than that of heterozygous mutations, with a SNP ratio of ≥ 99.75% for homozygous mutations.

Screening and functional analysis of candidate genes

The detected SNPs were analyzed using the SNP index method, the SNP indices of two subsequent libraries were calculated the differences in SNP indices of different chromosome segments in the disease resistant and susceptible libraries were compared, and genetic testing was conducted based on G statistics. When both the p-value and the false discovery rate (FDR) are less than 0.001, it was found that the genes controlling disease resistance are mainly concentrated in the 38,524,128–39,849,988 bp fragment (1.33 Mb in length, including 111 genes) on chromosome 2. A total of 741 genes were predicted in the 27,821,012–29,612,574 bp region of chromosome 3 (1.79 Mb in length, containing 92 genes), the 308–348,214 bp region of chromosome 10 (0.35 Mb in length, containing 34 genes), and the 53,867,581–58,017,852 bp region of chromosome 18 (4.15 Mb in length, containing 504 genes) (Fig. 2). A total of 741 genes that may related to the disease resistance traits of soybean cyst nematodes were identified during four specific genomic region. These four genomic regions had been identified as candidate regions for genes related to soybean cyst nematode resistance.

Based on the comparison results and gene functional annotation data in the NCBI database, 351 genes with functional annotations were predicted among four candidate fragments. Through gene functional clustering analysis, it was found that these 351 genes were annotated into 10 categories (Table 4), including RNA biosynthesis, RNA processing, protein biosynthesis, carbohydrate metabolism, lipid metabolism, enzyme classification, plant hormone action, isoflavone metabolism, chalcone metabolism, and chloroplast redox homeostasis. For instance, Glyma.18G287000, located in Chr18, was found to be associated with disease-resistance RPP13. Glyma.18G252800, also located in Chr18, was found to be associated with pentatricopeptide repeat-containing protein. Glyma.18G267900, also located in Chr18, was found to be associated with isoflavone7-O-methyltransferase-like protein. Glyma.18G287400, also located in Chr18, was found to be associated with 2-oxoglutaratedehydrogenase. Glyma.18G297200 and Glyma.18G298200, both located in Chr18, were found to be associated with Cytokininhydroxylase and 2-alkenalreductase (NADP(+) - dependent), respectively. Glyma.18G281700 was predicted to be associated with probable disease resistance protein. Glyma.18G252300 was predicted to be associated with ethylene-responsive element binding protein. Glyma.18G287100 was predicted to be associated with NBS-LRR disease-resistance protein. Glyma.18G267500 was predicted to be associated with isoflavone7-O-methyltransferase. Glyma.18G269500 was predicted to be associated with disease resistance protein(CC-NBS-LRRclass)family. The last three identified genes in Chr 18 were Glyma.18G267800, Glyma.18G281500 and Glyma.18G285800, which were predicted to be associated with O-methyltransferase family protein, Disease resistance protein and chalconereductase, respectively. The only identified gene located in Chr02 was Glyma.02G211400, which was predicted to be associated with serine/threonine-proteinkinase.

Based on gene functional annotation, the metabolic pathways involved in candidate genes were analyzed, and combined with relevant literature reports, it was found that five genes may be related to soybean disease resistance, namely the protein 2-oxoglutarate dihydrogenase containing pentapeptide repeat sequences, mitochondria, 2-alkaline reductase (NADP (+) dependent), and chalcone reductase serine/threonine protein kinase. The impact of these candidate genes on soybean disease resistance will be verified in subsequent work. Using vector that can overexpress or silence specific gene to infect plant to verify the effect of gene on soybean disease resistance. However, the resistance might be results from interaction between multiple genes; therefore, more extensive studies need to be designed and conducted conducted to verify our results.

RT-PCR

Selected the screened genes, design primers, extract RNA, and perform reverse transcription. Using ACTIN7 as the internal reference gene, perform relative quantitative analysis of the data by 2^−ΔΔCt approach. Relative level of maternal parent is shown in Fig. 3, with the maternal parent being 1. The expression levels of selected genes such as pentapeptide repeat protein, 2-hydroxyglutarate dihydrogenase, 2-olefin reductase (NADP (+) dependent), chalcone reductase, and serine/threonine protein kinase increase with the development of soybean cyst nematodes. Based on gene annotation and analysis results, the above five candidate genes exhibit non synonymous mutations in exons at the DNA level, and are highly expressed at the RNA level in roots. Therefore, there is a strong belief that this genes are closely related to the disease resistance of soybean cyst nematodes, especially with the increase of expression levels during soybean cyst nematode infection and development, which can be used for gene cloning, further confirming their biological functions in the development process of soybean cyst nematodes (Fig. 3, Table S1).

Preliminary gene mapping

The combination of BSA analysis method and high-throughput sequencing technology has shown great advantages in gene localization, and has become the main research strategy for fine mapping of new generation genes. Whole genome re sequencing is carried out on DNA mixed pools constructed by individuals with extreme phenotypes, and fine mapping of target genes is carried out by screening SNP, InDel and other markers closely associated with target traits; the method is extensively used for mining candidate genes for target traits. Zhang et al. (2020) used the BSA-Seq method to locate the yellow 622 multi lobule gene in soybean within the 0–250,000 bp, 1,510,000–3,480,000 bp, and 5,570,000–6,710,000 bp segments on chromosome 11, containing a total of 690 genes and identifying 6 potential genes. Liu et al. (2020) used BSA-Seq method to discover that the candidate regions for controlling sunflower stem are 0.46 Mb on chromosome 3 and 5 kb on chromosome 14, and screened 8 candidate genes. Han et al. (2018) utilized the BSA-Seq method for locating the QTL related to aphid resistance on chromosome 5 of cucumber, containing 43 genes, of which 16 may be related to aphid resistance. Chen et al. (2021) used BSA-Seq method for locating two candidate regions for rice blast disease on chromosome 6 at 10,082–11,397 kb and chromosome 11 at 120–266 kb.

Ochar et al. (2022) conducted bulk segregant analysis (BSA) and whole-genome resequencing, and obtained 12 potential target trait-related genome regions (20.32 Mb long) by using euclidean distance (ED) correlation algorithm. After comparing sequencing data between mutant and wild-type pools, there was just one single nucleotide mutation (C:G >T:A) in exon 1 of Glyma.19G207100 linked to soybean crinkled leaf.

Chen et al. (2022) used the BSA-Seq method to annotate a total of 143 anthracnose resistant genes on sorghum chromosome 5 within the 1.30 Mb region, including 49 non synonymous mutated genes and 16 frameshift mutated genes. These studies indicate that combining whole genome resequencing with BSA methods can quickly and effectively locate and explore the agronomic traits and stress resistance genes of crops. This study used H-10 (resistant) and Chidou4 (susceptible) as parents for hybridization to construct an F₂ segregation population. Using whole genome resequencing combined with BSA method, the resistance genes controlling the third race of soybean cyst nematode disease were preliminarily located. Research has found that the genes controlling disease-related traits are mainly concentrated in the region of 38,524,128–39,849,988 bp on chromosome 2, 27,821,012–29,612,574 bp on chromosome 3, 308–348,214 bp on chromosome 10, and 53,867,581–58,017,852 bp on chromosome 18, totaling 741 genes. Through gene annotation and metabolic pathway analysis of candidate genes, combined with relevant literature reports, 5 genes have been identified, which are related to soybean disease resistance.

Function of candidate genes

There are five candidate genes (including pentapeptide repeat protein, 2-hydroxyglutarate dehydrogenase, 2-olefin reductase(NADP(+)dependent), chalcone reductase serine/threonine protein kinase) related to disease resistance within the localization interval of this study.

Serine/threonine protein kinases are an important class of signaling molecules. Upon stressful conditions like salt, pest feeding, drought, or trauma, as well as other cytokines and hormone stimuli, serine/threonine protein kinases quickly phosphorylate and activate at the serine and threonine residues,and further activate downstream signaling molecules through cascade phosphorylation, activating specific signaling pathways, Ultimately, external signals are transmitted to the nucleus, activating or inhibiting specific gene levels (Zeng et al., 2021a).

Soybean cyst nematode resistance accounts for the complicated quantitative trait adjusted via major and several minor genes. Previous studies on linkage analysis have shown that the main locus of soybean cyst nematode disease is located on chromosomes 08, 11, and 18. For example, the three genes Glyma18g02580 (amino acid transporter), Glyma18g02590 (solid NSF attachment protein, SNAP18), and Glyma18g02610 (wound induced protein WI12) on the highly resistant material Peking are on chromosome 18, whereas Rhg1 is on chromosome 18 (G-linked group),which controls the resistance of multiple SCN physiological races. Rgh1 is involved in regulating host plant specific defense responses, Causing obstruction(necrosis) of syncytial development, thereby affecting the development and reproduction of soybean cyst nematodes (Liu et al., 2012). The disease resistance gene Glyma08g11490 (SHMT) in PI88788 (resistance) germplasm was located at the Rhg4 locus and was involved in encoding amino acid transporters α-SNAP and WI12 proteins are related to the SCN resistance mechanism (Liu et al., 2017). Li et al. (2016) found through linkage and association analysis that a QTL SCN3-11 resistant to cyst nematode (race 3) can explain 5.8% of the resistance variation in RIL populations derived from Zhonghuang 13 (susceptible) and Zhongpin 03-5373 (resistant). Aggregating the rhg1 site can significantly improve resistance. Through gene structure and expression analysis, it was further demonstrated that Glyma.11G234500 at SCN3-11 locus was a new disease resistance gene GmSNAP11 homologous to GmSNAP18 at rhg1-b locus (Lakhssassi et al., 2017).

According to our findings, resistance site of soybean resistant material H-10 is located within the range of 38,524,128–39,849,988 bp on chromosome 02, while the resistance site on chromosome 18 was enriched within the range of 53,867,581–58,017,852 bp. This may be a new resistance gene enrichment region, and the excavation and functional study of the new genes will provide a basis for the soybean cyst nematode resistance mechanism, and the aggregation of new resistance genes. In the future, overexpressing and silencing vectors could be designed and transformed to soybean to confirm its potential role in plant and nematode interaction.

Researchers have hypothesized that mechanisms contributing to soybean cyst nematode (SCN) resistance are not only species-specific but also genotype-specific (Kofsky, Zhang & Song, 2021). Moreover, the culturing environment may also affect the development of SCN resistance. Therefore, under different circumstance such in real fields, different results might be obtained from plants.