Characterization of Insect Resistance Loci in the USDA Soybean Germplasm Collection Using Genome-Wide Association Studies

Management of insects that cause economic damage to yields of soybean mainly rely on insecticide applications. Sources of resistance in soybean plant introductions (PIs) to different insect pests have been reported, and some of these sources, like for the soybean aphid (SBA), have been used to develop resistant soybean cultivars. With the availability of SoySNP50K and the statistical power of genome-wide association studies, we integrated phenotypic data for beet armyworm, Mexican bean beetle (MBB), potato leafhopper (PLH), SBA, soybean looper (SBL), velvetbean caterpillar (VBC), and chewing damage caused by unspecified insects for a comprehensive understanding of insect resistance in the United States Department of Agriculture Soybean Germplasm Collection. We identified significant single nucleotide (SNP) polymorphic markers for MBB, PLH, SBL, and VBC, and we highlighted several leucine-rich repeat-containing genes and myeloblastosis transcription factors within the high linkage disequilibrium region surrounding significant SNP markers. Specifically for soybean resistance to PLH, we found the PLH locus is close but distinct to a locus for soybean pubescence density on chromosome 12. The results provide genetic support that pubescence density may not directly link to PLH resistance. This study offers a novel insight of soybean resistance to four insect pests and reviews resistance mapping studies for major soybean insects.


