Identification of Candidate Genes Controlling Soybean Cyst Nematode Resistance in “Handou 10” Based on Genome and Transcriptome Analyzes

Soybean cyst nematode (SCN; Heterodera glycines Ichinohe) is a highly destructive pathogen for soybean production worldwide. The use of resistant varieties is the most effective way of preventing yield loss. Handou 10 is a commercial soybean variety with desirable agronomic traits and SCN resistance, however genes underlying the SCN resistance in the variety are unknown. An F2:8 recombinant inbred line (RIL) population derived from a cross between Zheng 9525 (susceptible) and Handou 10 was developed and its resistance to SCN HG type 2.5.7 (race 1) and 1.2.5.7 (race 2) was identified. We identified seven quantitative trait loci (QTLs) with additive effects. Among these, three QTLs on Chromosomes 7, 8, and 18 were resistant to both races. These QTLs could explain 1.91–7.73% of the phenotypic variation of SCN’s female index. The QTLs on chromosomes 8 and 18 have already been reported and were most likely overlapped with rhg1 and Rhg4 loci, respectively. However, the QTL on chromosome 7 was novel. Candidate genes for the three QTLs were predicted through genes functional analysis and transcriptome analysis of infected roots of Handou 10 vs. Zheng 9525. Transcriptome analysis performed also indicated that the plant–pathogen interaction played an important role in the SCN resistance for Handou 10. The information will facilitate SCN–resistant gene cloning, and the novel resistant gene will be a source for improving soybeans’ resistance to SCN.


INTRODUCTION
Soybean cyst nematode (SCN, Heterodera glycines Ichinohe) is a devastating pest affecting soybean (Glycine max [L.] Merr.) production worldwide (Riggs, 1977;Koenning and Wrather, 2010;Tylka and Marett, 2014;Mitchum, 2016;Miraeiz et al., 2020;Peng et al., 2021). SCN has caused approximately 36% of yield losses in the total soybean production from 1996-2014 in the United States (Kim et al., 2016). In Chifeng area in China, the highest yield reduction rate reached 44.43% because of SCN during 2014-2019 . Moreover, a highly virulent of SCN has been observed in China (Lian et al., 2019). SCN is a soil-borne pathogen and pest management is difficult. SCN management includes crop rotation, pesticide application, biological control, pest-resistant varieties, etc.; however, breeding for resistant varieties is the most effective method (Usovsky et al., 2021). Riggs and Schmitt (1988) reported 16 physiological races that could be differentiated, and the pathogenicity of each race was different. As in Missouri in the United States, the dominant race in the Huanghuai Valley in China is race 2, which was evolved from race 1 (Lu et al., 2006;Mitchum et al., 2007;Lian et al., 2016;Howland et al., 2018). Accordingly, race 2 was used to screen new varieties for SCN resistance in the Huanghuai Valley, unfortunately, most varieties were susceptible to SCN race 2. 1 Most studies showed that SCN resistance is a quantitative trait and controlled by multiple genes (Concibido et al., 2004;Wang, 2019); more than 200 QTLs have been mapped on 20 chromosomes. 2 Two major QTLs, rhg1 (Peking-type rhg1-a and PI88788-type rhg1-b) and Rhg4 (GmSHMT08), were cloned and functionally analyzed (Cook et al., 2012;Liu et al., 2012;Liu et al., 2017). Three genes (Glyma.18G022500, Glyma.18G022500, and Glyma.18G022700) around the rhg1-b locus and the gene GmSHMT (Glyma.08G108900) around the Rhg4 locus were identified as major-effect gene to SCN resistance (Cook et al., 2012;Liu et al., 2012;Guo et al., 2020). The PI88788-type requires at least 5.6 copies of rhg1-b (Cook et al., 2012;Patil et al., 2019), whereas the Peking-type requires rhg1-a and Rhg4 for SCN resistance (Patil et al., 2019). In the United States, most SCN resistant cultivars are from PI88788 and this has reduced the effectiveness of SCN prevention (Mitchum, 2016). For breeding resistant cultivars, it is vital to identify new quantitative trait loci (QTL) and genes underlying resistance (Liu et al., 2019) and to broaden the genetic basis for improving soybeans' resistance to SCN. In recent years, Numerous researchers have identified genes (Guo et al., 2020) and reveled a complex regulatory network involved in SCN resistance (Guo et al., 2020;Kofsky et al., 2021;Shi et al., 2021).
Handou 10 was first identified as being resistant against SCN race 1 in routine cultivar testing in 2008 by our team and was registered by the Hebei Variety Approval Committee in 2011. The yield of the variety was significantly better than the control, having the desired agronomic traits in the Hebei uniform test and pre-releasing tests in 2007-2010. However, the origin and inheritance of the SCN resistance in Handou 10 was unknown. The objective of this study was to identify the QTL and candidate genes controlling the resistance and provide the basis for molecular marker development and markerassisted breeding.

