UBC Faculty Research and Publications

Transcriptome profiling of wheat glumes in wild emmer, hulled landraces and modern cultivars Zou, Hongda; Tzarfati, Raanan; Hübner, Sariel; Krugman, Tamar; Fahima, Tzion; Abbo, Shahal; Saranga, Yehoshua; Korol, Abraham B Oct 13, 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


52383-12864_2015_Article_1996.pdf [ 1.5MB ]
JSON: 52383-1.0307439.json
JSON-LD: 52383-1.0307439-ld.json
RDF/XML (Pretty): 52383-1.0307439-rdf.xml
RDF/JSON: 52383-1.0307439-rdf.json
Turtle: 52383-1.0307439-turtle.txt
N-Triples: 52383-1.0307439-rdf-ntriples.txt
Original Record: 52383-1.0307439-source.json
Full Text

Full Text

RESEARCH ARTICLE Open AccessTranscriptome profiling of wheat glumes inwild emmer, hulled landraces and moderncultivarsHongda Zou1†, Raanan Tzarfati1†, Sariel Hübner2, Tamar Krugman1, Tzion Fahima1, Shahal Abbo3,Yehoshua Saranga3 and Abraham B. Korol1*AbstractBackground: Wheat domestication is considered as one of the most important events in the development ofhuman civilization. Wheat spikelets have undergone significant changes during evolution under domestication,resulting in soft glumes and larger kernels that are released easily upon threshing. Our main goal was to explorechanges in transcriptome expression in glumes that accompanied wheat evolution under domestication.Methods: A total of six tetraploid wheat accessions were selected for transcriptome profiling based on their rachisbrittleness and glumes toughness. RNA pools from glumes of the central spikelet at heading time were used toconstruct cDNA libraries for sequencing. The trimmed reads from each library were separately aligned to the referencesub-genomes A and B, which were extracted from wheat survey sequence. Differentially expression analysis andfunctional annotation were performed between wild and domesticated wheat, to identity candidate genes associatedwith evolution under domestication. Selected candidate genes were validated using real time PCR.Results: Transcriptome profiles of wild emmer wheat, wheat landraces, and wheat cultivars were compared using nextgeneration sequencing (RNA-seq). We have found a total of 194,893 transcripts, of which 73,150 were shared betweenwild, landraces, and cultivars. From 781 differentially expressed genes (DEGs), 336 were down-regulated and 445 wereup-regulated in the domesticated compared to wild wheat genotypes. Gene Ontology (GO) annotation assigned 293DEGs (37.5 %) to GO term groups, of which 134 (17.1 %) were down-regulated and 159 (20.4 %) up-regulated in thedomesticated wheat. Some of the down-regulated DEGs in domesticated wheat are related to the biosyntheticpathways that eventually define the mechanical strength of the glumes, such as cell wall, lignin, pectin and waxbiosynthesis. The reduction in gene expression of such genes, may explain the softness of the glumes in thedomesticated forms. In addition, we have identified genes involved in nutrient remobilization that may affect grain sizeand other agronomic traits evolved under domestication.Conclusions: The comparison of RNA-seq profiles between glumes of wheat groups differing in glumes toughnessand rachis brittleness revealed a few DEGs that may be involved in glumes toughness and nutrient remobilization.These genes may be involved in processes of wheat improvement under domestication.Keywords: Tetraploid wheat, Domestication, Glumes, RNA-seq, Differentially expressed genes* Correspondence: korol@research.haifa.ac.il†Equal contributors1Department of Evolutionary and Environmental Biology, The Institute ofEvolution, Faculty of Natural Sciences, University of Haifa, Haifa 31905, IsraelFull list of author information is available at the end of the article© 2015 Zou et al. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link tothe Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.Zou et al. BMC Genomics  (2015) 16:777 DOI 10.1186/s12864-015-1996-0BackgroundDomestication of plants was a major event in the estab-lishment of agriculture and human civilization. Wheatwas among the first domesticated plant species and isconsidered as one of the most important crops in theworld. Comparative studies of domesticated wheat withits wild progenitors lead to insights about the geneticbasis of their adaptation which could be beneficial forfuture crop improvement. During domestication andsubsequent crop improvement under domestication, nu-merous morphological and physiological characteristicsof the wild progenitors were modified to meet humanneeds. The first and pristine domestication trait in ce-reals, non-brittle rachis, is related to the loss of kerneldispersal mechanisms. As a result, there was a transitionfrom shattering hulled forms of wild einkorn wheat (T.boeoticum L., AbAb) and wild emmer wheat (T. turgidumL. ssp. dicoccoides, AuAuBB, also known as T. dicoc-coides), to non-shattering hulled (as hard-threshing)forms in the diploid einkorn wheat (T. monococcum L.,AmAm) and tetraploid emmer wheat (T. turgidum L. ssp.dicoccum, AuAuBB), respectively. Later on, during theevolution under domestication, a variety of changes haveoccurred, related to the glumes toughness, proportion ofkernel weight in the whole spike weight, shape andcolour, seed dormancy, disease and pest resistance, andhigh productivity in a wide range of environment [1].The genome of tetraploid wheat originated about 0.5million years ago from an interspecific hybridizationevent between the T. urartu (AuAu) and an unknown Bgenome ancestor presumably related to Aegilops spel-toides. The genome of hexaploid wheat has resultedfrom a second inter-specific hybridization between do-mesticated tetraploid cultivated emmer T. dicoccum(AuAuBB) and Ae. tauschii (DD) followed by genome du-plication ~9,000 years ago [2]. Durum wheat (T. turgidumL. ssp. durum) is the predominant form that was selectedfrom emmer and has free-threshing grain. Thus, T. dicoc-coides is the progenitor of both durum and bread wheat,and is central to wheat domestication evolution [3, 4].The genetic basis of events involved in plant domesti-cation and the nature of selection in domesticated cropshave been subjected to intense molecular genetics andgenomics studies over the past two decades [5, 6]. Alarge number of wheat domestication-related genes havebeen identified through quantitative trait locus (QTL)mapping [7–11], genome-wide association studies [12],and cloning [13, 14]. QTL mapping was one of the majorapproaches in genetic studies of plant domestication evo-lution and improvement, as well as in unravelling theagronomic potential of their wild progenitors. Most QTLanalyses of wheat domestication and improvement fo-cused on spike traits, including brittle rachis (preventingseed shattering) [8, 15] and glumes toughness (ease ofthreshing) [9, 16]. Many QTL studies have demonstratedthat major key domestication traits are controlled by arelatively small proportion of the genome, implying thateither pleiotropy or tight linkage among several loci maybe an important attribute in the evolution of domesticatedcrops [8, 11, 17]. Nowadays, dense SNP genetic maps areavailable for the traditional QTL analysis of populationsderived from crosses of domesticated plants with theirwild progenitors [18] as well as for the genome-wide asso-ciation studies [19, 20]. Comparison of QTL map loca-tions with genome sequencing or genome-wide SNPscanning has also been used to identify candidate genomicregions involved in selection during domestication[21, 22]. Cavanagh et al. [6] developed a high-throughputarray to integrate 9 K gene-associated SNPs in a world-wide sample of 2994 accessions of hexaploid wheat in-cluding landraces and modern cultivars to characterizethe impact of crop improvement on genomic and geo-graphic patterns of genetic diversity. The results showedthat there are minor genetic differences between landracesand cultivars. In another study, a wheat genotyping arraywas developed with about 90 K gene-associated SNPs,which is an excellent resource for fine-scale genetic dissec-tion of domestication related traits [23].Additional attempts to illuminate the domesticationprocess by using functional genomics included expressedsequence tag (EST) sequencing, microarray and more re-cently, RNA-seq technologies. Ergen and Budak con-structed six subtractive cDNA libraries and sequencedover 13,000 ESTs using wild emmer wheat accessionsand modern wheat in order to analyse the expressionprofile of drought related genes [24]. The first micro-array comparison between developing spikes of tetra-ploid wild (T. dicoccoides) and domesticated wheat (T.dicoccum and T. durum) at the stage of one week afterpollination, identified 38 and 24 differentially up- ordown-regulated genes, respectively, out of 2493 cDNAclones on the array [25]. Most of the genes that werefound to be up-regulated in the domesticated wheatwere related to carbon metabolism, such as Rubiscolarge and small subunits and the sucrose synthase.Among down-regulated genes in domesticated wheat theauthors noted storage protein genes and genes associ-ated with abiotic and biotic stress responses. Althoughcomprehensive studies using the microarray hadachieved a better understanding of the wheat genomeexpression [26, 27], the microarray technology has somelimitations compared to RNA-seq. Microarray analysisrelies on hybridization between probes and targets. Mostmicroarray studies are based on commercial arrays suchas the Wheat Genome Array (Affymetrix), where targettranscripts were designed using EST libraries of culti-vated wheat. Nevertheless, since there is high sequencesimilarity between wild and cultivated wheat, it was alsoZou et al. BMC Genomics  (2015) 16:777 Page 2 of 14successfully used for expression studies of wild emmer[28–30]. Nowadays, the advanced technology of NextGeneration Sequencing (NGS), enabling to sequence thewhole transcriptome (RNA-seq), was proved as an excel-lent approach to study changes in domestication relatedgenes and expression networks underlying plant domesti-cation and crop improvement [31–33]. NGS has remark-able advantages over the microarray in the detection ofnovel transcripts, allele-specific expression and splicejunctions [34]. Hence, RNA-seq can expand our view andprovide new insights into plant domestication evolution atthe genomics level.Wheat glumes are an important part of the spikelet,which is the dispersal unit of the plant. Genes involvedin development and structure of the glumes and spikesare interesting from both theoretical and practical as-pects [35]. The glumes are the closest vegetative tissueto the grain. As part of their role in reproduction, theglumes serve as a ‘defense line' for the kernels, and acton nutrient allocation and photo-assimilates conversa-tion destined for the developing kernels [36]. Theglumes composition and structure can greatly impactplants performance and their interaction with environ-ment. Recently, it was suggested that glumes can serveas a photosynthetically active sinks adjusting for thechanging metabolism demand of the kernels [37].Glumes can also maintain their metabolic activity longerthan other vegetative organs and influence the final yieldand nitrogen cycling [38]. Moreover, there is indicationthat glume phenotype has a possible correlation withsome beneficial agronomic traits [39]. Genes affectingglumes, like Q in wheat and tga1 in maize, were involvedin key steps of domestication and are related to diversebiological functions, implying significant roles of theglumes [13, 40]. As noted above, wheat glumes haveundergone significant changes along evolution under do-mestication. The main outcome of this process was thereduction in glumes toughness and the increase of thekernels weight proportion in the total spike weight(SpHI, spike harvest index) [16].In the current study, we explored the evolutionarychanges of the tetraploid wheat transcriptome by com-parative RNA-seq analysis of three dissimilar genotypicgroups, wild emmer wheat, tetraploid landraces andmodern T. durum cultivars, representing three differenttime points in wheat domestication. We have identifiedlarge differences in gene expression between the wildand domesticated wheat. Among the differentiallyexpressed genes, we identified genes that may be in-volved in glumes toughness and threshability, nutrientremobilization and the proportion of kernels in thewhole spike weight and other agronomic traits evolvedunder domestication.MethodsPlant materialA total of six tetraploid wheat accessions were selected fortranscriptome profiling based on their rachis brittlenessand glumes toughness [16]. These included: (1) two wildemmer wheat T. dicoccoides (accessions Y12-3 and A24-39) characterized by brittle rachis and tough (hulled)glumes; (2) traditional landraces including T. dicoccum(G581) and T. ispahanicum (G805) characterized by non-brittle rachis and tough glumes; and (3) two modern culti-vars of T. durum (‘Inbar’ and ‘Svevo’), characterized bynon-brittle rachis and soft glumes (Table 1).Plants were grown in three biological replicates as de-scribed in [16]. Glumes of the central spikelet of each geno-type were sampled at its heading time (when the spike wasfully emerged). Each accession was sampled independently1 h after sunrise. Glumes were collected, placed immedi-ately in Eppendorf tubes with RNAlater (Qiagen, Hilden,Germany), and stored at −20 °C for RNA extraction.RNA extraction and sequencingRNA was extracted from glumes using the Plant MiniKit including a digestion step with DNase I (Qiagen,Standford, CA, USA) for removal of DNA traces. Highquality RNA was confirmed using Bioanalyzer 2100 withRNA 6000 Nano Labchips (Agilent, Santa Clara, CA,USA). RNA samples were pooled to three groups in ac-cordance with their level of domestication, i.e., wild,landraces and cultivars. As the main objective of thisstudy was to identify transcription differences along do-mestication “gradient”, pooling samples should givehigher credence to representative genes of groups ratherthan genotypes. Each of the pools contained 1 μg RNAof the two accessions (Table 1). For each RNA pool, twoindependent biological replicates (i.e., six pools) wereused to construct RNA-seq libraries, and a third repli-cate was reserved for QPCR validation. The cDNA li-braries were constructed using NEBNext UltraDirectional RNA Prep Kit (New England Biolabs, MA),following the manufacturer’s instructions. After verifyingTable 1 Wild, landrace and cultivar tetraploid wheat genotypesused in the studyGroup Species Accessions Rachis and glumescharacterizationWild T. turgidum L. subsp.dicoccoidesY12- 3 Brittle rachis, hard tothreshA24-39Landrace T. turgidum L. subsp.dicoccumG 805 Non-brittle rachis, hardto threshT. ispahanicum Heslot G 581Cultivar T. turgidum subsp.durum (Destf.)Inbar Non-brittle rachis, softglumesSvevoZou et al. BMC Genomics  (2015) 16:777 Page 3 of 14the quality of the libraries indexed with six-nucleotidebarcodes, sequencing was performed on the IlluminaHiseq2000 machine using multiplexing for generating2 × 101 bp paired end reads. Sequencing was carried outat the Technion Genome Center (Haifa, Israel).Data processing, mapping and SNPs discoveryA tetraploid reference genome was prepared in silico byextracting sequences assigned to the A and B genomesfrom the chromosome survey sequencing (CSS) data of theIWGSC (International Wheat Genome Sequencing Con-sortium, http://www.wheatgenome.org) [41]. Sequencesfrom each RNA-Seq pool were cleaned and trimmed by re-moving adaptor sequences and low-quality reads usingTrimmomatic software (version 0.32) [42] with the follow-ing parameters: phred64, LEADING: 3, TRAILING: 3, SLI-DINGWINDOW:4:20, MINLEN:40 (phred quality scoresQ ≥ 20, read length ≥ 40). Each cleaned library was alignedto each of the tetraploid reference subgenomes separately,using the Subjunc aligner in Subread package (version1.43) [43] with the following parameters: -d 0, -D1000, -u, -H, -I 16, -S fr. The -u option was used to reportuniquely mapped reads only, whereas -H option was usedto breaks ties using Hamming distance when there wasmore than one best mapping location for a read, whichwould give the most accurate mapping results with littleor no cost to the mapping percentage.Because it is not feasible yet to index a large genome(more than 4 Gbp) by Subread, we had to split the wheatreference genome AABB into sub-genome A and sub-genome B, and then combine the alignment results usingthe following method. After alignment, the sum of map-ping quality scores (MQS) for each mapped read wasused to determine to which sub-genome (A or B) theread should be assigned. For accurate alignment, theread pairs had priority over singletons (when only oneread of a pair was mapped) and uniquely mapped readshave priority over ambiguously mapped reads. When thesame read was mapped to the two genomes, the genomewith the higher MQS was accepted and the other onewas discarded. The read that had the same alignmentscore in the two genomes was discarded by the customscript (such reads comprised a very low percentage).This may be an applicative methodology whenever thegenome size exceeds the tools limitation that can helpus to further characterize homoeolog-specific reads.Genotype calling was carried out with the alignment filesusing SAMtools/BCFtools (version 0.1.19, http://samtools.sourceforge.net) with default parameters. All SNPs withmaximum read depth less than 100 were kept for subse-quent analysis. The relationship between the mapping ra-tios and genetic distance from reads to reference genomewere examined by Pearson correlation.Differential gene expression analysisWe further used featureCounts [44] in the Subread pack-age to quantify the level of expression for each gene basedon the associated gtf (Gene Transfer Format) file providedwith the survey sequence. In order to reveal differentiallyexpressed genes (DEGs) in domesticated vs. wild acces-sions, we considered the common part of two subsets:DEGs between cultivated vs. wild and DEGs betweenlandrace vs. wild. DEGs in each of these two comparisonswere identified using DESeq software (version 1.6.1) [45]at selection cutoff log2Foldchange ≥ 1 and 10 % FDR(False Discovery Rate), implying that p-values were ad-justed for multiple testing based on Benjamini-Hochbergapproach at a level below 0.1.Functional analysis of differentially expressed genesGene Ontology (GO) terms were searched with Blas-t2GO [46]. First, we extracted the sequences of theDEGs from reference genomes and gtf files with acustom script. Then the sequences of DEGs werecompared to the NCBI nr (non-redundant) databaseusing blastx with a cutoff e-value less than 1e-5 [47].The blastx output, generated in xml format, was usedfor Blas2GO analysis to annotate the DEGs. GO func-tional classification for DEGs was performed usingthe WEGO software [48].Unmapped reads processing, de novo assembly,differential gene expression analysis and functionalannotationReads that failed in the alignment procedure were ex-tracted from alignment files using SAMtools (version0.1.19); and assembled de novo in Trinity (version 2.06)[49] with default parameters. We further aligned the rawunmapped reads of each group to the assembled contigswith Bowtie [50] and estimated genes abundance usingRSEM [51]. Because most of the reads identified as un-mapped were essentially ambiguously mapped, assem-bled transcripts with high identity (>70 %) to theIWGSC reference genome found by blastn were dis-carded. The remaining transcripts were annotated withBlast2GO as mentioned above and used to identify dif-ferentially expressed transcripts (DETs) among the threegroups using edgeR with default parameters [52].Quantitative reverse transcription PCR (QRT-PCR)A few transcripts were selected for validation of RNA-seqresults by QRT-PCR as described in [29, 30]. Total RNA(1 μg) from relative samples were used to generate com-plementary DNA (cDNA) templates using the qScriptcDNA synthesis kit (Quanta BioSciences, Gaithersburg,MD, USA). Gene-specific primers and SYBR Green PCRmaster mix (Quanta BioSciences, Gaithersburg, MD,USA) were used for QPCR on StepOne System (AppliedZou et al. BMC Genomics  (2015) 16:777 Page 4 of 14Biosystems, Darmstadt, Germany), according to the manu-facturer’s instructions. Gene-specific primers were de-signed by Primer3 (http://primer3.ut.ee/primer3) (Table 2).PCRs included 2.5 μL of 1st-strand cDNA (1:4 dilutedcDNA), 7.5 μL of SYBR Mix and 300 nM of each primer ina final volume of 10 μL. The primers were designed innon-polymorphic regions across accessions. We tested theefficiency of the primers using four serial dilutions (of 1:4)for each of the wheat groups. The specificity of each primerpair was monitored by heat dissociation curve analysis ofthe amplicon as a final step of the PCR. For COBRA,CesA-1, and CCR genes we used cDNA samples of threebiological replicates per accession as a template for QPCR.For the other six genes, FLA, FST, CesA-2, 6-SFT, LAC andLAC16, we used the cDNA pools similarly to the RNA-seq pools rather than individual cDNA (due to lim-ited RNA with two to three technical replicates for eachbiological sample). The endogenous gene Ubiquitin wasused to normalise variations within well-to-well andacross plates with forward and reverse primers: 5’ -TTGACAACGTGAAGGCGAAG - 3’ and 5’ -GGCAAAG ATGAGACGCTGCT - 3’ respectively.Quantification of gene expression was carried out bythe relative quantification method (2-ΔΔCT method) [53]implemented in StepOne software v2.3. Since efficien-cies of the several primers were not identical in thethree wheat group, we corrected the relative quantifica-tion results according to the efficiency of each of thegenes in each of the groups by StepOne software v2.3,which is based on: Efficiency = 10 (−1/slope) for each pri-mer and RQ = Etarget- (Ct treatment - Ct control)/Endogenous-(Ct treatment - Ct control).ResultsAssembly of reads into transcripts of the A and BgenomesSequencing of mRNA from glume tissue of wild, landraceand modern tetraploid wheat pools generated 147.4 mil-lion pair-end reads of 101 bp length. The number of readsfrom each pool ranged from 17.7 to 29.6 million (Table 3).Table 2 Primers for QRT-PCRGene Gene ID Forward (5’ – 3’) Reverse (5’ – 3’)CCR Ta1alLoc003924.2 TGCCGTGAGAAGAAGGTGAT CTTCTCTGCCATCGTCTTGCCOBRA Ta5asLoc003744.1 CGTCGCCGTTGAAGTAGATCT TCATCGCAAGGATGATAGAACCesA-1 Ta5alLoc000723.1 GTCAAGCAAGAACAGCAT ACCAGCTACCCACAAGAGCAAFLA Ta3bLoc003710.1 CAGTACCCGCTCAACGTCAC CCGGTGTAGAGCGTGTTGTCFST Ta3bLoc056384.1 TGAACATGAGCAAGCTGGAG CCATTCTGTCTCCCATGTCCLAC16 Ta4asLoc013789.1 GTCCGATCTACCCGTCTGTT GTGTGTTTATGCAAACCAAAGGLAC Ta4blLoc021918.2 GCAGAAGGTGACACGGCTAT GCCTTCCCTTGTGACGATTCesA-2 Ta5alLoc000723.1 ATCAGGCTTTGATTTCAGCAA TGGATGTCATGTCAAGCAAGA6-SFT Ta6bsLoc005412.1 AGCTGTCAGTGAGGGTGCTT CGTTGGGTACACTCGTGATGTable 3 Summary of samples and RNA-seq dataGroups Total reads Clean reads Sub-genome MappedreadsMappingratio (%)Total mappingratio (%)Mappedreads -uMapping ratio(%) -uTotal mappingratio (%) -uWild-1 23,271,884 22,088,185 A 17,586,172 79.6 86.5 13,102,552 59.3 74.6B 17,813,886 80.6 13,140,929 59.5Wild-2 27,577,198 26,160,151 A 20,700,875 79.1 85.9 15,242,936 58.3 73.7B 21,023,950 80.4 15,242,734 58.3Cultivar-1 17,740,180 16,773,599 A 13,527,727 80.6 87.4 9,474,682 56.5 71.3B 13,783,371 82.2 9,535,891 56.9Cultivar-2 29,551,798 27,977,895 A 22,256,822 79.6 86.6 16,395,908 58.6 74.3B 22,663,043 81.0 16,290,987 58.2Landrace-1 23,082,003 21,898,663 A 16,563,967 75.6 85.3 14,141,135 64.6 83.1B 16,971,148 77.5 14,275,683 65.2Landrace-2 26,163,934 24,724,657 A 18,580,284 75.1 84.4 15,687,853 63.5 81.9B 19,196,718 77.6 15,843,276 64.1Total 147,386,997 139,623,150Zou et al. BMC Genomics  (2015) 16:777 Page 5 of 14After removal of ambiguous nucleotides, low-quality se-quences (phred quality scores <20), and adapter se-quences, a total of 139.6 million cleaned reads were usedto quantify the level of expression by mapping to thetetraploid reference sequences extracted from the wheatchromosome survey sequencing (CSS) data. Withoutthe -u flag in Subjunc, the mapping ratios of sequencescorresponding to the A and B genomes of wheat rangedbetween 75.1 and 79.6 % in the A genome and from 77.5to 82.2 % in the B genome. With the -u flag, these ratiosranged from 56.5 to 64.6 % in the A genome and from56.9 to 65.2 % in the B genome. After assigning each readto the corresponding genome based on mapping qualityscore, the mapping ratio over all remaining reads reacheda range of 84.4–87.4 % without -u option and 71.3–83.1 %with -u option. Therefore, we used the results obtained byusing the -u flag to maintain a conservative approach fordownstream analysis.To test whether genetic similarity between each pooland the reference sequences has an effect on mapping ra-tio we compared the genetic distance calculated from vari-ant called among pools and mapping ratio. No correlationwas found between genetic similarity and mapping ratioeither with (r = 0.12, p = 0.701) or without (r = 0.41,p = 0.182) the -u flag (Additional file 1: Figure S1).These results further support our analytical approachand corroborate the downstream expression results.Differentially expressed genes (DEGs) betweendomesticated and wild wheatA total of 194,893 transcripts were found expressed inwild, landrace and cultivar pools, out of which 73,150transcripts were commonly expressed in the threegroups (Fig. 1). Differential expression was first com-pared between the wild wheat and each of the twogroups of domesticated wheat (landraces or modern cul-tivars) and then between wild and domesticated (land-race + modern cultivars). We found a higher number ofDEGs (2193 DEGs) in the comparison of modern vs.wild wheat than in the comparison between landracesvs. wild wheat (1662 DEGs). In the modern cultivar vs.wild, 1035 DEGs were down-regulated and 1158 DEGswere up-regulated in modern cultivars as compared tothe wild progenitor (Additional file 2: Table S1). A totalof 746 DEGs were down-regulated and 916 DEGs wereup-regulated in landraces as compared with the wildprogenitor (Additional file 3: Table S2). The comparisonbetween the domesticated (landraces + modern cultivars)vs. wild accessions identified only 781 DEGs, of which 445genes had higher expression in the domesticated and 336DEGs had higher expression in the wild wheat (Figs. 2 and3; Additional file 4: Table S3). A heat-map of 781 signifi-cant DEGs between wild and domesticated pools was cre-ated using DESeq (Fig. 4), demonstrating clustering thatdistinguished between the domesticated (landraces andmodern cultivars) and wild pools.Functional analysis of DEGs between wild anddomesticated wheatIn order to investigate transcriptome changes in glumesevolution under domestication, we assessed the expres-sion patterns of the DEGs in domesticated (landraces +cultivars) vs. wild wheat (Additional file 4: Table S3). Toannotate the DEGs in wild and domesticated groups,sequences were searched against the NCBI non-redundant (nr) protein database by blastx using a cut-off e-value of 10−5. GO terms were subsequentlyassigned to DEGs based on the blastx results. Out of781 DEGs, only 293 DEGs (37.5 %) were assigned toGO-term groups, including 134 (17.1 %) DEGs down-regulated and 159 (20.4 %) DEGs up-regulated in thedomesticated wheats compared to the wild accessions(Fig. 5). The DEGs were categorized into 29 groupsbased on GO annotation. The categories ‘cell, cellparts and organelles’, ‘binding and catalytic’ and ‘cellu-lar process and metabolic process’ showed highestnumbers of GO terms for the ‘cellular component’,‘molecular functionality’ and ‘biological process’ cat-egories, respectively. Interestingly, structural moleculeand transcriptional regulator (in ‘molecular’ GO cat-egory) and growth (in ‘biological process’ GOcategory) were found only among the DEGs down-Fig. 1 Proportional Venn diagram of transcripts among wild wheatgenotypes, cultivars and landracesZou et al. BMC Genomics  (2015) 16:777 Page 6 of 14regulated in the domesticated genotypes compared towild wheats. Molecular transducer and transitionregulator in the ‘molecular function’ GO categorywere found only in DEGs up-regulated in the domes-ticated compared to wild wheats.Based on gene annotation, we selected DEGs thatcould be regarded as candidate genes for domesticationprocess in two possible directions of evolutionarychanges (i.e., up- or down-regulated in the domesticatedwheat compared to its progenitor). However, we selectedthe candidate genes based on uni-direction expressiondifference (either up- or down-regulated in the domesti-cated accessions). For example, three cellulose synthasegenes were found only in the down-regulated DEGs inthe domesticated wheat, whereas three amino acid per-mease genes were found only among the up-regulatedDEGs in the domesticated wheat. A total of 22 DEGswere down-regulated in domesticated compared to wildwheat (Table 4). Using the available gene annotation[41], we found that many of these DEGs are related tocell wall organization or biogenesis; phenylpropanoidmetabolism; and carbohydrate metabolism and transporta-tion. For example, genes encoding for cellulose synthase(CesA), fasciclin-like arabinogalactan (FLA), trichomebirefringence-like (TBL), fiber protein, pectin lyase-likeprotein and pectin acetylesterase family protein, and CER1are related to cell wall organization and biogenesis. Genesencoding for phenylalanine ammonia lyase (PAL), cinna-moyl CoA reductase (CCR), flavonol 4-sulfotransferase(FST) and 4-coumarate:CoA ligase (4CL) are involved inphenylpropanoid metabolism and lignin biosynthesis; whilegenes encoding for sucrose synthase 2 (SUS2) and sucro-se:fructan-6-fructosyltransferase (6-SFT) are responsiblefor carbohydrate metabolism and transportation. We fur-ther identified transcription factor genes down-regulated indomesticated wheat, such as genes encoding plant-specifictranscriptional regulator NAC domain protein. Notably,COBRA genes have shown significant down-regulation inmodern cultivar compared to wild wheat (see Table S1 andDiscussion).In contrast to above down-regulated genes under do-mestication, 17 DEGs were found to be significantly up-regulated in domesticated wheat compared to the wildprogenitor (Table 5). The most abundant groups of up-regulated DEGs in the domesticated pools included 3-ketoacyl-CoA synthase (KCS) gene and Chalcone synthase(CHS) gene. Besides that, we also found amino acid per-mease (AAP) gene and silicon transporter (SIT) gene wereup-regulated in domesticated wheat.Unmapped reads processing, de novo assembly,differentially expression analysis and functionalannotationThe unmapped reads extracted from six libraries werepooled together and de novo assembled using Trinity togenerate a set of transcrips absent from the referencegenome. From the unmapped reads, 64,316 contigs wereassembled with length ranging from 224 bp to 24,492 bpand N50 of 494 bp. After removing transcripts that hadhigh identity (>70 %) to the IWGSC reference genome,7264 contigs ranging between 224 bp and 4296 bp werekept. These contigs are considered as novel transcripts.A total of 2761 novel transcripts had significant hit in1387 856806Cultivar vs. Wild Landrace vs. Wild822 410336 590 471445Cultivar < Wild Landrace < Wild Cultivar > Wild Landrace > WildFig. 2 Proportional Venn diagrams of DEGs in domesticated compared to wild wheat. a Total DEGs. b DEGs down-regulated in domesticatedwheat. c DEGs up-regulated in domesticated wheatFig. 3 Histograms of DEGs in cultivar, landrace and domesticatedcompared to wild wheat. a DEGs in cultivar genotypes. b DEGs inlandraces. c DEGs in domesticated wheatZou et al. BMC Genomics  (2015) 16:777 Page 7 of 14searches against the nr database using blastx with cutoff1e-5. GO analysis was conducted by Blast2go and GOterms were assigned to 1622 transcripts (Additionalfile 5: Figure S2). Differentially expressed transcriptswere validated based on the protocol for downstreamanalyses of de novo assemble using Trinity (see sec-tion Methods). We found 110 DETs in modern vs. wildwheat, out of which 67 were down-regulated and 43were up-regulated in modern cultivars as compared tothe wild progenitor. We also found 111 DETs in land-race vs. wild wheat, out of which 68 were down-regulated and 43 up-regulated in landrace as comparedto the wild progenitor. The comparison between the do-mesticated vs. wild accessions detected only 59 DETs, ofwhich 52 had higher expression in the domesticated and7 had higher expression in the wild wheat. It should benoticed that the overwhelming majority of these DETshave no known function and missing information abouttheir sub-genome location (Additional file 6: Table S4).Quantitative reverse transcription PCR (QRT-PCR)To validate the RNA-seq results, QPCR was performedfor nine selected DEGs that appeared to have higher ex-pression in the glumes of the wild accessions. Moreover,based on gene annotation, these genes can be consideredas interesting candidates for glumes evolution under do-mestication. PCR products of these genes’ primers amp-lified for each of the wheat groups (wild, landraces andcultivars) were specifically indicated by the single-peakmelting curves. Due to the limited RNA, we tested theW1 W2 C1 C2 L1 L2Fig. 4 Heat map of DEGs in glumes of domesticated vs. wild wheat.The heat map represents the genes expression level of the 781significant DEGs between wild and domesticated wheat (cultivar pluslandrace) from all the six groups (log2Foldchange ≥ 1 and FDR ≤ 0.1).Blue color indicates gene expression level. Wild is abbreviated to W,Cultivar is abbreviated to C, and Landrace is abbreviated to LFig. 5 Comparison of Gene Ontology classifications of DEGs in domesticated vs. wild wheat. Blue color indicates down-regulated DEGs in domesticatedcompared to wild wheat, red color indicates up-regulated DEGs in domesticated compared to wild wheat. All of DEGs are categorized into 29 functionalgroups based on GO classificationZou et al. BMC Genomics  (2015) 16:777 Page 8 of 14Table 4 DEGs down-regulated in glumes of domesticated wheat compared to wild progenitorid baseMean baseMean baseMean log2Fold Padj log2Fold Padj Putative annotation(Wild) (Cultiviar) (Landrace) Change(C/W) (C/W) Change(L/W) (L/W)Ta7alLoc001275.1 1093.35 416.13 388.49 −1.39 7.96E-03 −1.49 3.19E-03 4-coumarate:CoA ligaseTa1blLoc007155.1 5137.83 1403.42 2048.36 −1.87 1.00E-05 −1.33 6.17E-03 Cellulose synthaseTa5alLoc000723.1 1933.05 441.04 647.98 −2.13 7.97E-07 −1.58 6.55E-04 Cellulose synthaseTa3bLoc028980.1 4678.25 1223.81 1929.83 −1.93 3.84E-04 −1.28 8.50E-02 Cellulose synthaseTa4alLoc006547.1 177.41 0.00 0.00 -Inf 2.47E-15 -Inf 1.30E-17 CER1 proteinTa4alLoc026069.1 550.15 0.00 0.00 -Inf 8.29E-26 -Inf 4.65E-29 CER1 proteinTa1alLoc003924.2 40.98 1.82 0.39 −4.49 5.99E-02 −6.71 1.80E-03 Cinnamoyl CoA reductaseTa3bLoc003710.1 5610.65 1515.59 1849.21 −1.89 0.0001 −1.60 0.0016 Fasciclin-like arabinogalactan protein 7Ta3bLoc056384.1 25.40 0.00 0.00 -Inf 6.75E-02 -Inf 3.53E-02 Flavonol 4-sulfotransferaseTa4alLoc006913.1 299.00 4.85 7.32 −5.95 1.64E-16 −5.35 1.08E-16 Flavonol 4-sulfotransferaseTa5blLoc013288.1 944.48 181.67 269.42 −2.38 1.04E-07 −1.81 1.66E-04 NAC domain-containing protein 18Ta6blLoc001596.1 1507.23 377.78 671.34 −2.00 8.20E-06 −1.17 5.04E-02 Pectin lyase-like proteinTa3bLoc036242.1 251.80 34.64 32.65 −2.86 1.16E-03 −2.95 1.53E-04 Pectinacetylesterase family proteinTa3bLoc019897.1 219.61 36.13 30.72 −2.60 1.40E-02 −2.84 4.83E-03 Pectinacetylesterase family proteinTa2blLoc014498.1 719.00 118.04 306.22 −2.61 4.62E-08 −1.23 6.04E-02 Phenylalanine ammonia-lyaseTa3bLoc024051.1 342.10 34.22 106.23 −3.32 4.41E-09 −1.69 1.19E-02 Phenylalanine ammonia-lyaseTa7asLoc021287.1 3214.99 932.09 1145.19 −1.79 3.10E-02 −1.49 8.21E-02 Sucrose synthase 2, putative, expressedTa6bsLoc005412.1 839.45 149.96 376.09 −2.48 3.87E-07 −1.16 8.73E-02 Sucrose:fructan-6-fructosyltransferaseTa4alLoc019947.1 168.96 45.07 38.20 −1.91 1.72E-02 −2.14 5.10E-03 Fiber protein Fb34Ta7bsLoc005648.1 161.29 39.24 43.95 −2.04 2.28E-02 −1.88 3.37E-02 TRICHOME BIREFRIGENE like 22Ta4blLoc021918.2 235.61 19.18 49.36 −3.62 4.92E-08 −2.26 8.61E-04 LaccaseTa4asLoc013789.1 939.60 192.99 304.20 −2.28 5.70E-07 −1.63 1.16E-03 laccase 16 LENGTH=523Table 5 DEGs highly up-regulated in glumes of domesticated wheat compared to wild progenitorid baseMean baseMean baseMean log2Fold Padj log2Fold Padj Putative annotation(Wild) (Cultiviar) (Landrace) Change(C/W) (C/W) Change(L/W) (L/W)Ta7asLoc021951.1 42.47 1189.91 717.98 4.81 1.80E-04 4.08 2.50E-03 3-ketoacyl- synthase 12-likeTa7bsLoc002749.1 2.79 306.18 109.26 6.78 5.12E-05 5.29 4.28E-03 3-ketoacyl- synthase 12-likeTa6asLoc018551.1 0.00 53.86 18.54 Inf 4.91E-04 Inf 2.85E-02 3-ketoacyl-CoA synthaseTa7asLoc001384.1 20.78 591.48 521.17 4.83 1.59E-02 4.65 2.39E-02 3-ketoacyl-CoA synthaseTa7bsLoc001848.1 14.15 363.40 330.26 4.68 5.37E-02 4.54 8.58E-02 3-ketoacyl-CoA synthaseTa7bsLoc002750.1 14.74 443.97 186.17 4.91 3.14E-03 3.66 7.64E-02 3-ketoacyl-CoA synthaseTa7bsLoc005172.1 0.46 89.24 55.43 7.59 1.12E-02 6.90 3.76E-02 3-ketoacyl-CoA synthaseTa7bsLoc011512.1 67.99 919.81 753.44 3.76 3.74E-03 3.47 1.19E-02 3-ketoacyl-CoA synthaseTa6bsLoc008917.1 0.50 195.60 217.19 8.62 4.50E-10 8.77 3.18E-10 Chalcone synthaseTa2bsLoc009111.1 0.46 1663.14 453.70 11.81 2.19E-16 9.93 5.54E-12 Chalcone synthase 8, putativeTa3bLoc000987.1 0.93 1583.66 2159.97 10.74 4.24E-22 11.18 2.26E-23 Chalcone synthase 8, putativeTa4bsLoc019947.2 0.00 51.71 97.56 Inf 8.90E-04 Inf 1.74E-04 Chalcone synthase 8, putativeTa6bsLoc002330.1 0.00 306.31 386.67 Inf 3.33E-19 Inf 2.60E-20 Chalcone synthase 8, putativeTa2alLoc009166.1 71.52 237.64 320.20 1.73 0.0077 2.16 6.86E-05 Amino acid permease 6Ta2alLoc010251.2 11.86 139.49 88.59 3.56 0.00 2.90 6.08E-03 Amino acid permease-like proteinTa5asLoc003267.1 28.62 147.37 107.71 2.36 0.00 1.91 2.31E-02 Amino acid permeaseTa1alLoc016727.3 36.90 268.40 132.25 2.86 2.30E-07 1.84 1.04E-02 Silicon transporterZou et al. BMC Genomics  (2015) 16:777 Page 9 of 14efficiency of the primers of gene FLA, FST, CesA-2, 6-SFT, LAC and LAC16, using four serial dilutions (of 1:4)for each of the wheat groups. Amplification efficienciesbased on slopes of standard curve showed that theseprimers had high efficiencies ranged from 98.7 to 107.0 %,except FST gene amplified in wild pool (84.7 %) and 6-SFTgene amplified in cultivar pool (113.3 %). Meanwhile, theR2 values ranged from 0.994 to 0.999, except CesA-2 amp-lified in landrace pool (0.794) (Additional file 7: Table S5).As shown in Fig. 6, the fold-changes in gene expressiondetermined by QPCR were quite consistent with theirnormalized read counts (expression level) determined byRNA-seq. Namely, all nine genes selected for QRT-PCRvalidation exhibited a considerable reduction in their ex-pression level from wild to domesticated groups.DiscussionPlant domestication has fascinated scientists interestedin the evolutionary process ever since Darwin. Primaryefforts were aimed to discover the wild progenitors ofdomesticated plants using classical taxonomy and genet-ics. Subsequently, phylogenetic distances between wildand domesticated plants were established by DNAmarkers including RFLP, SSR, AFLP, DArT and SNP[54]. Ayal et al. [25] were the first to address the ques-tions related to wheat domestication by studying alter-ations in the transcriptome using cDNA microarray.They found 63 up- or down-regulated genes betweenwild and domesticated wheat. With the development ofNGS technology, there was tremendous progress in theevolutionary studies aimed at unravelling the molecularbasis of domestication using RNA-seq that can detect ex-pression changes in thousands of genes. To the best ofour knowledge, our study is the first that used RNA-seq tocompare domesticated and wild tetraploid wheat glumes.The transition from brittle rachis to non-brittle rachiswas probably the first (pristine) domestication event.After the domestication episode, wheat glumes weresubject to selection that made them more suitable forhuman needs. Some of the consequences were the emer-gence of easier to thresh spikes, which have a lower per-centage of chaff, i.e., an increased proportion of the totalkernel weight in the spike weight compared to the wildwheat. The wild and the landrace accessions of tetra-ploid wheat selected for this study have tough glumesand hulled seeds, which are non-free threshing. In con-trast, the modern cultivars are free threshing (have softglumes and non-hulled seeds). In our previous study re-lated to threshing time, the three studied groups showeda pattern of gradual decrease, consistent with thechronological time frame from wild to landrace andfrom landraces to modern cultivars [16]. To some ex-tent, the noted phenotypic difference could be caused bythe observed lower expression level of genes related tothe cell wall composition and glumes toughness (e.g.,genes in the lignin biosynthesis pathway including PAL,4CL and CCR) in the domesticated genotypes. Further-more, there was a significant increase in the SpHI inlandraces compared to the wild wheat accessions and aslight improvement in the assayed modern cultivarscompared to the landraces. This increase in the SpHIcould be a consequence of the finer glumes and up-regulation of genes involved in the transport of aminoacids (e.g., amino acid permease), which can facilitate inN retranslocation and grain filling [55].We selected hulled-glume wild and landrace acces-sions for comparison with free-threshing modernFig. 6 QRT-PCR validation of RNA-seq results for DEGsZou et al. BMC Genomics  (2015) 16:777 Page 10 of 14cultivars, in order to search for DEGs that may be asso-ciated with evolution under domestication. Since thewheat genome is not completely sequenced yet, we usedthe wheat survey sequence [41] that provides the infor-mation needed for phasing homeologs of the A and Bgenomes. Until now, there is no reliable draft genomesequence in tetraploid wheat. However, the sequences ofchromosome 5B, which is the first genomic survey se-quence in wild emmer wheat, has been published re-cently [56] Our results detected 123,370 transcripts inthe cultivar pools which are slightly lower than in thepublished Triticum turgidum transcriptome (140,118)built by the de novo assembly method [57]. The corres-pondence between the two studies is very good, despitethe fact that we analysed only glumes at heading timewhile Krasileva et al. [57] analysed young roots, youngshoots, spikes and grains. The possibility that there maybe less expressed transcripts in glumes than in other or-gans is consistent with a previous RNA-seq study of dif-ferent tissues in barley [55].The identified DEGs may be sought as genes that wereeither preferred or rejected not by the early farmers anddue to their association with traits were subject to selec-tion efforts during improvement evolution under domesti-cation. Yet, the possibility that some changes inexpression patterns was a result of correlated responses toselection caused by tight linkage or linkage disequilibriumof corresponding genes with agriculturally beneficial al-leles rather than directional selection should not be over-looked. Since there is a correlation between glumes shapeand some agronomic traits [39], it could be speculatedthat at some time point(s) during evolution under domes-tication, the shape of glumes served as an indication/marker for the presence or absence of specific traits ofinterest.Candidate genes for wheat evolution underdomesticationTo understand changes in gene expression that occurredduring evolution under domestication of tetraploidwheat, we selected 39 candidate DEGs in glumes for fur-ther characterization. Of these genes, 22 DEGs hadlower expression in domesticated wheat; some are re-lated to cell wall organization or biogenesis. In general,major components of plant cell wall are cellulose, hemi-cellulose, pectin, lignin and protein. However, we are notaware of other studies of genome expression in theglumes in the context of wheat domestication. Amongthe 22 down regulated DEGs we identified the followingcell wall related genes: CesA genes are responsible forcellulose synthesis, and evolved in primary and second-ary cell wall development of wheat [58]. FLA, a subset ofarabinogalactan protein (AGP), has both an AGP-likeglycosylated region and a putative fasciclin domain,which may contribute to cell adhesion, communicationand cell wall architecture in Arabidopsis, rice and wheat[59, 60]. TBL is a protein family containing a plant-specific DUF231 domain and may be involved in biosyn-thesis and deposition of secondary wall cellulose inArabidopsis [61]. Pectin is also one of the most import-ant components of the primary cell wall in plants. Wealso found DEGs related to pectin metabolism, such asgenes encoding pectin lyase-like protein and pectin acet-ylesterase family protein, which were down-regulated indomesticated compared to wild wheats. The lignin isconsidered as a major component of the secondary cellwall, providing the strength in plants. We have identifieda series of DEGs in the pathway of lignin biosynthesis,including PAL, CCR, FST and 4CL, which is in agree-ment with previous studies of cotton [32]. Likewise, twogenes encoding for laccases (LAC), which may be in-volved in lignin polymerization [62], were also down-regulated in domesticated wheat. All these genes weredown-regulated in the glume of domesticated wheat,suggesting that cell wall synthesis in glumes has under-gone a kind of loss/reduction of function during evolu-tion under domestication. In this study, we alsoobserved that CER1 (eceriferum) genes, which are associ-ated with plant cuticular wax production [63], were sig-nificantly down-regulated in domesticated wheat. Thesefindings are in agreement with higher wax content inthe surface of glumes in wild tetraploid wheat genotypes[64].In addition to the genes that are typically involved incell wall composition, we identified a COBRA gene thatwas expressed only in the glumes of wild emmer wheat(i.e., was down-regulated in domesticated wheat). COBRAencodes for a plant-specific glycosylphosphatidylinositol(GPI)-anchored protein with ω-attachment site at the Cterminus, a hydrophilic central region, a CCVS domain, apotential N-glycosylation site, N-terminal secretion signalsequence, and a predicted cellulose binding site. Extensivestudies have demonstrated that COBRA is critical for bio-synthesis of cell wall constituents comprising structuraltissues of roots, stalks, leaves and other vegetative organs[65]. Likewise, it was suggested recently that genes fromthe COBRA family were involved in deposition ofcrystalline cellulose into different secondary cell wallstructures [66].Among the 22 down-regulated DEGs in the domesti-cated accessions we identified one transcription factorcontaining a NAC domain protein gene. NAC (NAM,ATAF1/2 and CUC2) domain proteins are plant-specifictranscription factors known to play diverse roles in vari-ous plant developmental processes. The NAC domaingene, which was cloned from wild emmer wheat, accel-erates senescence and could enhance nutrient remobili-zation to the developing kernels, thereby improvingZou et al. BMC Genomics  (2015) 16:777 Page 11 of 14their nutritional content [67]. It is noteworthy that inbarley, regulation of gene expression in glumes develop-ment may have direct connection with remobilizationand accumulation of nitrogen in seeds, as was recentlyshown with respect to HvAAP genes [37, 55]. It wasdemonstrated that the shattering genes with a NAC do-main, which functionally activates secondary wall bio-synthesis and promotes the significant thickening ofsecondary walls by its high expression level, are presentin Arabidopsis, rice and soybean genomes [68]. This sug-gests that NAC domain protein may be related to thecontrol of the wheat shattering glumes and may haveplayed a role in cereals and legumes domestication. Ac-cording to our findings on DEGs down-regulated in theglumes of domesticated accessions compared to the wildprogenitor, we can speculate that higher expression ofcell wall controlling genes in wild wheat plays an im-portant role in its glumes toughness.Among the 17 DEGs that were up-regulated in glumesof domesticated wheat compared to the wild progenitor,we identified genes related to fatty acid elongation, fla-vonoid biosynthesis and amino acid transport. The mostabundant up-regulated DEGs in domesticated wheatwere KCS gene family. The KCS gene, a fatty acid elon-gase, determines fatty acid chain length and substratespecificity of the condensation reaction, a rate limitingstep, and its subsequent elongated products like alkanes,aldehydes, primary alcohols, secondary alcohols, ketonesand wax esters [69]. Another example of up-regulatedgenes in domesticated wheat was five CHS genes in-volved in the initial step of flavonoid biosynthesis, in thephenylpropaoid pathway, in pigments production, andplant resistance to biotic and abiotic stresses [70]. Inaddition, we found higher expression of a silicon trans-porter gene in the domesticated wheat which may be re-lated to Si element uptake and distribution [71].As mentioned above, regulation of AAP genes’ expres-sion in barley glumes may play a role in nitrogen remo-bilization and accumulation in seeds [37, 55]. Based onthe over-expression of AAP genes in glumes and in-creased SpHI in domesticated wheat compared to wildprogenitor, we could speculate that dry matter allocationfrom the glumes to grain filling has increased duringwheat evolution under domestication.ConclusionsIn the current study we employed a comparative tran-scriptome profiling of wheat glumes in wild emmer,hulled landraces and modern cultivars. We have identi-fied a few genes showing differential elevated expressionlevels at heading time that may be related to glumestoughness and could probably be involved in wheat evo-lution under domestication. Interestingly, we did notfind any significant differentially expressed genes withAP2 domain similar to Q genes. It is considered that thewheat Q gene confers soft glumes and influences a seriesof traits involved in the control of domestication relatedtraits such as brittle rachis, spike architecture and flow-ering time [14]. Likewise, we did not find differential ex-pression in the Tg that confers glumes toughness. Thisfact may be considered as indirect evidence that thesegenes, start to elevate their expression level after headingtime and culminate before ripening.The advance in new genomic approaches providesnew insight into domestication and evolution under do-mestication. It can facilitate the understanding of theorigin of agriculture, mobilization of the adaptive poten-tial of the wild and landrace germplasms, and finally, forthe rethinking on breeding strategies for the acceleratedimprovement under domestication. Our results showthat in addition to the classical domestication genes,there are many other genes differentially expressed be-tween the wild genotypes, landraces and modern culti-vars, which may be involved in control of agriculturallyimportant traits and basic biological processes, plant de-velopment, cell wall composition, stress tolerance, andpigmentation. The major advantages of RNA-seq tech-nology is that it can assist in unravelling candidategenomic/genetic targets of domestication and improve-ment selection even if nothing is known about the causalselected phenotype and it is not only limited to measur-able phenotypic traits.Availability of supporting dataRaw reads of transcriptome have been deposited into theNCBI Short Read Archive (SRA, http://www.ncbi.nlm.nih.gov/sra/) under the accession numbers: SRR2084071,SRR2084163, SRR2084091, SRR2084165, SRR2084092,and SRR2084160.Additional filesAdditional file 1: Figure S1. Pearson correlation between mappingratios and genetic distance. (PDF 7 kb)Additional file 2: Table S1. DEGs in cultivar compared to wild wheat.(XLSX 246 kb)Additional file 3: Table S2. DEGs in landrace compared to wild wheat.(XLSX 184 kb)Additional file 4: Table S3. DEGs in domesticated compared to wildwheat. (XLSX 90 kb)Additional file 5: Figure S2. Gene Ontology classifications of noveltranscripts. (PDF 97 kb)Additional file 6: Table S4. Differentially expressed transcripts.(XLSX 53 kb)Additional file 7: Table S5. Amplification efficiency of QRT-PCR primers.(DOCX 15 kb)Competing interestsThe authors declare that they have no competing interests.Zou et al. BMC Genomics  (2015) 16:777 Page 12 of 14Authors’ contributionsHZ conducted the data analysis with support from SH. RT and TKcontributed to the field experiment, RNA extraction and RT-PCR analysis. SAand YS supervised the field experiment. HZ and RT wrote the first draft ofthe manuscript. HS, TK, TF, SA, YS and AK participated in preparation of themanuscript. AK conceived and designed the study. All authors have read andapproved the final manuscript.AcknowledgementsThis study was supported by Israel Science Foundation grant # 800–2010. HZis grateful to the Israeli Council for Higher Education and University of Haifafor the postdoctoral fellowship and R.T. is thankful to the Matanel and WolfFoundation for awarding a PhD fellowship. We acknowledge with thank thehelp of Noa Sher from University of Haifa Bioinformatics Service Unit forpreparation of barcoded cDNA libraries.Author details1Department of Evolutionary and Environmental Biology, The Institute ofEvolution, Faculty of Natural Sciences, University of Haifa, Haifa 31905, Israel.2Department of Botany, University of British Columbia, Vancouver, BC V6T1Z4, Canada. 3Robert H. Smith Institute of Plant Sciences and Genetics inAgriculture, The Hebrew University of Jerusalem, Rehovot 7610001, Israel.Received: 2 July 2015 Accepted: 3 October 2015References1. Abbo S, Pinhasi van-Oss R, Gopher A, Saranga Y, Ofner I, Peleg Z. Plantdomestication versus crop evolution: a conceptual framework for cerealsand grain legumes. Trends Plant Sci. 2014;19(6):351–60.2. Dubcovsky J, Dvorak J. Genome plasticity a key factor in the success ofpolyploidy wheat under domestication. Science. 2007;316:1862–6.3. Nevo E, Korol AB, Beiles A, Fahima T. Evolution of wild emmer and wheatimprovement. Population genetics, genetic resources, and genome organizationof wheats progenitor, Triticum dicoccoides. Berlin: Springer; 2002. p. 364.4. Dvorak J, Akhunov ED, Akhunov AR, Deal KR, Luo MC. Molecularcharacterization of a diagnostic DNA marker for domesticated tetraploidwheat provides evidence for gene flow from wild tetraploid wheat tohexaploid wheat. Mol Biol Evol. 2006;23:1386–96.5. Lenser T, Theißen G. Molecular mechanisms involved in convergent cropdomestication. Trends Plant Sci. 2013;18:704–14.6. Cavanagh CR, Chao S, Wang S, Huang BE, Stephen S, Kiani S, et al. Genome-wide comparative diversity uncovers multiple targets of selection forimprovement in hexaploid wheat landraces and cultivars. Proc Natl Acad SciU S A. 2013;110(20):8057–62.7. Peng JH, Ronin YI, Fahima T, Röder MS, Li YC, Nevo E, et al. Domesticationquantitative trait loci in Triticum dicoccoides, the progenitor of wheat. ProcNatl Acad Sci U S A. 2003;100:2489–94.8. Nalam VJ, Vales MI, Watson CJW, Kianian SF, Riera-Lizarazu O. Map-basedanalysis of genes affecting the brittle rachis character in tetraploid wheat(Triticum turgidum L.). Theor Appl Genet. 2006;112(2):373–81.9. Sood S, Kuraparthy V, Bai G, Gill BS. The major threshability genes softglume (sog) and tenacious glume (Tg), of diploid and polyploid wheat,trace their origin to independent mutations at non-orthologous loci. TheorAppl Genet. 2009;119:341–51.10. Peleg Z, Fahima T, Korol AB, Abbo S, Saranga Y. Genetic analysis of wheatdomestication and evolution under domestication. J Exp Bot. 2011;62:5051–61.11. Tzarfati R, Barak V, Fahima T, Abbo S, Saranga Y, Korol AB. Novel quantitativetrait loci underlying major domestication traits in tetraploid wheat. MolBreeding. 2014;34:1613–28.12. Huang X, Han B. Natural variations and genome-wide association studies incrop plants. Annu Rev Plant Biol. 2014;65:531–51.13. Simons KJ, Fellers JP, Trick HN, Zhang Z, Tai YS, Gill BS, et al. Molecularcharacterization of the major wheat domestication gene Q. Genetics.2006;172(1):547–55.14. Zhang Z, Belcram H, Gornicki P, Charles M, Just J, Huneau C, et al.Duplication and partitioning in evolution and function of homoeologous Qloci governing domestication characters in polyploid wheat. Proc Natl AcadSci U S A. 2011;108(46):18737–42.15. Onishi I, Hongo A, Sasakuma T, Kawahara T, Kato K, Miura H. Variation andsegregation for rachis fragility in spelt wheat, Triticum spelta L. GenetResour Crop Evol. 2006;53:985–92.16. Tzarfati R, Saranga Y, Barak V, Gopher A, Korol AB, Abbo S. Threshing efficiency asan incentive for rapid domestication of emmer wheat. Ann Bot. 2013;112:829–37.17. Harlan JR, De Wet JMJ, Price EG. Comparative evolution of cereals.Evolution. 1973;27:311–25.18. Davey JW, Hohenlohe PA, Etter PD, Boone JQ, Catchen JM, Blaxter ML.Genome-wide genetic marker discovery and genotyping using next-generation sequencing. Nat Rev Genet. 2011;12(7):499–510.19. Harper AL, Trick M, Higgins J, Fraser F, Clissold L, Wells R, et al. Associativetranscriptomics of traits in the polyploid crop species Brassica napus. NatBiotechnol. 2012;30(8):798–802.20. Huang X, Zhao Y, Wei X, Li C, Wang A, Zhao Q, et al. Genome-wideassociation study of flowering time and grain yield traits in a worldwidecollection of rice germplasm. Nat Genet. 2012;44:32–9.21. Hufford MB, Xu X, van Heerwaarden J, Pyhäjärvi T, Chia JM, Cartwright RA,et al. Comparative population genomics of maize domestication andimprovement. Nat Genet. 2012;44(7):808–11.22. Li YH, Zhao SC, Ma JX, Li D, Yan L, Li J, et al. Molecular footprints ofdomestication and improvement in soybean revealed by whole genomere-sequencing. BMC Genomics. 2013;14:579.23. Wang S, Wong D, Forrest K, Allen A, Chao S, Huang BE, et al. Characterizationof polyploid wheat genomic diversity using a high-density 90 000 singlenucleotide polymorphism array. Plant Biotechnol J. 2014;12(6):787–96.24. Ergen NZ, Budak H. Sequencing over 13000 expressed sequence tags fromsix subtractive cDNA libraries of wild and modern wheats following slowdrought stress. Plant Cell Environ. 2009;32(3):220–36.25. Ayal S, Ophir R, Levy AA. Genomics of tetraploid wheat domestication. In:Tsunewaki K, editor. Frontiers of Wheat Bioscience, the 100th MemorialIssue of Wheat Information Service. Yokohama: Kihara Memorial Foundationfor the Advancement of Life Sciences; 2005. p. 185–203.26. Stamova BS, Laudencia-Chingcuanco D, Beckles DM. Transcriptomic analysisof starch biosynthesis in the developing grain of hexaploid wheat. Int JPlant Sci. 2009;2009:407426. doi:10.1155/2009/407426.27. Singh A, Mantri S, Sharma M, Chaudhury A, Tuli R, Roy J. Genome-widetranscriptome study in wheat identified candidate genes related to processingquality, majority of them showing interaction (quality x development) andhaving temporal and spatial distributions. BMC Genomics. 2014;15:29.28. Ergen NZ, Thimmapuram J, Bohnert HJ, Budak H. Transcriptome pathwaysunique to dehydration tolerant relatives of modern wheat. Funct IntegrGenomics. 2009;9(3):377–96.29. Krugman T, Chagué V, Peleg Z, Balzergue S, Just J, Korol AB, et al.Multilevel regulation and signalling processes associated withadaptation to terminal drought in wild emmer wheat. Funct IntegrGenomics. 2010;10:167–86.30. Krugman T, Peleg Z, Quansah L, Chagué V, Korol AB, Nevo E, et al.Alteration in expression of hormone-related genes in wild emmer wheatroots associated with drought adaptation mechanisms. Funct IntegrGenomics. 2011;11:565–83.31. Swanson-Wagner R, Briskine R, Schaefer R, Hufford MB, Ross-Ibarra J, MyersCL, et al. Reshaping of the maize transcriptome by domestication. Proc NatlAcad Sci U S A. 2012;109:11878–83.32. Yoo MJ, Wendel JF. Comparative evolutionary and developmental dynamicsof the cotton (Gossypium hirsutum) fiber transcriptome. PLoS Genet.2014;10(1), e1004073.33. Bellucci E, Bitocchi E, Ferrarini A, Benazzo A, Biagetti E, Klie S, et al.Decreased nucleotide and expression diversity and modified coexpressionpatterns characterize domestication in the common Bean. Plant Cell.2014;26(5):1901–12.34. Wang Z, Gerstein M, Snyder M. RNA-Seq. a revolutionary tool fortranscriptomics. Nat Rev Genet. 2009;10(1):57–63.35. Faris JD, Zhang Z, Chao S. Map-based analysis of the tenacious glume geneTg-B1 of wild emmer and its role in wheat domestication. Gene.2014;542(2):198–208.36. Wang ZM, Wei AL, Zheng DM. Photosynthetic characteristics of non-leaforgans of winter wheat cultivars differing in ear type and their relationshipwith grain mass per ear. Photosynthetica. 2001;39(2):239–44.37. Kohl S, Hollmann J, Erban A, Kopka J, Riewe D, Weschke W, et al. Metabolic andtranscriptional transitions in barley glumes reveal a role as transitory resourcebuffers during endosperm filling. J Exp Bot. 2015. doi:10.1093/jxb/eru492.Zou et al. BMC Genomics  (2015) 16:777 Page 13 of 1438. Simpson RJ, Lambers H, Dalling MJ. Nitrogen redistribution during grain growthin wheat (Triticum aestivum L.) IV. Development of a quantitative model of thetranslocation of nitrogen to the grain. Plant Physiol. 1983;71(1):7–14.39. Okamoto Y, Takumi S. Pleiotropic effects of the elongated glume gene P1on grain and spikelet shape-related traits in tetraploid wheat. Euphytica.2013;194:207–18.40. Wang H, Nussbaum-Wagler T, Li B, Zhao Q, Vigouroux Y, Faller M, et al. Theorigin of the naked grains of maize. Nature. 2005;436(7051):714–9.41. Brenchley R, Spannagl M, Pfeifer M, Barker GL, D'Amore R, Allen AM, et al.Analysis of the bread wheat genome using whole-genome shotgunsequencing. Nature. 2012;491(7426):705–10.42. Bolger AM, Lohse M, Usadel B. Trimmomatic: aflexible trimmer for Illuminasequence data. Bioinformatics. 2014. doi:10.1093/bioinformatics/btu170.43. Liao Y, Smyth GK, Shi W. The Subread aligner: fast, accurate and scalableread mapping by seed-and-vote. Nucleic Acids Res. 2013;41(10), e108.44. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general-purpose program forassigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.45. Anders S, Huber W. Differential expression analysis for sequence count data.Genome Biol. 2010;11(10):R106.46. Conesa A, Götz S, García-Gómez JM, Perol J, Talón M, Robles M. Blast2GO: auniversal tool for annotation, visualization and analysis in functionalgenomics research. Bioinformatics. 2005;21:3674–6.47. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, et al.Gapped BLAST and PSI-BLAST: a new generation of protein database searchprograms. Nucleic Acids Res. 1997;25:3389–402.48. Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z, et al. WEGO: a web toolfor plotting GO annotations. Nucleic Acids Res. 2006;34:293–7.49. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al.De novo transcript sequence reconstruction from RNA-seq using the Trinityplatform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512.50. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficientalignment of short DNA sequences to the human genome. Genome Biol.2009;10:R25.51. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq datawith or without a reference genome. BMC Bioinformatics. 2011;12:323.52. Robinson MD, McCarthy DJ, Smyth GK. EdgeR: a Bioconductor package fordifferential expression analysis of digital gene expression data.Bioinformatics. 2010;26(1):139–40.53. Livak KJ, Schmittgen TD. Analysis of relative gene expression data usingreal-time quantitative PCR and the 2(−Delta Delta C(T)) Method. Methods.2001;25(4):402–8.54. Gepts P. The contribution of genetic and genomic approaches to plantdomestication studies. Curr Opin Plant Biol. 2014;18:51–9.55. Kohl S, Hollmann J, Blattner FR, Radchuk V, Andersch F, Steuernagel B, et al.A putative role for amino acid permeases in sink-source communication ofbarley tissues uncovered by RNA-seq. BMC Plant Biol. 2012;12:154.56. Akpinar BA, Yuce M, Lucas S, Vrána J, Burešová V, Doležel J, et al. Molecularorganization and comparative analysis of chromosome 5B of the wildwheat ancestor Triticum dicoccoides. Sci Rep. 2015;5:10763.57. Krasileva KV, Buffalo V, Bailey P, Pearce S, Ayling S, Tabbita F, et al.Separating homeologs by phasing in the tetraploid wheat transcriptome.Genome Biol. 2013;14(6):R66.58. Kaur S, Dhugga K, Gill K, Singh J. Functional Informatics of cellulose synthasegenes in wheat. Plant & Animal Genome XXIII, San Diego, CA; 2015; P0015.59. Faik A, Abouzouhair J, Sarhan F. Putative fasciclin-like arabinogalactan-proteins (FLA) in wheat (Triticum aestivum) and rice (Oryza sativa):identification and bioinformatic analyses. Mol Genet Genomics.2006;276(5):478–94.60. MacMillan CP, Mansfield SD, Stachurski ZH, Evans R, Southerton SG.Fasciclin-like arabinogalactan proteins: specialization for stem biomechanicsand cell wall architecture in Arabidopsis and Eucalyptus. Plant J.2010;62(4):689–703.61. Bischoff V, Nita S, Neumetzler L, Schindelasch D, Urbain A, Eshed R, et al.TRICHOME BIREFRINGENCE and its homolog AT5G01360 encode plant-specific DUF231 proteins required for cellulose biosynthesis in Arabidopsis.Plant Physiol. 2010;153(2):590–602.62. Zhao Q, Nakashima J, Chen F, Yin Y, Fu C, Yun J, et al. Laccase is necessaryand nonredundant with peroxidase for lignin polymerization duringvascular development in Arabidopsis. Plant Cell. 2013;25(10):3976–87.63. Hu X, Zhang Z, Li W, Fu Z, Zhang S, Xu P. cDNA cloning and expressionanalysis of a putative decarbonylase TaCer1 from wheat (Triticum aestivumL.). Acta Physiol Plant. 2009;31:1111–8.64. Wang J, Li W, Wang W. Fine mapping and metabolic and physiologicalcharacterization of the glume glaucousness inhibitor locus Iw3 derived fromwild wheat. Theor Appl Genet. 2014;127(4):831–41.65. Cao Y, Tang X, Giovannoni J, Xiao F, Liu Y. Functional characterization of atomato COBRA-like gene functioning in fruit development and ripening.BMC Plant Biol. 2012;12:211.66. Ben-Tov D, Abraham Y, Stav S, Thompson K, Loraine A, Elbaum R, et al.COBRA-LIKE2, a Member of the Glycosylphosphatidylinositol-AnchoredCOBRA-LIKE Family, Plays a Role in Cellulose Deposition in Arabidopsis SeedCoat Mucilage Secretory Cells. Plant Physiol. 2015;167(3):711–24.67. Uauy C, Distelfeld A, Fahima T, Blechl A, Dubcovsky J. A NAC Generegulating senescence improves grain protein, zinc, and iron content inwheat. Science. 2006;314(5803):1298–301.68. Dong Y, Yang X, Liu J, Wang BH, Liu BL, Wang YZ. Pod shattering resistanceassociated with domestication is mediated by a NAC gene in soybean. NatCommun. 2014;5:3352.69. Lokesh U, Kiranmai K, Pandurangaiah M, Sudhakarbabu O, Nareshkumar A,Sudhakar C. Role of plant fatty acid elongase (3 keto acyl-CoA synthase) genein cuticular wax biosynthesis. Res Rev: J Agric Allied Sci. 2013;2(4):35–42.70. Trojann V, Musilováa M, Vyhnáneka T, Klejdusb B, Hanáčeka P, Havela L.Chalcone synthase expression and pigments deposition in wheat withpurple and blue colored caryopsis. J Cereal Sci. 2014;1:48–55.71. Yamaji N, Chiba Y, Mitani-Ueno N, Feng Ma J. Functional characterization ofa silicon transporter gene implicated in silicon distribution in barley. PlantPhysiol. 2012;160(3):1491–7.Submit your next manuscript to BioMed Centraland take full advantage of: • Convenient online submission• Thorough peer review• No space constraints or color figure charges• Immediate publication on acceptance• Inclusion in PubMed, CAS, Scopus and Google Scholar• Research which is freely available for redistributionSubmit your manuscript at www.biomedcentral.com/submitZou et al. BMC Genomics  (2015) 16:777 Page 14 of 14


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items