Identification of SNPs in RNA-seq data of two cultivars of Glycine max (soybean) differing in drought resistance

The legume Glycine max (soybean) plays an important economic role in the international commodities market, with a world production of almost 260 million tons for the 2009/2010 harvest. The increase in drought events in the last decade has caused production losses in recent harvests. This fact compels us to understand the drought tolerance mechanisms in soybean, taking into account its variability among commercial and developing cultivars. In order to identify single nucleotide polymorphisms (SNPs) in genes up-regulated during drought stress, we evaluated suppression subtractive libraries (SSH) from two contrasting cultivars upon water deprivation: sensitive (BR 16) and tolerant (Embrapa 48). A total of 2,222 soybean genes were up-regulated in both cultivars. Our method identified more than 6,000 SNPs in tolerant and sensitive Brazilian cultivars in those drought stress related genes. Among these SNPs, 165 (in 127 genes) are positioned at soybean chromosome ends, including transcription factors (MYB, WRKY) related to tolerance to abiotic stress.


Introduction
Soybean (Glycine max) is a legume crop that plays an important economic role in the international market, with a world production of almost 260 million tons for the 2009/2010 harvest. Brazil ranks as the world's second largest producer and exporter, with about 25% of the production. Soybean production is influenced by weather oscillations, especially long periods of drought. In the Brazilian soybean culture, the frequency of drought events has increased in recent decades, probably associated with the weather changes in the world (Stokstad, 2004;Schiermeier, 2006). For example, in the states of the south of Brazil, responsible for 40% of domestic production, losses have been as much as 25% of production in recent harvests. During this time some drought tolerant and sensitive cultivars were isolated from these regions. An understanding of the molecular mechanisms governing such contrasting pheno-types upon water deprivation could provide insights for the creation of new cultivars and help in assisted selection.
Single nucleotide polymorphisms (SNPs) are single base differences between DNA sequences of individuals or lines. They can be assayed and exploited as highthroughput molecular markers. SNP markers have the potential for use in association genetics approaches (Cardon and Bell, 2001;Flint-Garcia et al., 2003). The recent availability of high-throughput DNA sequencing data has enabled studies based on highly informative SNPs. The evaluation of SNPs in large EST sequence data sets from agricultural crops has been employed to generate highdensity genetic maps and identify variable genomic regions (Du et al., 2003;Choi et al., 2007;Novaes et al., 2009;Pindo et al., 2008;Duran et al., 2009). The scalability and availability of SNPs in highly automated genotyping assays has made this molecular marker the first choice in genetic linkage and association studies in a variety of species.
( Barbazuk et al., 2007), Illumina Solexa (Van Tassell et al., 2008), and SOLiD (Melum et al., 2010) on different individuals, cultivars or even species are prerequisites for the traditional method of whole genome SNP discovery. In this context, genomic sequences of different individuals are aligned to a reference genome and nucleotide variation is detected.
This study used high-throughput mRNA sequencing (RNA-seq) reads derived from suppression subtractive libraries (SSH) to identify SNPs in up-regulated genes from tolerant (Embrapa 48) and sensitive (BR 16) cultivars submitted to drought stress. Such SNPs can be useful for assisted selection of soybean varieties with higher drought tolerance.

Construction of cDNA libraries and sequencing
Soybean genotype Embrapa 48 (tolerant) and BR 16 (sensitive) were analyzed under two experimental conditions: drought stress and normal irrigation (for further details see Rodrigues et al., 2012, this issue). Leaves and roots from stressed and control plants were collected at three different times (25-50, 75-100, 125-150 min). RNA was isolated from tissue samples and used to construct suppression subtractive hybridization (SSH) cDNA libraries (Rodrigues et al., 2012). These three libraries, enriched in genes up-regulated during drought stress, were sequenced using Illumina/Solexa technology. The reads corresponding to genes enriched in such libraries were used in SNP mapping (see below).

Gene identification from RNA-seq data
The reference transcriptome assembly was constructed using 1,276,813 soybean ESTs available at NCBI. The bdtrimmer software (Baudet and Dias, 2007) was used to exclude ribosomal, vector, low quality and short (less than 100 bp) sequences (using default parameters). The remaining sequences were assembled with the CAP3 program (Huang and Madan, 1999) using default parameters, generating 60,747 unigenes (30,809 contigs and 29,938 singlets) (Nascimento et al., 2012, this issue).
The RNA-seq reads from SSH libraries from tolerant and sensitive cultivars were submitted to quality filtration considering bases greater than Q20 and merged in one single file. The information about the cultivars was included on the read ID to facilitate SNP genotyping. Bowtie software (Langmead et al., 2009), considering default parameters (maximum of 2 mismatches), was used to map the reads against the reference transcriptome and the output file was saved in SAM format ( Figure 1).