Plant Materials
A recombinant inbred line (RIL) mapping population of 392 F 2:8 lines was developed by single seed descent (SSD) from the cross between Zheng 9525 and Handou 10. Zheng 9525 was cultivated by Henan Academy of Agricultural Sciences, whereas Handou 10 was cultivated by Henan Jintun Seed Industry Co., Ltd., and Handan Academy of Agricultural Sciences. The seeds of nine differential cultivars, PI88788, Peking, PI437654, PI209332, PI548316, PI89772, PI90763, Pickett, and Lee, were obtained from Henan Academy of Agricultural Sciences.

Soybean Cyst Nematode Resistance Identification
Handou 10, Zheng 9525, and nine differential cultivars were evaluated for the resistance to SCN [HG types 2.5.7 (race 1), 1.2.5.7 (race 2), and 1.2.3.5.6.7 (race 4)] in a climate room in Henan Academy of Agricultural Sciences. Plastic cups (Ø 6 cm × h 12 cm) were filled with soil infected by SCN [HG types 2.5.7 (race 1), 1.2.5.7 (race 2), and 1.2.3.5.6.7 (race 4)], respectively. Lees were planted in SCN-infested soil, and after 30 days, we collected cysts from the roots using a 710-250 µm sieve tower. Cysts were collected from the 250 µm sieve and rinsed. The eggs were collected by breaking open cysts with a rubber stopper and collecting the eggs on a sieve stack consisting of 250 µm -75 µm -25 µm sieves. The mixture from the 25 µm sieve was backwashed into a 50 mL plastic conical tube. A 40% sucrose solution was added to the tubes, stirred, and centrifuged at 2,000 rpm for 5 min. Eggs in the middle layer or supernatant were then collected over a 25 µm sieve. Handou 10, Zheng 9525, and nine differential cultivars were transplanted with five replicates. Each replicate was one plant in a plastic cup (Ø 6 cm × h 12 cm). Five days after transplantation, seedlings were inoculated with about 4,000 eggs per cup. The plants grew at 70 to 80% relative humidity, 28-24 (L/D), and a photoperiod of 16 h: 8 h (L:D) and were watered daily.
The SCN resistance of the 392 RILs and the parents (Zheng 9525 and Handou 10) were evaluated for SCN resistance against HG type 1.2.5.7 (race 2) and HG type 2.5.7 (race 1) in a climate room. Five days after sowing, two plants of each line were transplanted into a plastic cup (Ø 6 cm × h 12 cm) with three replicates, each replicate had two plants. Five days after transplantation, seedlings were inoculated with about 4,000 eggs per cup. The growth conditions were the same as the above.
Thirty days after inoculation, nematode cysts were collected from the roots of each replicate and counted by an image analysis software (Wang et al., 2014). A female index (FI) was calculated as follows: FI (%) = (average number of cysts on each line/average number of cysts on Lee) × 100. FI was used as phenotype data for QTL analysis. The lines were rated as resistant (FI < 10), moderately resistant (10 ≤ FI < 30), moderately susceptible (30 ≤ FI < 60), or susceptible (FI ≥ 60) to classify the response to SCN.

