Whole Genome Re-sequencing of Soybean Accession EC241780 Providing Genomic Landscape of Candidate Genes Involved in Rust Resistance.

BACKGROUND: In this study, whole genome re-sequencing of rust resistant soybean genotype EC241780 was performed to understand the genomic landscape involved in the resistance mechanism.
METHODS: A total of 374 million raw reads were obtained with paired-end sequencing performed with Illumina HiSeq 2500 instrument, out of which 287.3 million high quality reads were mapped to Williams 82 reference genome. Comparative sequence analysis of EC241780 with rust susceptible cultivars Williams 82 and JS 335 was performed to identify sequence variation and to prioritise the candidate genes.
RESULTS: Comparative analysis indicates that genotype EC241780 has high sequence similarity with rust resistant genotype PI 200492 and the resistance in EC241780 is conferred by the Rpp1 locus. Based on the sequence variations and functional annotations, three genes Glyma18G51715, Glyma18G51741 and Glyma18G51765 encoding for NBS-LRR family protein were identified as the most prominent candidate for Rpp1 locus.
CONCLUSION: The study provides insights of genome-wide sequence variation more particularly at Rpp1 loci which will help to develop rust resistant soybean cultivars through efficient exploration of the genomic resource.

INTRODUCTION

Soybean, one of the most important oilseed crop worldwide and is severely affected by Soybean rust. Soybean rust is caused by a devastative pathogen Phakopsora pachyrhizi Syd. & P. Syd and spread of the disease occurs through wind disseminated urediniospores. Uredinospores can develop rapidly, causing loss of foliar area and ultimately leading to a severe reduction in seed yield. Soybean rust was reported for the first time in 1902 in Japan, and subsequently in other parts of Asia and Australia during 1934 [1]. Later on, by 2005, soybean rust spread to most of the soybean growing area worldwide, including north and south America [2-4]. Soybean rust is now a global problem that requires urgent attention to secure oil seed production under continuously increasing market demand.

Development of resistant cultivars is considered as most effective approach for the soybean rust control. Resistance to rust pathogen was found in wild species of the genera Glycine [5, 6] and resistance is mediated by a hypersensitive reaction (HR). This is a common type of response triggered when plant resistance genes (R-genes) are challenged by pathogen avirulence genes (Avr-genes) [7]. Resistant genotypes show red–brown (RB) lesions and a spectrum ranging from no to high levels of sporulation depending on the genotype. Susceptible genotypes are characterized mostly by light brown (TAN) lesions with profuse sporulation.

Host plant resistance mediated by Rpp (Resistance to P. pachyrhizi) genes has been found in numerous soybean accessions. About 10 Rpp genes along with different alleles have been mapped to seven genetic loci namely Rpp1, Rpp2, Rpp3, Rpp4, Rpp5, Rpp6 and Rpp7 [8-12]. Garcia et al. also discovered a recessive resistance allele (rpp5) in PI 200456 [13]. Monteros et al. mapped a resistance gene (designated Rpp?Hyuuga) at the Rpp3 locus in the Japanese cultivar Hyuuga (PI 506764) [14], and Kendrick et al. later found a second resistance gene in Hyuuga at the Rpp5 locus [10]. Chakraborty et al. identified Rpp1-b, a different dominant resistance allele at the Rpp1 locus in PI 594538A [15], Ray et al. discovered additional resistance alleles at the Rpp1 locus in PI 587886 and PI 587880A [9], and Hossain et al. also reported other alleles at the Rpp1 locus in PI 594767A and PI 587905 [16].

Recent development in Next-generation Sequencing (NGS) technologies provides opportunities to perform whole genome re-sequencing in a cost effective manner [17, 18]. One of the many applications of whole genome re-sequencing is to understand the genetic basis of phenotypic differences by comparing the genomic sequences. In this regard, in the present investigation, whole genome re-sequencing of rust resistant soybean genotype EC241780 and large scale comparative sequence analysis with susceptible cultivars Williams 82 and JS 335 was performed. The sequences were analysed to discover a large number of DNA polymorphisms including SNPs and InDels in the two parental lines. Further, the functional relevance of SNPs and InDels were analysed based on their presence in intergenic regions, coding and non-coding regions. Functional impact of the SNP on the functionality of the protein was predicted using computational approaches. To our knowledge, this is the first report of whole genome re-sequencing of soybean rust resistant line with sequence depth of ~30X. The results provided here can be used in future studies on high-throughput genotyping and molecular breeding for resistance to rust in soybean.

MATERIALS AND METHODS