SNP detection
SNP detection was performed with the SAMtools pileup program (Li et al., 2009), which found the variations in the SAM file, followed by VarScan (Koboldt et al., 2009), which identifies and filters variants based on read counts, base quality and allele frequency ( Figure 1).
We developed a Perl script to compare the filtered SNP lists generated by the pipeline described above for the two datasets. Putative SNPs were tagged if the reads involved were mapped unambiguously on the reference transcriptome and the minor allele appeared at least 10 times. The SNPs were discarded if the depth was less than 20 and the frequency of one allele was more than 80%. For each candidate SNP, the algorithm accessed the reads over that position in the SAM file, and observed if all variations occurred only between the two cultivars.
The unigenes identified through mapping of RNAseq reads that presented polymorphic sites between both cultivars were annotated and grouped into Gene Ontology classification (GO terms) using blast2go software (Conesa and Götz, 2008). In order to identify the chromosome position of polymorphic genes, the unigenes were mapped to the soybean genome assembly (Schmutz et al., 2010) using Exonerate software (Slater and Birney, 2005).

Results and Discussion
Alignment of RNA-seq reads to reference transcriptome The soybean data available for bioinformatics analysis comprise the reference genome (Willians 82 cultivar; Schmutz et al., 2010), gene models and the EST assembly described in Material and Methods. We chose this assembly as reference to the mapping in order to avoid (1) splicing alignment problems and (2) the absence of the untranslated regions. 332 Vidal et al. After evaluating the SSH libraries, a total of 12,285,871 reads of 45 bp and 30,326,963 reads of 76 bp were obtained for sensitive and tolerant cultivars, respectively. For the sensitive cultivar, 6,317,010 reads were aligned to 7,039 contigs and 2,659 singlets, providing an average depth of coverage of 651 reads by reference sequence. For the tolerant cultivar, 6,120,258 reads were aligned to 15,667 contigs and 6,284 singlets, providing an average depth of coverage of 279 reads by reference sequence. We searched for SNPs in contigs that mapped into both cultivars. After mapping, 7,897 contigs have reads assigned to both cultivars, and were used in SNP identification.

Polymorphism detection
The identification of sequence polymorphisms was performed using SAMtools pileup and VarScan software. We identified a total of 44,510 variations (in 11,000 unigenes) that come from SNPs within each cultivar, SNPs between the cultivars and reference assembly, and SNPs between the cultivars.
To identify inter-cultivar SNPs in the soybean reference sequences and to identify the SNPs between both cultivars (i.e. non-allelic SNPs), we used an in-house SNP filter developed to identify the positions of robust candidate sequence polymorphisms relative to the aligned Solexa reads from each cultivar. We applied another filter requiring at least 10 reads on each cultivar and a maximum of 80% difference between the major and minor allele.
The 6,698 putative polymorphisms were identified in over 2,222 transcripts in tolerant and sensitive cultivars (~3 SNPs/gene). The majority of these polymorphisms represent allelic SNPs (intra-cultivar SNPs) and are not useful as molecular markers for soybean breeding. Nevertheless, we found 165 putative SNPs between tolerant and sensitive cultivars in a total of 127 unigenes (Supplementary Material Table S1). Figure 2 summarizes the GO annotation. As expected, the GO analysis of these 127 genes revealed that many are related to stress response and other ontologies possibly related to stress. Among such genes are a series of signal transducers (calmodulin, ankyrin, GTP binding proteins) and transcription factors (i.e., WRKY, MYB, Zinc Finger, Homeodomain-ZIP), all of which were cataloged in a set of soybean transcription factors (TFs) database focusing on abiotic stress responsive transcription factors (Mochida et al., 2009). Among these soybean TFs, the WRKY and MYB genes were experimentally evaluated. Zhou et al. (2008) verified that soybean WRKY genes provide tolerance to abiotic stresses in transgenic Arabidopsis plants. Soybean MYB genes were considered part of the stress tolerance apparatus based on salt and freezing stress assays in transgenic plants (Liao et al., 2008). The presence of variability in genes that likely coordinate the first steps of stress signaling, thus controlling a series of protective proteins against drought effects, denote such genes as possible markers for assisted selection. We consider all of these polymorphisms as strong candidates for molecular markers for selective breeding.

Mapping SNPs in soybean chromosomes
Putative SNPs detected in uniquely mapped reference sequence unigenes were plotted along the soybean chromosomes. Alignment with the soybean genome showed that the genes with these identified putative SNPs were not distributed uniformly across the genome. Instead, they are more prominent at the chromosome ends ( Figure 3). Wu et al. (2010) detected that many SNPs were clustered in gene-rich, high-recombination euchromatic regions in soybean chromosomes. This positional trait of SNPs may be a consequence of intense recombination/mutation events that are crucial for increased variability in autogamous species such as soybean.

Conclusion
The identification of SNPs in contrasting cultivars for drought stress may be useful to breeders in Marker Assisted Selection (MAS) or even in Genome Wide Selection SNPs contrasting drought resistant cultivars 333  (GWS). They can be added to the upcoming markers derived from high-throughput gene sequencing. We believe that the presence of SNPS in transcription factors is outstanding. Such proteins could be controlling a series of genes responsive to stress, making them good candidates for transgenesis and as starting points to understand drought tolerance mechanisms in soybean.

Supplementary Material
The following online material is available for this article: Table S1 -All genes with intra-cultivar SNP annotations.
This material is available as part of the online article from http://www.scielo.br/gmb.

License information:
This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.