DNA Preparation and Whole-Genome Re-sequencing
Leaf samples of each progeny line and their parents were collected at the seedling stage. DNA was extracted by the plant genomic DNA kit [TIANGEN Biotech (Beijing) Co., Ltd.]. Zheng 9525 and Handou 10 were sequenced by Illumina HiSeq 4000, whereas 392 RIL 2:8 lines from Zheng 9525 × Handou 10 were sequenced by HiSeq X Ten at Huada Gene Technology Co., Ltd. (Shenzhen, China). The sequences of the parents were aligned with the reference genome (Gmax_275_Wm82.a2) using SOAP2 (Li et al., 2009a) and single-nucleotide polymorphisms (SNPs) between the parents were detected by SOAPsnp (Li et al., 2009b). The SNPs between the parents were identified and filtered: (1) Mass value is greater than 20; (2) At least three reads are supported; (3) Heterozygous sites are removed. The pseudomolecules of parental genome sequences were obtained by SNPs between parents with the reference genome. The reads of the progeny population were compared with pseudomolecules of the parental genome sequences using SOAPaligner.

Genetic Map Construction
Instead of using the SNPs as such for linkage mapping, a sliding window-based approach was used to identify bin markers where consecutive SNPs were merged into one bin. We extract the bin area on each chromosome as bin marker, only the bin markers without segregation distortion are selected to construct a map. A genetic map with the bin markers was constructed with JoinMap 4.0 using the maximum likelihood mapping algorithm (Van Ooijen, 2006). Groups were created depending on LOD scores ≥ 3.0 and a maximum distance of 50 cM. The resulting linkage groups were assigned to specific chromosomes according to the reference genome (Gmax_275_Wm82.a2). Regression mapping was used as the mapping algorithm with Kosambi's mapping function to convert recombination frequency into map distance.

Quantitative Trait Loci Mapping
Additive QTLs of SCN resistance with FI to HG types 2.5.7 (race 1) and 1.2.5.7 (race 2), were analyzed using the composite interval mapping (CIM) method in WinQTLCart 2.5 software . The walking speed for CIM was 1 cM and the LOD threshold at the 5% probability level was determined by a 1,000 permutation test.

Transcriptome Sequencing
The seedlings with consistent growth were selected to transplant into a plastic cup with sterile soil after Handou 10 and Zheng 9525 were sown in vermiculite, with one plant per cup. Five days later, well-developed plants of the same size were selected to inoculate 4,000 or 0 SCN HG type (1.2.5.7) (race 2) eggs per cup, with three replicates per genotype of each treatment, and five plants for each replicate. Ten days after inoculation, roots were collected, washed, frozen in liquid nitrogen and stored at -80 • C until use. Total RNA was extracted by TRIzol (Invitrogen). Each biological replicate contained pooled roots from five individual plants. The RNA transcriptome sequencing and preliminary data analyzes were carried out by Shenzhen Huada Gene Technology Co., Ltd. (China). For each replicate, a mRNA library was constructed and sequenced using the DNBSEQ platform. Adaptor reads, reads with an unknown base N greater than 5%, and low-quality reads (bases with a quality value of less than 15 account for more than 20% of the total bases in the reads) were filtered out of the raw data to obtain high-quality (clean) reads using SOAPnuke (v1.4.0) (Chen et al., 2018). Clean reads were mapped to the soybean reference genome (G.max Wm82.a2.v1) by HISAT (Hierarchical Indexing for Spliced Alignment of Transcripts, V2.1.0) (Kim et al., 2015) and aligned using Bowtie 2 (v2.2.5) (Langmead and Salzberg, 2012). The genes and transcripts were calculated using RSEM (v1.2.8) (Li and Dewey, 2011). Significant differentially expressed genes (DEGs) were obtained with false discovery rates (FDRs) of ≤ 0.05.

Functional Annotation and Pathway Enrichment
Significant DEGs were annotated with Kyoto Encyclopedia of Genes and Genomes (KEGG). 3 The statistical enrichment of DEGs in KEGG pathways was accomplished with the R (version 3.1.1) function phyper. 4 Then FDR was performed on p-value, and q-values of ≤ 0.05 was considered to be significantly enriched.

