pea ( pisum sativum l. ) , an important crop cultivated worldwide , is a valuable model system in plant genetics . since gregor mendel 's famous experiments , several scientific discoveries have occurred in modern pea genetics ; these new insights include information regarding genetic control of compound leaf development [ 2 , 3 ] and the molecular basis of symbiotic interactions with beneficial nitrogen - fixing bacteria ( rhizobia ) [ 46 ] .
the study of pea gene polymorphism in relation to agronomically important traits is essential to both basic and applied research on this crop plant [ 7 , 8 ] .
unfortunately , the large size ( more than 4  gb ) ( http://data.kew.org/cvalues/ ) and complexity   of the pea genome hamper its sequencing as well as the discovery of this crop plant 's genetic resources , both of which are desperately needed for molecular and genomics - assisted breeding [ 8 , 10 ] . as an alternative to
whole genome sequencing , analysis of transcriptomes by rna sequencing can provide extensive information about expressed genes [ 11 , 12 ] .
because next - generation sequencing technologies are applicable to all organisms , including those for which information about genome organization is insufficient or lacking , considerable progress in pea transcriptome sequencing has been achieved over the last few years .
massive amounts of transcriptomic data have been obtained in the form of high - quality sequence reads that have been used for molecular marker creation , whole genome map construction [ 1317 ] , and characterization of host - pathogen interactions ( pea - sclerotinia sclerotiorum ) .
all these data ( as well as additional unpublished pea transcriptome sequencing results ) have been uploaded to the ncbi sequence read archives ( sra ) ( http://www.ncbi.nlm.nih.gov/sra/ ) as raw reads .
assemblies have been created for some of these data and deposited in the ncbi transcriptome shotgun assembly ( tsa ) database , allowing users to perform data mining ( e.g. , blast searching ) and to study pea gene polymorphism .
several genes have tissue - specific expression , however , and can therefore only be studied through analysis of the appropriate tissue .
one such example involves symbiotic genes necessary for the establishment and development of nitrogen - fixing nodules which are predominantly expressed in those temporary plant organs . to date
, only a few samples from pea nodules have been sequenced ( available as raw sra archives ) , and only one assembly built from a mixture of sequencing reads from different organs ( including nodules ) is present in the tsa database ( see table 1 for available pea nodule transcriptome sequencing results ) .
because this assembly was based on nodules harvested at a very late stage of symbiotic nodule development ( 3-month - old plants ) , it presumably contains insufficient information on nodule - specific transcripts .
the aim of our work was thus to sequence the transcriptome of young pea nodules , construct an assembly , and analyze the resulting assembly for unique sequences .
along with nodules , we harvested root tips to analyze their transcriptome content as well .
seeds of pea laboratory line sge   were surface - sterilized with concentrated sulfuric acid ( 98% ) ( 15  min on a shaker ) , washed 10 times with autoclaved distilled water , and germinated on petri dishes containing sterile vermiculite for 3 days .
the germinated seeds were then planted individually into 200  ml ceramic pots containing quartz sand , watered with 100  ml of 2x nitrogen - free mineral nutrition solution , and inoculated with an aqueous suspension of rhizobium leguminosarum bv .
plants were harvested 12 days after inoculation ; nodules and root tips ( 5  mm distal portion of the root ) were placed in liquid nitrogen , ground into powder , and stored at 80c .
total rna was extracted from 100  g of material using an rneasy plant mini kit ( qiagen , hilden , germany ) .
cdna libraries were constructed and sequenced according to the instructions provided with the genome analyzer iix platform ( illumina , san diego , ca , usa ) .
after total rna extraction and dnase - i treatment , mrnas were captured using oligo ( dt ) magnetic beads and fragmented .
first - strand cdna was synthesized from these fragments using random hexamer primers ; double - stranded cdna was then generated , purified with magnetic beads , and subjected to end reparation and 3 single adenylation .
sequencing adaptors were ligated to the adenylated fragments , and dna fragments having adapter molecules on both ends were then amplified .
after a quality control step performed on a 2100 bioanalyzer ( agilent technologies , santa clara , ca , usa ) , the cdna library products were sequenced in a single - read run with 75  bp length reads on an illumina genome analyzer iix platform .
preliminary quality control of the raw sequencing data was performed via the fastqc v.0.11.3 application ( http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ ) , which indicated that the reads were of acceptable quality . for adapter removing cutadapt version 1.8.1  
low - quality read removal and trimming were then performed with the assembly program trinity v.2.0.6   using the 
next , contig assembly was performed by trinity with default assembly parameters , including kmer = 25 . as a result ,
two fasta files were obtained , one for the nodule sample and one for the root tip sample .
statistical parameters for the assemblies were obtained by running the trinitystats.pl script included in the trinity package . as a step in assessment of assembly quality
, the generated reads were mapped to the assemblies with the bowtie2 program v. 2.2.5 .
the contigs were grouped with known pea sequences ( obtained from http://www.ncbi.nlm.nih.gov/ ) using cd - hit_est from the cd - hit package ( http://cd-hit.org/ )   with parameters -c 0.80 -n 6 . in order to distinguish the transcripts enriched in nodules as compared to root tips , the reads of both libraries
the differential expression was calculated using edger package   under a negative binomial model , with biological coefficient of variation 0.2 and fdr cutoff value 0.001 .
the resulting blast output was processed using publicly available blast2go software ( v.2.5.0 ) ( biobam bioinformatics sl , valencia , spain )   to retrieve associated gene ontology ( go ) terms describing biological processes , molecular functions , and cellular components . to detect transcripts containing reliable full - length cds regions
two approaches based on similarities of either nucleotide or amino acid sequences were used . for each transcript , blast search against ncbi refseqgene database was performed in order to find orthologous sequences , and then these sequences were aligned by smith - waterman algorithm   with a  5 - 0  substitution matrix . also , as an alternative approach , we used transdecoder software   for cds region prediction based on homology search against swiss - prot protein database . to extend the annotation of the full - length nodule - specific transcripts , the nucleotide sequences were converted into amino acid sequences and then mapped to the kyoto encyclopedia of genes and genomes ( kegg ) web server ( http://www.genome.jp/kegg/ ) .
direct sequencing of pcr fragments was performed on an abi prism 3500  xl system ( applied biosystems , usa ) at the genomic technologies , proteomics , and cell biology core center of all - russia research institute for agricultural microbiology ( arriam , saint petersburg , russia ) .
sequencing generated 52,021,865 reads from the  nodules  library and 17,684,604 reads from the  root tips  library .
after removal of adapter and index sequences , 75.1% reads were equal to or longer than 70  bp and 10.2% of reads were shorter than 10  bp for  nodules  and 74.5% and 11.2% of reads , accordingly , for  root tips . 
a total of 58,397 contigs belonging to 48,628 genes ( as termed by trinity ) were constructed from the nodule set , with a mean contig length of 880.81 , a median contig length of 620 , and an n50 of 1,282 .
root tip reads were assembled into 37,287 contigs of 35,081 genes , with a mean contig length of 841.14 , a median contig length of 558 , an n50 of 1,260 , and 1,055 total isoforms .
length distributions of contigs in the two assemblies are presented in figure 1 , and distributions of isoform numbers are shown in figure 2 . to our knowledge , the appearance of isoforms can be due to either alternative splicing or the presence of paralogous sequences expressed in the tissue .
as an illustration of the first case , we were able to detect two splice variants in the  nodules 
assembly for transcripts of the symbiotic gene ign1 , an ortholog of ign1 ( ineffective greenish nodules 1 ) of lotus japonicus ( regel . )
k. larsen   ( genbank accession number kr047192 ; tr2831|c0_g2_i1 and tr2831|c0_g2_i4 in the  nodules  assembly ) .
the longer transcript retained the first intron and , according to prediction by netstart 1.0 , could be translated into a protein variant lacking the first 21 n - terminal amino acids .
an example of the second case involved enod6 ( for early nodulin 6 ; genbank accession number x63700 ) , which encodes a short protein belonging to a group of nodule - specific cysteine - rich ( ncr ) peptides [ 37 , 38 ] .
assembly ) derived from paralogous genes encoding cysteine - rich peptides specific for nodules ( it should be noted that some of these isoforms could be artificial chimeras containing parts of different transcripts that have extensive segments sharing 100% similarity ) . evaluating the quality of a de novo transcriptome assembly without a reference genome is challenging .
we therefore implemented three approaches previously recommended for managing this task [ 39 , 40 ] .
first , we used the pea expressed - sequence tag ( est ) sequences represented in genbank as a standard to estimate assembly quality .
of these , 2,571 ests ( 13.8% ) shared no similarity with any contigs of the assembly .
furthermore , 102 ests were filtered out on the basis of an e - value cutoff of 1  10 . from the remaining 15,903 ests ( 85.6% ) , we chose hits with maximal coverage of the est ( one per est ) and evaluated their coverage and identity distributions . among these est - contig pairs ,
following the second recommended approach , we mapped back all reads to contigs in both assemblies ; as a result , 89% and 91% of reads in 
nodules  and  root tips  assemblies were , respectively , aligned back to the contigs , demonstrating that our assemblies were of acceptable quality .
efd ( ethylene response factor required for nodule differentiation ) , vpy ( vapyrin ) , and nsp1 ( nodulation signaling pathway 1 )   of medicago truncatula gaertn . and sen1 ( stationary endosymbiont nodule 1 )   of l. japonicus .
long transcripts with high identity were found for all four genes ( table 2 ) ; these transcripts allowed us to design primers flanking coding sequence ( cds ) regions and to amplify the corresponding regions in cdna synthesized from 4-week - old nodules of pea genotypes sge and finale . except for allelic variations of vpy , sen1 , and nsp1 that were found between sge transcriptome and finale cdna genotypes , we observed complete sequence correspondence for all four genes , thereby demonstrating the satisfactory quality of the created assembly . to evaluate the  root tips  assembly , we selected genes involved in glutathione biosynthesis : gsh1 ( gamma - glutamylcysteine synthetase precursor [ af128455.1 ] ) , gshs ( glutathione synthetase precursor [ af231137.1 ] ) , and hgshs ( putative homoglutathione synthetase [ af258319.1 ] )
assembly identified one contig ( tr9283|c0_g1_i1 ) completely identical to full - length gsh1 and one contig ( tr8244|c0_g1_i1 ) completely identical to full - length gshs except for a single nucleotide polymorphism in the 3 untranslated region .
this search also revealed four contigs ( tr2244|c0_g1_i1 , tr3920|c0_g1_i1 , tr13401|c0_g1_i1 , and tr25033|c0_g1_i1 ) representing portions of hgshs ( each with 100% identity ) that had not been assembled into a contig , probably because of insufficient overlapping of reads due to the low expression levels of this gene in pea root tips .
consequently , despite the good quality of the  root tips  assembly , its coverage was insufficient for finding full sequences of rare transcripts ; nevertheless , the discovery of partial sequences allows primers to be designed for whole - transcript pcr amplification and transcript - end amplification by rapid amplification of cdna ends ( race ) methodology . to obtain information on nodule - specific genes
, we attempted to select portions of previously unknown sequences of the  nodules  assembly by clustering them together with pea sequences produced from nonnodular tissues .
in addition to our  root tips  assembly , the pea transcriptome assemblies created by franssen et al .   and
kaur et al .   currently meet this requirement . using cd - hit software , 277,211 sequences of these four pea transcriptome assemblies
were grouped into 61,521 clusters ( where a cluster is defined according to cd - hit as a set of similar sequences created to reduce sequence redundancy and to improve the performance of other sequence analyses ) . among these clusters , 10,391
the 13,305 nodule - specific clusters included 14,998 contigs belonging to 14,171 genes ( i.e. , without gene isoforms ) .
the transcripts were first aligned against plant protein sequences in the ncbi nr protein database ( 24.04.15 release ) .
the following parameters were used : an e - value cutoff of 1  10 , the same alignment direction for all high - scoring segment pairs ( hsps ) in a hit , and 20 ( or more , if of the same e - value ) hits for a query .
only about 30% of sequences less than 1,000  bp long in our  nodules  assembly were successfully annotated .
one - third of nodule - specific contigs ( 5,940 ) were associated with plant proteins in the nr database , with 565,464 total hits . of these ,
 nitrogen compound metabolic process  ) along with  single - organism cellular process  , 
,  establishment of localization  , and  single - organism localization  ( figure 4(a ) ) .
this distribution reflects the processes occurring in nodules , such as microsymbiont ( rhizobia ) hosting within cells and nitrogen compound metabolism . within the molecular function category ,
contigs were assigned to  heterocyclic compound binding  ,  organic cyclic compound binding  ,  ion binding  ,  small molecule binding  ,  transferase activity  , and  carbohydrate derivative binding  ( figure 4(b ) ) .
these terms may be related to metabolite exchange between plants and bacteria , including exchanges with signal molecules .
regarding cellular components , the major go terms were  cell part  ,  membrane - bounded organelle 
( figure 4(c ) ) , which are similarly concerned with the formation and functioning of symbiotic compartments in nodule cells .
it seems also valuable to distinguish the transcripts that are preferentially expressed in nodules as compared to root tips . by mapping the reads of both libraries (  nodules  and  root tips  ) to  nodules  assembly via bowtie2 v. 2.2.5 and calculating the differential expression via edger package with 0.001 fdr cutoff we selected 1081 contigs that represent genes with significantly higher expression level in nodules ( supplementary file 1 in supplementary material available online at http://dx.doi.org/10.1155/2015/695947 ) .
still , more detailed analysis , including verification of the results of such  digital expression  analysis by real - time pcr , is needed , along with addition of more time points to the experiment .
to identify sequences of novel unreported , highly reliable transcripts of pea , we analyzed the 14,998 sequences of the  nodule only  clusters using transdecoder software .
as a result , 593 putative full - length orfs were found in 536 contigs ( supplementary file 2 ) . as an alternative approach to the identification of full - length transcripts
the remaining sequences were filtered according to the following criteria : minimal query coverage of 0.8 , maximum e - value of 1  10 , and the same direction for all hsps .
we aligned each pair using the smith - waterman algorithm   with a  5 - 0  substitution matrix and identified aligned fragments corresponding to the cds regions of hit sequences ( as determined by genbank ) .
we selected 427 alignments ( comprising 153 unique pea sequences , some of which aligned to multiple genbank accessions representing gene isoforms or paralogs ) with the following characteristics : ( 1 ) full coverage of the hit cds region by a contig ; ( 2 ) identity higher than 0.8 ; and ( 3 ) the contig having possible start and stop codons within a 50  bp region .
, we identified 581 novel sequences containing putative full - length cds in pea nodule transcriptome , among which 536 were found by transdecoder and additional 45 were detected by alternative approach based on blast against known plant mrna sequences .
for annotation of these 581 sequences , homologous genes were found by blastn search in medicago truncatula genome ( ver .
ko ( kegg orthology ) identifiers were assigned to the novel sequences , and 109 entries out of 581 ( 18.8% ) were successfully annotated ( supplementary file 2 ) . in our opinion , our generated
assembly adds valuable information , especially with respect to nodule - specific sequences , to the existing knowledge about pea transcriptome : some unique sequences of pea symbiosis - related genes can be identified only in our assembly .
an example of this case involves cle genes , some of which were shown to participate in systemic regulation of nodule formation in several legumes such as m. truncatula , l. japonicus , and glycine max [ 4648 ] .
the clavata3/embryo surrounding region - related ( cle ) gene family is composed of numerous genes that contain conserved cle domains in various plant species and encode short regulatory peptides ( cle - peptides ) ( for review see ) . in m. truncatula ,
two cle genes , mtcle12 and mtcle13 , have nodulation - related expression patterns that are linked to proliferation and differentiation . in pea , sequences of cle genes are not known , but it was shown that overexpression of mtcle13 gene leads to similar effects ( severe reduction in nodulation ) in both pea and m. truncatula , proving that mtcle13 is functional in pea .
so we sought for the sequences homologous to mtcle12 and mtcle13 ( medtr4g079630.1 and medtr4g079610.1 , resp . ) in our nodule transcriptome assembly .
blastn search using the medtr4g079610.1 transcript sequence ( encoding cle13 peptide ) as a query against our  nodules  assembly retrieved the contig tr8317|c0_g1_i1 , which contains a full open reading frame ( orf ) corresponding to cle13 of p. sativum ; the same search against  organism pisum sativum ( taxid : 3880 )  in the ncbi tsa database returned two partial transcripts : ( 1 ) gb|gcmk01019899.1| ( tsa :  pisum sativum ps_029064 transcribed rna sequence  ) , containing only part of the cle13 orf , and ( 2 ) the apparently chimeric gb|gcmo01040960.1| ( tsa :  pisum sativum ps_150017 transcribed rna sequence  ) containing a portion of the ankyrin repeat gene ( similar to glycine max ankyrin repeat - containing protein at5g02620-like [ loc100812799 ] ) as well as the 5-part of the cle13 transcript .
a blastn search using the medtr4g079630.1 transcript sequence ( encoding cle12 peptide ) as a query found two contigs : tr116|c0_g1_i1 , containing the full orf of p. sativum cle12 , and tr23484|c0_g1_i1 , containing the full orf of an unknown protein similar to p. sativum and m. truncatula cle12 and therefore presumably a paralog of p. sativum cle12 .
we thus tentatively designated these genes as cle12a and cle12b , respectively ( table 3 ) . at the same time , the search against the ncbi tsa database retrieved no significant homologs of cle12 in pea .
the pea transcriptome assembled after illumina sequencing is thus a good resource for the study of pea transcripts related to nodulation .
such potential research targets include genes encoding nodule - specific peptides such as ncr- and glycine - rich protein peptides that have been exhaustively described in m. truncatula [ 38 , 51 ] but not p. sativum , as well as other symbiotic genes expressed in nodules , including cle peptide - encoding genes . also , the present
assembly is a convenient tool that can facilitate study of transcription changes in nodules of symbiotic mutants . in pea ,
several mutant lines with impaired nodule formation were obtained and phenotypically characterized [ 5254 ] .
rna sequencing of the whole transcriptome from mutant nodules is considered to be a reasonable approach for further characterization of genes and gene networks that operate during nodule development . in this regard
, the present assembly of pea nodule transcriptome can be used as a reference for mapping reads and differential expression analysis . also , this reference transcriptome is indispensable for annotation of short contigs obtained according to the mace ( massive analysis of cdna ends ) protocol , which implies sequencing of 3-part of each transcript instead of the whole mrna .
direct annotation of 3-parts of transcripts by gene ontology or kegg is often not successful because of dissimilarity of these regions between different species , and the present assembly containing significant number of pea nodule - specific transcripts can therefore serve as a reference for annotation of differentially expressed transcripts revealed by mace technology .
the aim of the present study was the acquisition of nodule - specific transcript sequences via next - generation sequencing . using an illumina platform
, we obtained 52 million reads from a sample derived from young pea nodules and more than 17 million reads from a root tip sample .
we constructed the assemblies ( more than 58,000 and 37,000 contigs from nodules and root tips , resp . ) and analyzed the  nodules 
we identified approximately 15,000 nodule - specific contigs associated with different go biological function terms . of these
, 581 sequences were found to possess full cdss and could thus be considered as new nodule - specific transcripts of pea . because the ability of pea plants to form symbiotic nodules is an agronomically important trait , information about pea nodule - specific gene sequences can be applied by scientists and breeders for primer design , gene - based marker creation , polymorphism studies , and real - time pcr .
the next challenge for researchers is characterization of pea transcripts specific to another symbiosis formed by leguminous plants : arbuscular mycorrhiza , which is also of great importance to both fundamental science and contemporary sustainable agriculture .
 the large size and complexity of the garden pea ( pisum sativum l. ) genome hamper its sequencing and the discovery of pea gene resources . 
 although transcriptome sequencing provides extensive information about expressed genes , some tissue - specific transcripts can only be identified from particular organs under appropriate conditions . in this study , we performed rna sequencing of polyadenylated transcripts from young pea nodules and root tips on an illumina gaiix system , followed by de novo transcriptome assembly using the trinity program . 
 we obtained more than 58,000 and 37,000 contigs from  
 nodules  and  root tips  assemblies , respectively . 
 the quality of the assemblies was assessed by comparison with pea expressed sequence tags and transcriptome sequencing project data available from ncbi website . 
 the  nodules  
 assembly was compared with the  root tips  assembly and with pea transcriptome sequencing data from projects indicating tissue specificity . as a result , approximately 13,000 nodule - specific contigs were found and annotated by alignment to known plant protein - coding sequences and by gene ontology searching . of these 
 , 581 sequences were found to possess full cdss and could thus be considered as novel nodule - specific transcripts of pea . 
 the information about pea nodule - specific gene sequences can be applied for gene - based markers creation , polymorphism studies , and real - time pcr . 
1. Introduction
2. Materials and Methods
3. Results and Discussion
4. Conclusions