Plant Sample Preparation and Sequencing

For whole genome re-sequencing, genomic DNA from the leaves of 22 days old soybean genotype EC241780 was extracted using Qiagen DNeasy kit (Qiagen, GmbH, Hilden). The quality and quantity of the extracted DNA were determined using Bioanalyzer 2100 (Agilent Technologies, Singapore) and Qubit 2.0 Fluorometer (Invitrogen Life Technologies). The library preparation was performed according to the manufacturer’s instructions. Paired-end sequencing of the prepared libraries was performed using Illumina HiSeq 2500 platform (Illumina Technologies) in order to generate 100-bp paired-end reads. Only high quality sequence reads were employed for downstream analysis by trimming low quality reads.

Genomic Resources Used for Comparative Genomics

For comparative analysis, whole genome re-sequenced data of rust susceptible cultivar JS 335 reported by Yadav et al. was obtained from the National Institute of Plant Genome Research (NIPGR), New Delhi, India [19]. In addition, Williams 82 reference genome was also used for the comparative genomics.

Read Mapping

The high quality paired-end reads obtained from whole genome sequencing were aligned to the Williams 82 reference genome (V1.0.26) retrieved from the Ensembl database (http://plants.ensembl.org/) [20]. The sequence alignment was performed using the Burrows-Wheeler Alignment Tool (BWA software available at http://bio-bwa.sourcefourge.net) [21] with default parameters [22]. The mapping output of BWA program was further filtered to retain only those reads that mapped to the reference genome at only one position. Base quality score recalibration and estimates on the coverage of the reference genome were performed using FastQC v0.10.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) [23] and SAMtools v0.1.18 (https://sourceforge.net/projects/samtools/files/samtools/) [24], respectively. The aligned reads were initially sorted using Picard tool and then the read duplicates were removed.

Identification and Functional Annotations of SNPs and InDels

The high-quality reads were used for the identification of SNPs and InDels by individually comparing the EC241780 genome with the Williams 82 reference genomes and JS 335. Variant calling was performed using Genome Analysis TKLite-2.3-9 (GATK) toolkit Unified Genotyper to identify SNPs and InDels. The identified variants were further filtered based on default parameters including depth, and variant score along with slight modifications to retain quality variants. The stringent criteria used for variant filtering were minimum variant frequency of ≥ 90%, average quality of the SNP base ≥ 30 and minimum read depth of 10. Genomic distribution of variant calling was performed by studying the frequency and position of SNPs and InDels for each chromosome.

Annotation of Genomic and Genic SNPs and InDels

Soybean annotation files were obtained from the Ensembl Database (http://plants.ensembl.org/) [20] in order to identify gene positions, functional descriptions, locations of introns, exons and UTRs. SNPs and InDels identified from EC241780 were annotated using data from the Ensemble database and using in-house bioinformatics pipeline (VariMAT – Variation and Mutation Annotation Toolkit). Similarly, functional annotation of the SNPs was performed to identify synonymous and non-synonymous substitutions.

Haplotype Analysis

Publically available 50 thousand SNP data for 19,652 Glycine max accessions [25] were retrieved to perform haplotype analysis. Using the 50K SNP data, a unique haplotype was identified that spans the rust resistance locus identified in EC241780. Phylogenetic analysis of disease resistance genes was performed using MEGA7 software tool [26]. The evolutionary history was inferred using the UPGMA method.

RESULTS

Mapping of Illumina Reads

Whole genome re-sequencing of rust resistant line EC241780 was performed using Illumina HiSeq 2500 platform and quality reads were obtained after FastQ analysis. A total of 374 million raw reads were obtained and adaptor trimming removed 12792 reads (0.003%). High quality reads were mapped to soybean reference genome Williams 82. Overall 287.3 million reads (76.64%) of total reads were mapped to the reference genome. Out of the total aligned reads, 203.5 million reads (70.84%) mapped uniquely to the reference genome (Table ).

Detection and Distribution of Variations

Using minimum default parameters implemented in GATK, initially, we identified a total of 2,082,134 SNPs and 311,361 InDels between the reference Williams 82 genome and soybean accession EC241780. Stringent parameters were further used to minimize the detection of false positive SNPs and InDels. Finally, only those SNPs and InDels with minimum read depth of 10, average base quality of 30 and polymorphism call rate of ≥90% were filtered. Based on these parameters, a total of 2,024,915 polymorphic loci were detected out of which 1,799,343 SNPs and 225,572 InDels were obtained. The distribution of DNA polymorphisms detected for Williams 82/EC241780 across the twenty soybean chromosomes were analysed. The number of SNPs and InDels varied across the chromosome (Fig. ). The total number of SNPs and InDels within a chromosome was found to be directly proportional to the chromosome length. The largest number of SNPs and InDels were detected in chromosome 18 and smallest number in chromosome 12 (Supplementary Table ).

Furthermore, SNPs and InDels were identified by comparing sequence of EC241780 and soybean cultivar JS 335, which is most widely grown soybean cultivar in India and has been commonly used as one of the parent for the development of mapping populations. Initially, 484,704 SNPs and 73,124 InDels were identified between JS 335 and EC241780. Subsequently, with the stringent parameters 414,924 SNPs and 48,962 InDels distributed on 20 chromosomes were identified (Supplementary Table ). Between JS 335 and EC241780 the largest number of SNPs (68,428) and InDels (5008) were detected in chromosome 14 and chromosome 18 respectively, whereas smallest numbers of SNPs (5106) and InDels (1006) were detected in chromosome 13 and in chromosome 4 respectively.

Gene Annotation

The annotation of soybean reference genome was used to reveal the distribution of SNPs and InDels within various genomic regions such as intergenic and intragenic. SNPs and InDels identified from EC241780 were annotated using gene model annotation information for Glycine max present at Ensembl Database. Genomic annotation of variants (SNPs and InDels) in EC241780 revealed that the majority of the variants (>64% of total variants) fall in the intergenic region based on the current gene model. Only 2.9% of SNPs and 1.03% of the InDels were identified within the coding region of the genome. A total of 1,751,039 SNPs and 275,713 InDels were identified in the intergenic region. Genic region comprised of 311,779 SNPs and 65,756 InDels. Within coding region, 75,905 SNPs and 4,347 InDels were identified, whereas, 191,705 SNPs and 46,673 InDels were present in the intronic region. The 5’ UTRs and 3’ UTRs regions also showed the presence of SNPs (44,169) and InDels (14,736) (Supplementary Table ). Functional annotation of SNPs in EC241780 revealed that mostly the variants resulted in missense (53% of total variants) and silent (40% of variants) changes in the coding regions (Supplementary Table ). To validate the variants, over 69 representative SNPs were selected from different positions and were validated using mass array genotyping approach (Supplementary Table ).

Soybean Accession EC241780 Shows High Sequence Similarity With PI 200492 at Rpp1 Locus on Chromosome 18

To identify rust resistance gene, comparative analysis was performed between accession EC241780 and other soybean accessions. Single nucleotide polymorphisms (SNPs) from EC241780 were compared with SoySNP50K data of various rust resistant and susceptible accessions [25]. Interestingly, SNPs from chromosome 18 showed high sequence similarity between accessions EC241780 and PI 200492 near Rpp1 locus (Supplementary Table ). Accession PI 200492 has Rpp1 gene and shows resistance to soybean rust. Phylogenetic tree was further constructed using SNP data from 81 soybean accessions. These 81 accessions include various rust resistant soybean lines which are known source of Rpp1, Rpp2, Rpp3, Rpp4, and Rpp5. All the 81 genotypes grouped into seven different clades (Fig. ). Soybean accession EC241780 grouped with PI 200492 in clade VI along with PI 547875, PI 594177 and PI 368039, which are known to have Rpp1 rust resistance gene. These results clearly indicate that

soybean genotype EC241780 has the same Rpp1 rust resistance gene as in PI 200492. To identify the putative Rpp1 rust resistant gene in EC241780, NBS-LRR genes were identified from over 12 Mb genomic regions from position 48753972 – 61504294 on chromosome 18 and were further characterized (Fig. ).

Comparative Analysis Identified Three Candidate NBS-LRR Disease Resistance Genes at Rpp1 Locus

To characterize the rust resistance genes at Rpp1 locus, we identified nine NBS-LRR genes (Glyma18G46875, Glyma18G51533, Glyma18G51546, Glyma18G51715, Glyma18G51741, Glyma18G51765, Glyma18G51930, Glyma18g51950 and Glyma18G51960) on chromosome 18 between genomic region 56545509 to 6069629 bp (Supplementary Fig. ). Three NBS-LRR genes viz. Glyma18G51715, Glyma18G51741 and Glyma18G51765 are located within the Rpp1 genomic region previously mapped precisely on Gm18 from 60463046 to 60612672 bp [27]. These three genes were shortlisted as potential candidates of rust resistance gene. Multiple sequence alignment of three NBS-LRR proteins showed high protein sequence similarity (Supplementary Fig. ). Comparative analysis was further performed to identify the genetic variation in these NBS-LRR genes from EC2417080 and JS 335. Interestingly, there were larger numbers of SNPs and InDels within the NBS-LRR genes at Rpp1 locus. In soybean accession EC241780, gene Glyma18G51741 has 59 SNPs and one InDel while Glyma18G51765 has 50 SNPs and one InDel in the coding region. Three InDels and 34 SNPs were also identified in Glyma18G51546. Total number of synoyomous and non-synoymous SNPs in Glyma18G51715, Glyma18G51741 and Glyma18G51765 and other NBS-LRR genes near Rpp1 locus are listed in Supplementary Table .

In addition to Rpp1, several studies have demonstrated the presence of different Rpp1 alleles in other soybean accessions. The other soybean lines carrying Rpp1 alleles includes Huuang Dou (PI 587880A), Bai Dou (PI 587886), Xiao Huang Dou (PI 587905), Zhao Ping Hei Dou (PI 594767A), Xiao Jing Huang, Himeshirazu (PI 594177) and PI 561356. Allele Rpp1? in PI 561356 was mapped to the region of Gm18:60518978-60613377 [28]; and Rpp1? in PI 587880A and PI 587886 was mapped to the region of Gm18:58722971-60612672 [15]. We identified five NBS-LRR genes viz. Glyma18G51533, Glyma18G51546, Glyma18G51930, Glyma18g51950 and Glyma18G51960 on chromosome 18 between genomic region 59,283,411 to 60,910,083 which are potential candidates for Rpp1? alleles. Comparative analysis identified several SNPs and InDels in the putative NBS-LRR genes near Rpp1? locus. The largest numbers of 69 SNPs were found in the gene Glyma18g51950 in EC241780 with a total of 18 synonymous and 51 non-synonymous mutations (Supplementary Table ). Gene Glyma18g51546 has 3 InDels with two in-frame mutation and one frame-shift mutation. Large number of SNPs and Indels in the NBS-LRR genes near Rpp1 locus caused significant changes in the NBS-LRR gene sequences in EC241780 as compared to Williams 82 (Supplementary Table ). The functional effect for the amino acid changes located in the candidate genes predicted using PROVEAN tool identified two variations in Glyma18G51741 and five in Glyma18G51765 with deleterious effect (Table , Supplementary Table ). Whereas, none of the variations in gene Glyma18G51715 showed a deleterious effect.

DISCUSSION

This study was initiated with a broad objective of identifying genomic variations between the rust resistant line EC241780 and susceptible cultivar JS 335, which are used as parental lines in soybean breeding programme. Whole genome sequence information of parental lines provides an opportunity to identify informative markers useful for subsequent breeding applications, to prioritise candidate genes at known loci, and most importantly to get allelic information of the candidate genes. Over a million of markers identified with whole genome re-sequencing facilitates high resolution mapping and provides tightly linked or gene specific markers. Apart from the obvious applications, whole genome sequence analysis of resistance and susceptible soybean genotypes performed in the present study helped us to understand the genomic landscape involved in the rust resistance mechanism. Earlier, similar efforts have been made in chickpea, where whole genome re-sequencing of the parental lines have been successfully used to identify candidate genes and to develop genomic resources for the breeding applications [29].

Required sequence depth for whole genome re-sequencing applications largely depends on the biological question. Another major concern is higher cost involved in the high depth sequencing. Earlier, whole genome re-sequencing effort by Zhou et al. achieved an average 11x depth in soybean genotypes [30]. Similarly, Lam et al. performed whole genome re-sequencing of 32 diverse soybean genotypes with about 5x sequence depth [31]. Another significant study by Valliyodan et al. has provided whole genome re-sequencing of 106 soybean lines at 17x depth [32]. All these studies have great importance and provide valuable resource to the soybean community. However, because of low depth the data has limitations in identifying reliable sequence variation in high copy genes and gene family members sharing high level of sequence similarity. In this study, we generated high quality re-sequencing data at 30x depth which is the gold standard as per the present genomics standard. The 30x sequence data helped to identify sequence variation in NBS-LRR gene family with more confidence. We have identified all possible polymorphisms between EC241780 and JS 335 and selected highly informative variations. Our data allowed us to develop genetic markers covering the whole genome which will be highly informative for genetic mapping of rust resistance and may be utilized in a breeding scheme to develop a rust resistant cultivar.

According to the previous studies, soybean lines PI 200492, PI 368039, PI 587886, Xiao Jing Huang, and Himeshirazu may carry the original Rpp1 as identified by Hartwig & Bromfield [33]. Bhor et al. had indicated the presence of Rpp1 gene in genotype EC241780 using SSR markers [34]. In this study, three NBS-LRR genes were identified within the Rpp1 genomic region in soybean accession EC241780. The SNPs, and InDels identified in the putative rust resistance genes will be valuable for genotyping the mapping populations originating from a combination of EC241780 and JS 335 and for marker assisted selection. Further, to facilitate genomics assisted resistance breeding efforts in soybean, chromosome-wise list of useful SNPs have been identified. Thus, we predict that our analyses will be valuable, not only for the fundamental analysis of EC241780/JS 335 derived populations, but also for enhancing our general knowledge of variation in the soybean genome. Identification of molecular markers at the genome wide scale would not only identify markers linked to currently prevalent rust races but also would result in the characterization of novel genomic variants.

The functional annotation of SNPs in NBS-LRR candidate genes showed higher non-synonymous SNPs resulting into amino acid changes. Previous reports in soybean have suggested two models of NBS-LRR gene evolution, first where the genes were tandemly duplicated and then undergone through selective pressures and then copied to different locations [35]. The second model suggest duplication of a single NBS-LRR gene first and then the tandemly duplications as suggested in the first model. Disease resistance genes are known to have a positive selection and more particularly NBS-LRR gene family [36-38]. In Arabidopsis, Mondragón-Palomino et al. have reported strong positive selection in many members of the NBS-LRR gene family [39].

SSR and SNP markers play a great role in the construction of high density linkage map in crop species. The availability of a large number of SSRs, SNPs and InDels from parental lines increases the utility of such markers in genomics assisted breeding. Hence, the development of novel genomic resources for mapping rust-resistant loci in soybean is imperative. This study offers high-density coverage of DNA polymorphisms across the entire soybean genome, which could be utilized for high-resolution genotyping and in molecular breeding programs. In this study, we obtained high quality homozygous SNPs from JS 335 and EC241780. Such genomic variation will greatly help in enhancing molecular markers repository, fine mapping investigations, and molecular genetics research for soybean rust.

Recently, Pedley et al. have identified genes governing rust resistance through extensive VIGS mediated screening [40]. This study identified Rpp1 and Rpp1b on the same chromosomal region suggesting the complexity of the locus with possible multiple resistance genes with high allelic variation. Almost a decade before Pedley et al. [40], Chakraborty et al. have shown the co-localization of Rpp1 and Rpp1b using QTL mapping approach [15]. Several similar examples where reported candidate genes does not have high confidence or miss identified due to the localization of different resistance genes have been reported in soybean [41]. A high number of alleles, the drastic difference between allelic responses, and co-localization of resistance genes makes it difficult to select the appropriate candidate. The identification of the correct candidate is indeed highly important for molecular understanding but does not limit the exploration of loci in the breeding program.

CONCLUSION

The sequence variations identified in the present study not only provide the promising candidates but also facilitate marker development for breeding application. The information generated in the present study pave the path for the exploration of Rpp1 locus for the development of rust resistance soybean cultivar.

Table 1

Summary of sequence reads derived from rust resistance soybean genotypes EC241780 which were mapped onto the reference genome (Williams 82).

S. No.	Mapping details of EC241780	-
1.	Total number of quality reads	374,907,382
2.	Number of aligned reads	287,318,932
3.	Percentage of aligned reads	76.64
4.	Number of uniquely aligned reads	203,558,177
5.	Percentage of uniquely aligned reads	70.84

Table 2

Deleterious effects predicted by PROVEAN tool for the amino acid changes located in the candidate genes for rust resistance in soybean selected based on whole genome sequence information of genotype EC241780. Detailed result for each of the predicted amino acid change is provided in the Supplementary Table 6.

Variant#	Score	Prediction (Cutoff= -2.5)	-	Variant	Score	Prediction (Cutoff= -2.5)
Glyma18G51741
D543E	-2.583	Deleterious	-	E543D	2.616	Neutral
I522F	-2.558	Deleterious	-	F522I	2.591	Neutral
Glyma18G51765
A1163E	4.293	Neutral	-	E1163A	-4.193	Deleterious
N767H	5.053	Neutral	-	H767N	-5.453	Deleterious
W593C	8.376	Neutral	-	C593W	-8.409	Deleterious
I520T	-3.556	Deleterious	-	T520I	3.572	Neutral
D517E	-2.599	Deleterious	-	E517D	2.716	Neutral

# None of the variation observed in gene Glyma18G51715 found to have Deleterious effect.