Candidate Genes Analyzes of Major Quantitative Trait Loci Intervals
The QTLs that could be detected for the resistance of Handou 10 to SCN HG type 2.5.7 (race 1) and 1.2.5.7 (race 2) were considered major QTLs. The physical positions of the major QTL intervals were identified according to the above genetic map and the reference genome ("Williams 82.a2.v1") (Jiang et al., 2018). The DEGs between the infected roots of Handou 10 and Zheng 9525 in the major QTL intervals were the candidate genes controlling SCN resistance in Handou 10. and Rhg4 were referred to the reports from Shi et al. (2015) and Kadam et al. (2016). The PCR reaction mixture was prepared according to the instructions of the KBioscience (Herts, United Kingdom). The program was set to hold at 94 • C for 15 min, followed by 10 touch-down cycles of 20 s at 94 • C and 1 min at 65-59 • C (dropping 0.6 • C per cycle), and then 23 cycles of 20 s at 94 • C, 1 min at 57 • C. The PCR amplification product was read with a PHERAstar SNP typing detector, and SNP alleles were determined based on the ratio of fluorescence signals.

Analysis of Whole-Genome Re-sequencing
Zheng 9525 and Handou 10 were re-sequenced using an Illumina HiSeq4000 platform with sequencing average depths of 11.94× and 16.94×, respectively (Supplementary Table 1 Table 2). The sequencing data showed that the average mapped reads was 12.74 M and the average bases were 1.91 G. A total of 607,635 SNPs were identified in the RIL population and 8,593 breaking points were detected using the sliding window approach (Supplementary Table 3).

Construction of High-Density Genetic Map
A high-density genetic map was constructed by joining 5,233 bin markers (missing < 20%, Supplementary Table 4). The map comprised 4763.23 cM with an average distance of 0.91 cm between adjacent markers (Table 3). Chromosome 18 had the highest number of bin markers (405 markers). The linkage map length was the largest for chromosome 5 (325.35 cM), and the smallest was for chromosome 4 (89.47 cM) ( Table 3).

Quantitative Trait Loci Mapping
Seven QTLs for SCN resistance with FI were detected on six chromosomes by CIM, which were designated as SCN_7_1, SCN_8_2, SCN_12_3, SCN_15_4, SCN_18_5, SCN_18_6, and SCN_20_7, respectively ( Table 4). The resistance alleles of five QTLs (SCN_7_1, SCN_8_2, SCN_18_5, SCN_18_6, and SCN_20_7) were derived from the resistant parent Handou 10, and the resistance alleles of the remaining two QTLs (SCN_12_3 and SCN_15_4) were derived from the susceptible parent. Three QTLs (SCN_7_1, SCN_8_2, and SCN_18_6) were detected in race 1 and race 2. The percentage of the explained variance of the identified QTLs varied from 1.91 to 7.73%. According to their physical position in the genome, the QTL intervals of SCN_8_2 (6.75-10.35 Mb) and SCN_18_6 (0.05-3.16 Mb) overlapped with the rhg1 and Rhg4 loci (see Footnote 2), respectively.

Transcriptomic Analysis
We analyzed the transcriptomes of the resistant and susceptible parents infected with race 2. Each sample produced 6.45 G data on average (Supplementary Table 5). The average comparison rate between each sample and reference genome was 80.70%. A total of 912 and 981 significant DEGs were detected between the infected and uninfected roots in Handou 10 and Zheng 9525, respectively (Figure 2). In Handou 10, 312 and 600 upand down-regulated genes were detected, respectively. In Zheng 9525, 395 and 586 up-and down-regulated genes were detected, respectively. A total of 4424 DEGs were detected between the infected roots of Handou 10 and Zheng 9525.
The KEGG analysis showed that 2100 DEGs were assigned to 132 pathways, and the plant-pathogen interaction pathway (ko04626) was the top enriched pathway (Figure 3) with 213 DEGs between the SCN infected roots of Handou 10 and Zheng 9525. The results suggested that the plant-pathogen interaction played the most important role in the SCN resistance for Handou 10. MPK signaling-plant (ko04016), phenylpropanoid biosynthesis (ko00940), plant hormone signal transduction (ko04075) and other pathways also played important roles in the resistance to SCN (Figure 3).