INTRODUCTION
Insect damage is one of the major limiting factors for soybean (Glycine max (L.) Merr.) production as insects vector viruses and cause damage by feeding on foliage, vascular sap, stems, roots, pods, and seeds (Steffey, 2015). In order to manage yield losses, application of insecticides has often been the first tool in management. For example, insecticide usage increased in response to soybean aphid (SBA, Aphis glycines Matsumura) dissemination in the north central region of the USA (Coupe and Capel, 2015), where 80% of soybean production occurs. Insect damage also devastated production in the southern parts of the USA where 234 million dollars of soybean losses were reported despite investing 279 million dollars on insect management (Musser et al., 2016;Ortega et al., 2016). Soybean is one of the most important field crops for providing dietary protein and oil worldwide, and it is estimated that a 50% increase in soybean production is needed to meet the population expansion by 2030 (Hartman et al., 2011;Ainsworth et al., 2012). In addition to insecticidal use, an alternative way to control insect damage is through genetic resistance. Insect resistance has been widely studied and applied in breeding program for several crops (Edwards and Singh, 2006;de Morals and Pinheiro, 2012). Characterization of novel resistance to different insects may be essential to sustain soybean productivity.
Plants defend themselves against insects in a number of ways and mechanisms of resistance have been described as antibiosis, non-preference, and tolerance or the interaction of these three factors (Painter, 1941). While antibiosis describes the ability of resistant plants to restrict insect growth or propagation, nonpreference or later referred to as antixenosis, deters insects from feeding and propagating on resistant plants . Both antibiosis and antixenosis resistance could be genetically governed by the same locus in a resistant genotype. For example, one of the major insect resistance quantitative trait loci (QTL) in soybean is the QTL-E (SoyBase QTL name: corn earworm 8-1) on chromosome (chr) 15, on linkage group E (LG-E) (Terry et al., 2000;Hulburt et al., 2004). QTL-E contributes 26% of the antibiotic effect and 20% of the antixenotic effect of resistance to corn earworm (CEW, Helicoverpa zea Boddie) (Boerma and Walker, 2005). QTL-E has been mapped to the same location to the Pb locus that determines the tip phenotype of pubescence, sharp (Pb) or blunt (pb) (Palmer and Xu, 2008;Parrott et al., 2008;Ortega et al., 2016). It has been shown that the Pb locus provides antixenotic resistance by discouraging insect feeding on soybeans with sharp-tipped pubescence (Hulburt et al., 2004), a trait which is rare in domesticated soybean but common in wild soybean (G. soja) (Broich and Palmer, 1981). Another insect resistance QTL is QTL-M (SoyBase QTL name: corn earworm 1-1) on chr 7 (LG-M) which accounts for 22% of the antibiotic effect and 37% of the antixenotic effect of resistance to CEW (Rector et al., 1998(Rector et al., , 1999(Rector et al., , 2000, as well as resistance to other insects including several lepidopteran insects such as soybean looper (SBL, Pseudoplusia includens Walker), velvetbean caterpillar (VBC, Anticarsia gemmatalis Hübner) and a coleopteran insect, Mexican bean beetle (MBB, Epilachna varivestis Mulsant). QTL-M also exhibits synergistic epistasis to QTL-G on chr 18 (LG-G, SoyBase QTL name: corn earworm 6-1) that conditions antibiosis and to QTL-H on chr 12 (LG-H, SoyBase QTL name: corn earworm 1-2) that conditions antixenosis (Rector et al., 1998(Rector et al., , 2000Zhu et al., 2008). Accordingly, QTL-M has become a major breeding target for soybean, and a follow up study that fine-mapped QTL-M between Sat_258 and Satt702 in a 0.25 cM interval .
Along with resistance to defoliator insects, resistance to piercing-sucking insects in soybean has also been well documented. Aphid resistance in soybean was found to be antibiotic and antixenotic in varieties "Dowling" and "Jackson, " but only antixenotic in the variety "PI 71506" (Hill et al., 2004).
The antibiotic QTL in "Dowling" and "Jackson" were mapped approximately to QTL-M, and were assigned as Rag1 (Resistance to Aphis glycines 1) and Rag, respectively (Li et al., 2007). Other SBA resistance QTL include antibiotic Rag2 on chr 13 (Mian et al., 2008;Hill et al., 2009), antixenotic Rag3 and Rag3b on chr 16 (Zhang et al., 2010, and two antibiotic recessive QTL (rag1c and rag4) on chr 7 and 13, respectively . Fine mapping efforts on Rag1 discovered two nucleotide binding leucine-rich repeat (NBS-LRR) genes among 13 candidate genes in a 115 Kb interval on chr 7, and fine mapping on Rag2 identified one NBS-LRR gene along with seven candidate genes in a 54 Kb region on chr 13 (Kim et al., 2010a,b). However, instead of the NBS-LRR gene (Glyma13g26000) within the Rag2 region, another NBS-LRR (Glyma13g25970) that locates on the border of Rag2 was proposed to be the candidate resistance gene based on differential expression analyses (Brechenmacher et al., 2015).
The resolution of linkage mapping has been limited by the density of traditional DNA markers, such as simple sequence repeat (SSR), and by the limited recombination that occurs during bi-parental crossing. In addition, another disadvantage of linkage mapping is the monotonous resistance source from one of the parents. The majority of resistance to defoliator insects (CEW, MBB, SBL, and VBC) was derived from three Japanese accessions (PI 171451, PI 227687, and PI 229358)  and resistances to SBA were reported in a few accessions (Dowling, PI 200538, PI 567324, PI 567537, PI 567541B, PI 567543C, PI 567597C, and PI 587732) (Hill et al., 2004;Bales et al., 2013;Jun et al., 2013;Kim et al., 2014). However, SBA biotypes that overcome different Rag resistance were reported in the field (Hill et al., 2010), and thus identification of additional resistance sources is important. With the advance in biotechnology, methods such as Affymetrix GeneChip and genotyping-bysequencing enabled single nucleotide polymorphisms (SNPs) have become useful as novel genetic markers (Barabaschia et al., 2016). High-density and high-quality SNPs across the entire genome empowers the genome-wide association study (GWAS), which relies on linkage disequilibrium (LD) that has remained through historical recombination in a diverse population. Accordingly, GWAS provides better mapping resolution and also detects multiple genetic sources in a germplasm collection. The power of GWAS has been well recognized in numerous studies, including those on soybean agronomic traits (Zhou et al., 2015) and on soybean disease resistance (Chang et al., 2016a,b,c). There are fewer GWAS focusing on insect resistance in soybean Liu et al., 2016).
The goal of this study was to integrate the insect resistance phenotypes in the Germplasm Resources Information Network (GRIN) (www.ars-grin.gov) and the SNPs in the Soybean Germplasm Collection maintained by the United States Department of Agriculture Agricultural Research Service (USDA-ARS). A comprehensive GWAS for six soybean insects, including beet armyworm (BAW, Spodoptera exigua Hübner), MBB, potato leafhopper (PLH, Empoasca fabae Harris), SBA, SBL, and VBC, was performed to discover novel insect resistance sources. Candidate resistance genes were highlighted for each significant result.

Phenotypic and Genotypic Data Preparation
There were eight cases of insect-related phenotypic data in the GRIN database. Data were categorical records and were transferred to ordinal scales. Phenotypic data of BAW (data contributed by R. L. Nelson), SBL (data contributed by D. Gary, L. Lambert, and R. L. Nelson), and VBC (data contributed by R. L. Nelson) were recorded in five percentage ranges: (i) 0-20% defoliation, (ii) 21-40% defoliation, (iii) 41-60% defoliation, (iv) 61-80% defoliation, and (v) 81-100% defoliation. These five ranges were converted to numbers using 1-5, respectively. Phenotypic data of SBA (Hill et al., 2004;Mensah et al., 2005;Mian et al., 2008;Bhusal et al., 2013; and data contributed by K. Dashiell and L. Hesler) and the defoliation damage caused by unspecified chewing insects (data contributed by L. Hesler) were: (i) resistant, (ii) mostly resistant, and (iii) susceptible. These three categories were converted to 1 to 3, respectively. An original code ranging from 1 to 5, representing little feeding to severe feeding for PLH was used (data contributed by R. L. Nelson). Original values ranging from 1 to 5 (representing damage scales from 1.0 to 1.4, 1.8 to 2.2, 2.6 to 3.4, 3.8 to 4.2, and 4.6 to 5.0, respectively) for MBB were used (Nelson et al., 1987(Nelson et al., , 1988Juvik et al., 1989b;Bernard et al., 1998; and data contributed by T. Elden). The record for CEW contained only 27 resistant soybean accessions (Joshi, 1977), and thus the data of CEW was excluded from our study. One agronomic trait, pubescence density, in the GRIN database was also included in our study. The record for pubescence density has six categories from no pubescence to dense pubescence (glabrous, sparse, semisparse, normal, semi-dense, and dense), and these categories were converted to 1-6, respectively (Nelson et al., 1987(Nelson et al., , 1988Juvik et al., 1989a,b;Coble et al., 1991;Bernard et al., 1998;Hill et al., 2005Hill et al., , 2008Peregrine et al., 2008; and data contributed by R. L. Nelson). Because the categorical entity challenged the ordinal scale for normality transformation, the raw ordinal phenotypic data were used in this study. Genotypic data were the SNPs derived from SoySNP50K project . Preprocessing of the raw SoySNP50K data was identical to our previous pipeline (Chang et al., 2016a,b,c). Soybean accessions with both phenotypic and genotypic data were included in the association analysis. The number of SNPs with minor allele frequency (MAF) above 0.01 and the sample size (the total number of soybean entries) for each of the seven insects are listed in Table 1.

Genome-Wide Association Study (GWAS) and Linkage Disequilibrium (LD) Analyses
The R package "GAPIT2" version 2016.03.01 was applied for the association analyses (Lipka et al., 2012;Tang et al., 2016). Kinship was estimated using the VanRaden method (VanRaden, 2008) in a regular mixed linear model (MLM) (Yu et al., 2006). Principal component analysis and Bayesian information criterion (BIC)-based model selection were applied to determine how many principal components (PCs) should be included for additional population structure correction. A false-discovery rate (FDR) using the Benjamini-Hochberg procedure at 0.05 and a Bonferroni-corrected p-value at 0.05 were used to determine significance in the multiple association tests. If the first step GWAS identified a significant SNP, the SNP was assigned as a fixed covariate in the MLM for the second step of GWAS. The stepping procedure was stopped when no more significant SNPs were indentified, and SNPs used as fixed covariates were considered as true associated signals. Pairwise LD in the flanking region of the true associated signals were calculated using TASSEL5 with a slide window of 500 bp (Bradbury et al., 2007). Candidate resistance genes were examined within the high LD region centering the significant SNP based on the soybean genome assembly version Glyma.Wm82.a1.v1.1 in SoyBase (http://soybase.org).

RESULTS
After extracting consensus phenotypic and genotypic data from the GRIN database and the SoySNP50K, there were seven cases for association analyses (BAW, MBB, SBA, SBL, PLH, VBC, and the defoliation damage caused by unspecified chewing insects). For all the cases, a kinship matrix was included in the MLM. Based on the BIC-based model selection, no PC was included in the MLM for any case except for BAW that used three PCs. Among these seven cases, we identified significant signals for five cases, but not for BAW or the defoliation damage caused by unspecified chewing insects ( Table 1).

Mexican Bean Beetle (MBB)
There were four SNPs on chr 11 that had a FDR below 0.05 and passed the Bonferonni threshold in the first step GWAS (Figures 1A,B). After including the most significant SNP (ss715609849) as a fixed covariate, the second step GWAS identified no more significant SNPs for MBB (Table 1). LD analysis revealed a relatively high LG condition throughout a 14 MB region centering around ss715609849. Pairwise LD values for most SNPs in this region were between R 2 values of 0.2 to 0.4, although a region from 27,000,000 to 30,000,000 bp on chr 11 had higher R 2 values ( Figure 1C). There were two genes (Glyma11g29010 and Glyma11g29287) that contained a leucinerich repeat (LRR) domain within this 3 Mb genomic region.

Potato Leafhopper (PLH)
GWAS for PLH identified a locus on chr 12, where nine significant SNPs passed FDR and the Bonferroni threshold (Figures 2A,B). After fixing the most significant SNP (ss715612746) as a covariate, there was no more significant SNPs. This result indicated that the other eight SNPs might be in LD to ss715612746, and indeed the LD analysis found a very narrow interval (around 37,036,017-37,356,120 bp) flanked by these SNPs on chr 12 (Figure 2C), and within this region, there were two LRR domian-containing genes (Glyma12g33930 and Glyma12g34020) and one myeloblastosis (MYB) transcription factor gene (Glyma12g33911). In addition, the location of ss715612746 was within a pubescence density-related QTL 2-7 (qtuH-1) based on the SNP location and QTL 2-7 location from the genome assembly version Glyma.Wm82.a2 (Du and Fu, 2009). In order to understand the genetic architecture of pubescence density in soybean, the phenotype record in the GRIN database was used to run GWAS. We discovered a significant locus on chr 12 (Figures S1A,B); however, the locus for pubescence density ranged from 33 to 36 Mb whereas the locus for PLH was from 37.04 to 37.36 Mb (Figures 2C,D). This indicates that the locus for PLH resistance and the locus for soybean pubescence density do not overlap ( Table 1).

Soybean Looper (SBL)
Only one SNP (ss715592245) was significant on chr 5 for SBL, and there were no more signals after including ss715592245 as a fixed covariate in the second step of GWAS (Figures 3A,B). LD analysis for ss715592245 identified a high LD region starting around the SNP marker toward about a 3 Mb downstream genomic region ( Figure 3C). There were four LRR domain-containing genes (Glyma05g05730, Glyma05g06140, Glyma05g06231, and Glyma05g07050) and two MYB transcription factors (Glyma05g04900, Glyma05g06410) located in this high LD region from 4,180,000 to 7,180,000 bp.

Velvet Bean Caterpillar (VBC)
Two SNPs on chr 18 with a FDR below 0.05 and the significance was close but below the Bonferroni threshold (Figures 4A,B). After fixing the most significant SNP (ss715629577) as a covariate, there were no more significant signals in the second step of GWAS. LD analysis identified a region around 1,740,000-2,126,000 bp that may harbor resistance genes ( Figure 4C). There were four LRR domain-containing genes (Glyma18g02850, Glyma18g03040, Glyma18g03053, and Glyma18g03066) located within this 4 Mb LD region.

DISCUSSION
Soybean resistance mapping using GWAS has been demonstrated on a variety of diseases (Wen et al., 2014;Vuong et al., 2015;Chang et al., 2016c) and the resistance mechanisms have been detailed for many important soybean diseases (Whitham et al., 2016). However, there are just a few GWAS focusing on soybean insect resistance like that of Frontiers in Plant Science | www.frontiersin.org (C) Regional LD analysis that surrounds the most significant SNP (ss715609849), which is below the red triangle. The pairwise LD between SNPs (gray dots) to ss715609849 follows the right y axis. The orange and gray color for lines indicates SNPs with FDR below or above 0.05, respectively, following the left y axis. The red horizontal dash line represents the minimal significant level at the cutoff of FDR 0.05. The pink background highlights the high LD region where candidate genes were examined.
the common cutworm (CCW, Spodoptera litura Fabricius) Liu et al., 2016). With the availability of GRIN phenotypic data and SoySNP50K genotypic data, the aim of our study was to integrate the information for discovering novel insect resistance sources in the USDA Soybean Germplasm Collection and highlight candidate resistance genes for each insect. In the cases of MBB, PLH, SBL, and VBC, we successfully identified significant SNPs and high LD regions. Because mechanical studies for insect resistance discovered the importance of LRR domain-containing genes (Rossi et al., 1998; and MYB transcription factors (Misra et al., 2010), we reported these two groups of candidate genes within the high LD region. Nonetheless, the possibility that other candidate genes (such as metabolic enzymes or kinases) are involved in insect resistance should not be excluded because insect resistance is complicated and often controlled by multiple mechanisms (Mitchell et al., 2016;Schuman and Baldwin, 2016). Advanced analyses in terms of allelic variation, haplotype diversity, and genomic selection using current SNPs, re-genotyped SNPs with higher marker density, or whole genome re-sequencing data may also provide better understanding of insect resistance in soybean populations (Zhou et al., 2015;Patil et al., 2016). Essentially, molecular cloning and functional analysis for candidate genes are required to confirm their roles in insect resistance.
MBB is a chewing and defoliating insect that damages many legume crops including soybean. MBB is widely distributed in the USA and the southern part of Canada with damage estimates of up to 80% defoliation during the vegetative growth stages of soybean (Nottingham et al., 2016). The soybean cultivar "Davis" was reported to show antixenotic resistance to MBB through coumestrol, an isoflavonoid compound (Burden and Norris, 1992). In addition, two soybean accessions (PI 171451 and PI 229358) were reported to display antibiotic resistance dependent FIGURE 2 | GWAS for potato leafhopper (PLH) and pubescence density. (A) Quantile-quantile plot indicates the fitness of the regular MLM for PLH association analysis. (B) Manhattan plot identifies a significant locus for PLH resistance on chr 12. Green dots represent the nine SNPs with FDR below 0.05. Red horizontal line indicates the Bonferroni threshold. (C) Regional LD analysis for PLH locus that surrounds the most significant SNP (ss715612746), which is below the red triangle. The pairwise LD between SNPs (gray dots) to ss715612746 follows the right y axis. The orange and gray color for lines indicate SNPs with FDR below or above 0.05, respectively, following the left y axis. The red horizontal dash line represents the minimal significant level at the cutoff of FDR 0.05. The pink background highlights the high LD region where candidate genes were examined. (D) Regional LD analysis for pubescence density that surrounds the most significant SNP (ss715612489), shown in red triangle. The pairwise LD between SNPs (gray dots) to ss715612489 follows the right y axis. The orange and gray color for lines indicate SNPs with FDR below or above 0.05, respectively, following the left y axis. The red horizontal dashed line represents the minimal significant level at the cutoff of FDR 0.05. The pink background highlights the high LD region where candidate genes were examined. The genomic region for PLH resistance (C) and pubescence density (D) is exactly the same, but these two loci are not overlapping.
Frontiers in Plant Science | www.frontiersin.org (C) Regional LD analysis that surrounds the most significant SNP (ss715592245), which is below the red triangle. The pairwise LD between SNPs (gray dots) to ss715592245 follows the right y axis. The orange and gray color for lines indicate SNPs with FDR below or above 0.05, respectively, following the left y axis. The red horizontal dashed line represents the minimal significant level at the cutoff of FDR 0.05. The pink background highlights the high LD region where candidate genes were examined. on jasmonic acid regulation (Rufener et al., 1989;Iverson et al., 2001). To the best of our knowledge, there is no genetic mapping study for MBB resistance in soybean or in other host plants.
Our study discovered a significant SNP within a 3 Mb high LD genomic region on the chr 11.
The PLH is a major pest in the USA due to its broad host range including alfalfa, common bean, and soybean. Yield losses caused by the PLH were estimated around $66 per hectare for alfalfa (Lamp et al., 1991;Baker et al., 2015) and $2 million for common bean (Gonzales et al., 2002;Brisco et al., 2015). PLH prefers warmer and drier conditions, and a recent study suggested that PLH infestation would be enhanced with increased temperatures due to climate change (Baker et al., 2015). While insecticide application is one way to control PLH, both antibiotic and antixenotic resistance to PLH were found in common bean (Brisco et al., 2015). Some studies on soybean resistance to PLH proposed that antixenotic resistance is governed by pubescence density (Johnson and Hollowell, 1935;Elden and Lambert, 1992), but others suggested additional pubescence characteristics such as orientation (Boerma et al., 1972) or chemicals in the glandular trichomes might be more important than density (Elden and Lambert, 1992;Ranger and Hower, 2001;Peiffer et al., 2009). In the GWAS for PLH, we identified a SNP within the pubescence density QTL 2-7 (qtuH-1), which is close to pubescence density QTL 2-8 (qtuH-2) (Du and Fu, 2009). In addition to this SNP, two genetic markers for soybean CCW resistance were also located within qtuH-1 included Sat_218 found by linkage mapping (Oki et al., 2012) and a SNP (BARC-043061-08513) found by GWAS . These results indicated qtuH-1 on chr 12 might contribute to CCW and PLH resistance through an antixenotic mechanism that depends on pubescence density. On the other hand, genetic mapping for soybean pubescence density using GWAS in our study provided a better resolution than previous studies using linkage mapping. We showed that Regional LD analysis that surrounds the most significant SNP (ss715629577), which is below the red triangle. The pairwise LD between SNPs (gray dots) to ss715629577 follows the right y axis. The orange and gray color for lines indicate SNPs with FDR below or above 0.05, respectively, following the left y axis. The red horizontal dashed line represents the minimal significant level at the cutoff of FDR 0.05. The pink background highlights the high LD region where candidate genes were examined. the pubescence density locus is distant, about 3 Mb from the PLH locus, and the distance is larger than the average LD size of 0.25 Mb on chr 12 (Vuong et al., 2015). Our results supported the idea that pubescence density may not be the major effect for PLH and another leafhopper (Empoasca terminalis Distant) resistance in soybean (Boerma et al., 1972;Nasruddin and Melina, 2014).
SBL has become one of the most problematic pests for soybean production as the insect developed resistance to carbamates, cyclodienes, DDT, organophosphates, permethrin, and pyrethroids (Boethel et al., 1992). While soybean resistance to SBL has been shown to be antixenotic, the resistance depends on the sharp pubescent tip and the Pb locus on chr 15 (Hulburt et al., 2004). Another study mapped antixenotic resistance to SBL approximately to QTL-M on chr 7 . In addition, antibiotic resistance that relies on chemicals such as afrormosin and phaseol (Caballero and Smith, 1986), growth inhibitors (Smith, 1985;Beach and Todd, 1988), and phytoalexins such as coumestrol and glyceollin were also reported (Liu et al., 1992). Our results discovered a new SNP on chr 5 that may harbor resistance to SBL, and pointed out candidate LRR domain-containing genes.
VBC is another pest that acquired insecticide-resistance (Boethel et al., 1992) and it has been regarded as the most damaging insect for soybean in the southeastern USA (Barbara, 2014). To the best of our knowledge, there is no resistance mapping study for VBC in any crop. Our study identified a 4 Mb LD region on chr 18 with four LRR domain-containing genes. It has been shown that resistance of some soybean cultivars is likely antibiotic and two flavonoids (genistin and rutin) were reported to reduce VBC weights (Piubelli et al., 2005), but other types of candidate genes in the LD region should be considered.
Studies on soybean resistance to the SBA are more abundant than those on all the other insects, and the phenotypic data in the GRIN database contains thousands of accession responses to SBA. Unfortunately, the phenotypic data is unbalanced and only 11 soybean accessions in the GRIN database were classified resistant or mostly resistant. Under the circumstances, GWAS for SBA is problematic because any SNP that is completely shared by these 11 soybean accessions would be considered as a significant signal, and each of these SNPs will co-occur with another SNP. The GWAS result supported this and hundreds of SNPs (Figures S2A,B). When the most significant SNP (ss715596142) was fixed in the second step of GWAS, there were no more significant signals. This observation indicated these SNPs might be confounded with ss715596142 since they were highly correlated to ss715596142 or in LD with ss715596142. We found three LRR domain-containing genes (Glyma07g13440, Glyma07g14810, and Glyma07g14791) and one MYB transcription factor (Glyma07g14480) under a LD region of ss715596142 (from 12,000,000 to 14,800,000 bp) ( Figure S2C). However, the location of ss715596142 on chr 7 differs from Rag1 and the two proposed candidate LRRcontaining genes (Glyma07g06890 and Glyma07 g06920) in the Rag1 locus (Kim et al., 2010a). Instead, ss715596142 is close to Rag3-1, which is one of the resistance QTL reported in soybean variety PI 567541B . Additional evidence will be needed to support our preliminary analysis for SBA.

CONCLUSION
In this study, we report mapping results for the first time for soybean resistance to MBB and VBC located on chr 11 and 18, respectively. While we discovered a novel region for SBL, a locus for PLH that is close to but distinct from the pubescence density locus which also was found. Unfortunately, because of the complexity of insect resistance mechanism and the LD region for the four pests were still to broad, as dozens to hundreds of genes were found that have potential as candidate resistance genes within the LD region. We highlighted two groups of candidate resistance genes that were previously proved for their functions in insect resistance, but advanced studies using additional approaches such as molecular cloning or differential expression analysis may be needed to reduce the candidate resistance gene pool, as was demonstrated for the Rag2 locus of SBA resistance (Brechenmacher et al., 2015). Our study provided integrated phenotypic and genotypic data and provided a novel insight into soybean insect resistance that might prove useful for insect resistance breeding and reducing insecticide usage.

AUTHOR CONTRIBUTIONS
HC, Completed data analyses and wrote draft of manuscript. GH, Contributed to writing of manuscript.

FUNDING
Research reported in this publication was supported by the USDA Agricultural Research Service.

ACKNOWLEDGMENTS
We thank the contributors that supplied data to the Germplasm Resources Information Network (GRIN) (www.ars-grin.gov) and the SNPs of Soybean Germplasm Collection that maintained by the United States Department of Agriculture Agricultural Research Service (USDA-ARS). We thank T. Herman for her editorial assistance.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017. 00670/full#supplementary-material  Manhattan plot identifies a significant locus on chr 7. Green dots represent the 276 SNPs with FDR below 0.05. Red horizontal line indicates the Bonferroni threshold. (C) Regional LD analysis that surrounds the most significant SNP (ss715596142), shown in red triangle. The pairwise LD between SNPs (gray dots) to ss715596142 follows the right y axis. The orange and gray color for lines indicate SNPs with FDR below or above 0.05, respectively, following the left y axis. The red horizontal dashed line represents the minimal significant level at the cutoff of FDR 0.05. The pink background highlights the high LD region where candidate genes were examined.