Candidate Gene Analyzes of SCN_7_1, SCN_8_2, and SCN_18_6
According to the physical positions of the major QTL intervals, there were 6, 23, and 12 candidate genes from the DEGs between the infected roots of Handou 10 and Zheng 9525 for SCN_7_1, SCN_8_2, and SCN_18_6, respectively ( Table 5). The candidate gene Glyma.08G108900 at SCN_8_2 and candidate genes Glyma.18G022400, Glyma.18G022500, and Glyma.18G022700 at SCN_18_6 have been reported for SCN resistance (Guo et al., 2019;Cook et al., 2012;Liu et al., 2012). Glyma.08G108900 was Frontiers in Plant Science | www.frontiersin.org FIGURE 1 | Female index (FI) and rating distributions of SCN resistance to HG type 2.5.7 and HG type 1.2.5.7 for the RILs derived from Zheng 9525×Handou 10, respectively. (A) The FI distribution of SCN resistance to HG types 2.5.7; (B) the rating of SCN resistance to HG types 2.5.7; (C) the FI distribution of SCN resistance to HG types 1.2.5.7; (D) the rating of SCN resistance to HG types 1.2.5.7.   FI: female index. Physical position was based on the Glycine max genome assembly version Wm82.a2v1. The significant LOD thresholds were estimated by 1,000 permutations at the 5% significance level as follows: 1.79 for HG type 1.2.5.7 and 2.49 for HG type 2.5.7. R 2 : percentage of phenotypic variance explained by a QTL. Additive effect: The positive value implies that Handou 10 decreases the phenotypic value (FI); the negative value implies that Zheng 9525 decreases the phenotypic value. located in the Rhg4 locus. Glyma.18G022400, Glyma.18G022500, and Glyma.18G022700 were located in the rhg1 locus. These three genes (Glyma18g02580, Glyma18g02590, and Glyma18g02610 in Wm82.a1) in rhg1 were all expressed in the resistance of Handou 10 according to the transcriptome data ( Table 5).

DISCUSSION
Handou 10 is a variety with a broad-spectrum resistant to SCN. Our recent experiment result showed that Handou 10 was resistant to SCN HG type 7 (race 3). In this study, Handou 10 was identified as being resistant against HG type   DNA \Assay Rhg1-2 a Rhg1-5 a GSM381 b GSM383 b Rhg4-3 a Rhg4-5 a GSM191 b DD7_1 C DD7_5 C DD7_8 C DD7_9 C DD7_11 C DD7_16 C   Handou 10  GG  CC  GG  GG  TT  GG  GG  AA  TT  TT  AA  TT  AA   Zheng 9525  CC  GG  TT  CC  AA  CC  CC  CC  AA  AA  GG  CC  CC   PI437654  GG  CC  GG  GG  AA  GG  GG  CC  TT  TT  AA  TT  AA   Peking  GG  CC  GG  GG  TT  GG  GG  AA  AA  AA  GG  CC  AA 2.5.7 (race 1), and moderately resistant against HG type 1.2.5.7 (race 2); additionally, its resistance was identified as being Peking-type based on the analysis of the QTL (Table 4), KASP markers (Table 6) and the CNV at the rhg1/Rhg4 (Supplementary Figure 1). SCN_8_2 and SCN_18_6 overlapped with the resistance loci Rhg4 and rhg1-a, respectively. This two additive QTLs to SCN resistance with FI were also detected by IciMapping 4.2 software (Supplementary Table 7). SCN_8_2 and SCN_18_6 were interactive for SCN resistance to HG types 2.5.7 (race 2) (Supplementary Table 8). According to the CNV analyzes (Suvakov et al., 2021) for the whole-genome re-sequencing, the CNV at the rhg1 and Rhg4 for Handou 10 is 2.95 and 2.05, respectively (Supplementary Figure 1). The resistance to SCN and CNV at the rhg1/Rhg4 in Handou 10 is similar to PI404166, PI 437679, and PI089772 (Patil et al., 2019). The rhg1-a, sometimes in combination with the Rhg4, provides strong resistance to SCN (Guo et al., 2019). Three genes (Glyma.18G022400, Glyma.18G022500, and Glyma.18G022700) are responsible for the resistance provided by rhg1-b (Cook et al., 2012). In our study, these three candidate genes in rhg1-a might be regulate the SCN resistance of Handou 10. However, the mechanisms of rhg1-a and rhg1b to SCN resistance are not same (Cook et al., 2012;Patil et al., 2019). For rhg1-b, overexpression of the individual genes in roots was ineffective and SCN resistance is conferred by copy number variation. The Peking-type requires rhg1a and Rhg4 for SCN resistance (Patil et al., 2019). In this study, the expression of Glyma.08G108900 located in the Rhg4 locus was verified by qPCR (Supplementary Figure 2) and this was basically consistent with transcriptome analysis. Glyma.08G108900 was responsible for the resistance to SCN and the function has been verified. The resistant mechanism of Handou 10 is complicate and we are breeding new superior lines with it now. SCN_7_1 is an important QTL to SCN resistance for Handou 10. Some QTLs (Webb et al., 1995;Ferreira et al., 2011;Abdelmajid et al., 2014;Chang et al., 2016;Vuong et al., 2015; related to SCN resistance on chromosome 7 were reported at www.soybase.org. The physical location of SCN_7_1 was basically the same as the marker php02301a mapped in PI437654 (Webb et al., 1995), but far from the QTLs detected by other researchers (see Footnote 2). There are six DEGs in the SCN_7_1 region and these candidate genes are not previously described as SCN resistance genes. Among the DEGs at SCN_7_1, Glyma.07G139700 and Glyma.07G139800 both code Glutathione S-transferases (GSTs), which were up-regulated between infected and uninfected SCN in Handou 10 and down-regulated in Zheng 9525 according to transcriptome data. Recently, the expression of Glyma.07G139800 was verified by qPCR (Supplementary Figure 2) and this was basically consistent with transcriptome analysis in our study. GSTs are multifunctional enzymes which play a crucial role in cellular detoxification and oxidative stress tolerance (Rezaei et al., 2013). GST was elevated in the SCN-infected roots relative to uninoculated roots . The over-expression of a GST gene from wild soybean (Glycine soja) enhances drought and salt tolerance in transgenic tobacco (Ji et al., 2010). However, the relationship between glutathione metabolism and the disease resistance of Handou 10 still needs to be further studied.
Soybean resistance to SCN is regulated by multiple genes. In our study the plant-pathogen interaction pathway was the most enriched KEGG pathway between infected Handou 10 and Zheng 9525, and between infected Zheng 9525 and uninfected Zheng 9525, whereas the MAPK signaling pathway was the most enriched KEGG pathway between infected and uninfected Handou 10. Plants and pathogens should be studied together as an interacting system (Peyraud et al., 2017). Soybean root cells undergo dramatic morphological and biochemical changes after being infected with SCN . The development of cyst nematodes in the infected SCN roots in Handou 10 was slower than in Zheng 9525 after 10 days of inoculation (Supplementary Figure 3). Many DEGs may be related to the development of soybean cyst nematode. These DEGs genes (Glyma.07G139700, Glyma.07G139800, Glyma.08G108900, Glyma.08G118900, Glyma.08G119200, Glyma.18G022400, Glyma.18G022500, and Glyma.18G022700) might be the important candidate gene to SCN resistance according to KEGG analysis and genes function.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NGDC CAR005982.

AUTHOR CONTRIBUTIONS
WL, JW, and HW conceived the project and designed the experiments. HW, YL, JL, HL, YW, CL, SW, and HZ performed the experiments. HW, WL, JW, and QS participated in the data analysis and manuscript revised. All authors contributed to the article and approved the submitted version.