Genome-wide characterization and expression profiling of HD-Zip gene family related to abiotic ...

20 Pages·2017·2.66 MB·English

RESEARCHARTICLE
Genome-wide characterization and
expression profiling of HD-Zip gene family
related to abiotic stress in cassava
ZehongDing*☯,LiliFu☯,YanYan,WeiweiTie,ZhiqiangXia,WenquanWang,MingPeng,
WeiHu*,JiamingZhang*
KeyLaboratoryofBiologyandGeneticResourcesofTropicalCrops,InstituteofTropicalBioscienceand
Biotechnology,ChineseAcademyofTropicalAgriculturalSciences,Haikou,Hainan,China
☯Theseauthorscontributedequallytothiswork.
a1111111111
*[email protected](ZD);[email protected](WH);[email protected](JZ)
a1111111111
a1111111111
a1111111111
Abstract
a1111111111
Homeodomain-leucinezipper(HD-Zip)genefamilyplaysimportantrolesinvariousabiotic
stressesandhormonesignalinginplants.However,noinformationiscurrentlyavailable
regardingthisfamilyincassava(Manihotesculenta),animportantdrought-tolerantcropin
OPENACCESS tropicalandsub-tropicalareas.Here,57HD-Zipgenes(MeHDZ01-57)wereidentifiedinthe
Citation:DingZ,FuL,YanY,TieW,XiaZ,Wang cassavagenome,andtheywereclassifiedintofoursubfamiliesbasedonphylogeneticanal-
W,etal.(2017)Genome-widecharacterizationand ysis,whichwasfurthersupportedbytheirgenestructureandconservedmotifcharacteris-
expressionprofilingofHD-Zipgenefamilyrelated
tics.Ofwhichfivegenepairswereinvolvedinsegmentalduplicationbutnonefortandem
toabioticstressincassava.PLoSONE12(3):
e0173043.doi:10.1371/journal.pone.0173043 duplication,suggestingthatsegmentalduplicationwasthemaincausefortheexpansionof
MeHDZgenefamilyincassava.GlobalexpressionprofilesrevealedthatMeHDZgenes
Editor:Jin-SongZhang,InstituteofGeneticsand
DevelopmentalBiologyChineseAcademyof wereconstitutivelyexpressed,ornotexpressed,ortissue-specificexpressedinexamined
Sciences,CHINA tissuesinbothcultivatedandwildsubspecies.Transcriptomicanalysisofthreegenotypes
Received:September2,2016 showedthatmostofMeHDZgenesrespondeddifferentlytodroughtandpolyethylene
glycoltreatments.Subsequently,quantitativeRT-PCRanalysisrevealedcomprehensive
Accepted:February14,2017
responsesoftwelveselectedMeHDZgenestovariousstimuliincludingcold,salt,andABA
Published:March1,2017
treatments.ThesefindingswillincreaseourunderstandingofHD-Zipgenefamilyinvolved
Copyright:©2017Dingetal.Thisisanopen inabioticstressesandsignalingtransduction,andwillprovideasolidbaseforfurtherfunc-
accessarticledistributedunderthetermsofthe
tionalcharacterizationofMeHDZgenesincassava.
CreativeCommonsAttributionLicense,which
permitsunrestricteduse,distribution,and
reproductioninanymedium,providedtheoriginal
authorandsourcearecredited.
DataAvailabilityStatement:Allrelevantdataare
Introduction
withinthepaperanditsSupportingInformation
files.
Homeodomain-leucinezipper(HD-Zip)proteinsareoneofkeytranscriptionfactorsinplants.
Funding:Thisstudywassupportedbythe Basedonsequenceanalysis,HD-Zipproteinswerecharacterizedbythepresenceoftwofunc-
InternationalScienceandTechnologyCooperation tionaldomains:ahomeodomain(HD)responsibleforspecificbindingtoDNAandaleucine
ProgramofChina(2010DFA62040,
zipper(LZ)motifwhichcloselylinkedtoHDandactedasadimerizationmotif[1,2].To
2014DFA30680),theNaturalScienceFoundationof
date,manymembersofHD-Zipproteinshavebeenfoundinalargenumberofplantspecies,
China(31471561),andtheNaturalScience
includingArabidopsis[3],rice[4],maize[5],soybean[6],legume[7],banana[8],andsoon.
FoundationofHainanProvince(20163120,
20153048). Basedonthestructures,HD-Zipproteinsweremainlycategorizedintofourgroups:HD-ZipI
PLOSONE|DOI:10.1371/journal.pone.0173043 March1,2017 1/20
HD-Zipgenefamilyincassava
Competinginterests:Theauthorshavedeclared toHD-ZipIV[2].MembersofHD-ZipIandIIcontainedtheconservedHDandLZdomains
thatnocompetinginterestsexist. andboundtothesimilarsequenceCAAT-N-ATTG[9].Inaddition,HD-ZipIImemberscon-
tainedothermotifs,suchasCPSCE[2].BesidestheHDandLZdomains,membersofboth
HD-ZipIIIandIValsocontainedaSTARTdomainwithputativelipidbindingcapability
[10].However,aconservedC-terminalMEKHLAmotif,whichwasprobablyinvolvedinoxy-
genredoxandlightsignaling,waspresentinthemembersofHD-ZipIIIbutabsentinthe
HD-ZipIVsubfamily[11].
HD-Zipproteinshavebeendemonstratedtoparticipateinmanybiologicalprocesses
relatedtoplantgrowthanddevelopment[1,2].MembersofHD-ZipIsubfamilymainlypartic-
ipatedinabioticstressessuchasdrought,extremetemperatures,osmoticandlightstresses
[2,3].Forexample,AtHB5,-6,-7and-12,whichbelongedtoArabidopsisHD-ZipIsubfamily,
wereup-regulatedordown-regulatedbydroughtstressorexternallyappliedabscisicacid
(ABA)[12–15].Othermembersofthissubfamily,suchasAtHB52andAtHB53,wereup-regu-
latedindarkcondition[3].HD-ZipIIproteinsweremainlyinvolvedinshadeavoidanceand
auxinresponse[16].Forexamples,AtHB2inthissubfamilywasspecificallyandreversiblyreg-
ulatedbyredtofar-redratio,anditwasparticipatedintheregulationofArabidopsisshade
avoidanceresponse[17].HAT2, anothermemberofHD-ZipIIsubfamily,wasisolatedasan
auxin-induciblegenethroughDNAmicroarrayscreening[18].HD-ZipIIImembersmainly
participatedinthedevelopmentalregulationofapicalmeristem,vascularbundles,auxintrans-
portandlateralorganinitiation[2].Forexamples,AtHB8andAtHB15wereinvolvedinvascu-
lardevelopment[19,20].SeveralstudieshavedemonstratedthatREV,PHBandPHV,along
withKANADI,wereinvolvedincontrollingtheabaxial-adaxialpatterningoflateralorgans
[21].HD-ZipIVproteinsgenerallyplayedimportantrolesinanthocyaninaccumulation,epi-
dermalcelldifferentiation,trichomeformationandrootdevelopment[1,2].Forexamples,
atml1pdf2doublemutantsresultedinseveredefectsinshootepidermalcelldifferentiation
andfailedtosurviveaftergermination[22].ANL2affectedanthocyaninaccumulationinsub-
epidermaltissues[23],whileothermembersinthissubfamily,suchasHDG11,HDG12and
GL2,wereinvolvedintrichomedevelopment[2,24].Recently,manystudieshavedemon-
stratedthatHD-ZipIIIandIVsubfamilieswerealsoinvolvedinabioticstressessuchassalt,
droughtandABAtreatment[6,8,25].
Cassava(Manihotesculenta)isoneofthemostimportantrootcropsintropicalandsub-
tropicalareasandprovidesstaplefoodforover700millionpeoplearoundtheworld[26].Due
toitsgoodperformanceunderadverseconditionssuchasdroughtandlowfertilization,cas-
savahasbeenconsideredasanimportantfoodsecuritycropworldwide.However,asarepre-
sentativetropicalcrop,cassavaisverysensitivetocoldstress[27].Todate,themechanisms
underlyingitsabioticstressesremainelusive.Withtheavailabilityofthecassavagenome
sequences[28,29],itprovidesanexcellentopportunityforgenome-wideannotation,classifica-
tionandcomparativegenomeresearch.AlthoughHD-Zipgeneshavebeenextensivelycharac-
terizedinmanyspecies,asystemicstudyofHD-Zipgenefamily(especiallyforthestress-
responsivemembers)hasnotbeenperformedincassava.Inthisstudy,atotalof57HD-Zip
geneswereidentifiedandcharacterizedinthecassavagenome.Besidestheirbasiccharacteris-
ticsincludingphylogeneticrelationship,genestructure,conservedproteinmotifandgenome
location,theexpressionprofilesoftheseHD-Zipgeneswereinvestigatedindifferentorgansof
variousgenotypesandseveralstressconditions,includingdrought,polyethyleneglycol(PEG),
cold,salt,andABAtreatmentsthroughRNA-seqorquantitativeRT-PCR(qRT-PCR).These
findingswillprovideausefulfoundationforfurtherinvestigationsintothemolecularbiologi-
calfunctionsoftheHD-Zipgenesincassava.
PLOSONE|DOI:10.1371/journal.pone.0173043 March1,2017 2/20
HD-Zipgenefamilyincassava
Materialsandmethods
Plantmaterialsandstresstreatments
TorevealtheexpressionofHD-Zipgeneseitherindifferenttissuesorinresponsetodrought
stressincassava,severalpreviouslygeneratedRNA-seqdata[30,31]wereused.Intotalfour
cassavaaccessions,includingonewildsubspecies(W14,M.esculentasubsp.flabellifolia)and
threecultivatedvarieties(Arg7,Ku50andSC124),wereusedinthisstudy.Ofwhich,Arg7is
adaptedtohigh-latitudegeographicalregionofArgentinaandofmoderatedrought-tolerance;
SC124isawidelyplantedcultivarinChinaandcansurviveinseveredroughtcondition,while
Ku50isarepresentativecultivarwithhighrootyieldandhighstarchcontentinroottubers
andcangrowunderunfavorableenvironmentalconditions[29,32,33].
Incassavagrowingseason,Arg7,Ku50andW14wereplantedinthefarmofChineseAcad-
emyofTropicalAgriculturalSciences(Haikou,China)inthesummerof2013.Todetectthe
changesofgeneexpressionindifferentcassavatissues,samplesfromleaf(90daysafterplant-
ing),stem(90daysafterplanting),earlystorageroot(ESR,90daysafterplanting),middlestor-
ageroot(MSR,150daysafterplanting)andlaststorageroot(LSR,210daysafterplanting)
werecollectedforRNA-seqanalysis.
Toinvestigatethechangesofgeneexpressioninresponsetodroughtstress,stemsof
about15cminlengthandtwotothreebudswereselectedandplantedverticallyinpots
(sand:vermiculite=1:1;height×upperdiameter×bottomdiameter=18.8×18.5×14.8
cm),aspreviouslydescribed[34].About45dayslater,twoexperimentswereconducted:(1)
seedlingsofArg7,SC124andW14werewithheldwaterfor0day(control)and12days,and
theirleavesandrootswerecollectedforRNA-seq,respectively;(2)seedlingsofKu50were
wateredby20%polyethyleneglycol(PEG)6000solution,thentheleaves,includingfolded
leaf(FL),fullexpandedleaf(FEL)andbottomleaf(BL),aswellasroot(RT)werecollected
at0,3,24hafterthetreatmentforRNA-seq,respectively.Eachsamplewaspooledfrom
15plants.
Toexaminetheresponseofcassavatodifferentabioticstresses,includingcold,salt,and
ABAtreatment,2-month-oldcassavaseedlingsofArg7wereusedandtheirleaveswerecol-
lectedforqRT-PCR.Forsaltstress,cassavaseedlingswerewateredwith300mMNaClsolu-
tionfor24days.ForABAtreatment,cassavaseedlingsweresprayedwith100μMABAfor72
h.Forcoldstress,cassavaseedlingswereexposedto4˚Cfor48h,andthenreturnedtonormal
growthconditionsfor14daysofrecovery.
IdentificationofHD-Zipfamilygenesincassava
Thewholegenomeandproteinsequencesofcassavaweredownloadedfromphytozome
databasev9.0(https://phytozome.jgi.doe.gov/).HiddenMarkovModel(HMM)profilesof
homeobox(HD,PF00046)andhomeoboxassociatedleucinezipper(HALZ,PF02183)were
downloadedfromPFam(http://pfam.sanger.ac.uk/)andwereusedasqueriestoidentify
putativecassavaHD-ZipgenesbylocalHMM-basedsearchessettingE-values<0.01using
HMMER3[35].Additionally,proteinsequencesofHD-ZipgenesfromArabidopsisandrice
[36]wereusedinaBLASTsearchusingBLASTPagainstcassavaproteinstoexploremore
HD-ZipswhichmightbemissedbyHMMprofilesearching.Afterremovingredundant
sequences,eachHD-Zipcandidatewasfurtherconfirmedbyinvestigatingconserveddomains
usingdatabaseofCDD(http://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi)inauto-
maticmode(threshold=0.01,maximumhits=500)andPFamsettingE-valueequalto1.0,
respectively.ThemolecularweightandisoelectricpointsofeachHD-Zipproteinwerecalcu-
latedusingExPASytools(http://www.expasy.org/tools/).
PLOSONE|DOI:10.1371/journal.pone.0173043 March1,2017 3/20
HD-Zipgenefamilyincassava
Phylogeneticandsequenceanalysis
BesidesHD-ZipsfromArabidopsisandrice,proteinsequencesofpublishedHD-Zipgenes
fromVitisvinifera,Medicagotruncatula,Populustrichocarpa,andmaize[5,36]werealsoused.
ProteinsequencesofHD-ZipgeneswerealignedusingMUSCLE[37].Phylogenetictreewas
conductedinMEGA6[38]usingMaximum-LikelihoodmethodwithJones-Taylor-Thornton
(JTT)aminoacidsubstitutionmodel.Bootstrapanalysiswasperformedusing1000replicates
withthepartialdeletionmodelforgaps/missingdata.
HD-Zipexon-intronstructurewasperformedwithGeneStructureDisplayServer(GSDS,
http://gsds.cbi.pku.edu.cn/)[39]byaligningcDNAtotheircorrespondinggenomicDNA
sequences.ProteinmotifswerepredictedbyMEMEsoftware[40]withoptimummotifwidth
rangedfrom11to50andsettingthemaximumnumberofmotifsequalto15.Thepredicted
motifsofHD-ZipproteinswereannotatedbysearchingagainstInterProScandatabase(http://
www.ebi.ac.uk/Tools/pfa/iprscan/).
TodemonstratethepossibleregulatorymechanismofHD-Zipgenesinvariousstress
responses,1500bppromoterregionupstreamofthetranscriptionstartsite(ATG)wasusedto
analyzecis-elementsusingPlantCAREdatabase(http://bioinformatics.psb.ugent.be/webtools/
plantcare/html/).Theoverrepresentationanalysisofcis-elementswasperformedbyexact
binominaltestineachsubfamilyofHD-Zipgenes.
Chromosomallocationandgeneduplication
CassavaHD-Zip(referredtoMeHDZ)genesweremappedoncassavachromosomesaccording
totheirpositionsintheearlygenomeversion6.1fromphytozomedatabase(https://
phytozome.jgi.doe.gov/).TheschematicdiagramofMeHDZgeneswasdrawnbyMapInspect
(http://www.plantbreeding.wur.nl/UK/software_mapinspect.html).
Tandemduplicatesweredefinedasgenesthatoccurredwithina50kbregion.Segmental
duplicatesweredetectedbyBLASTN(score<1e-5)using100kb(containing50kbupstream
and50kbdownstream)sequencesflankingtothegenes,andwerefurtherconfirmedbychain-
ingalignments(alignmentlength>200bpandsequenceidentity>85%)[41].
TranscriptomeandqRT-PCRanalysis
TotalRNAwasextractedusingRNAplantreagentkits(TiangenCompany,Beijing,China).
EachRNA-SeqlibrarywasconstructedaspreviouslydescribedinYangetal[42].Then,librar-
ieswereindexed,pooledandsequencedonIlluminaHiSeq2000.Transcriptomeanalysiswas
performedaspreviouslydescribed[33,34].Forexamples,sequencequalitywasexamined
usingFastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/),cleanreadswere
mappedtocassavagenome(version4.1)usingTophatv2.0.13[43],andrawcountdatawere
obtainedandnormalizedbyCuffdiffembeddedinCufflinkspipelinev2.1.1[44].Thetran-
scriptomedataweresubmittedtoNCBIdatabaseandtheiraccessionswereshowninS1Table.
ExpressionofMeHDZgenesinresponsetovariousstimuli,includingcold,salt,andABA
treatments,wasexaminedbyqRT-PCRmethodusingSYBR-green(TaKaRaBiotechnology
Co.Ltd,Dalian,China)andStratageneMx3005Psystem(Stratagene,CA,USA)aspreviously
described[34].Thecassavaactingene[45]wasusedasaninternalcontrol.Thespecificprim-
ersoftargetMeHDZgenesweredesignedbyPrimer5.0software(S2Table).Subsequently,the
specificityofeachprimerpairwasexaminedbyqRT-PCRmeltingcurveanalysis,agarosegel
electrophoresis,andPCRproductssequencing[30].Foreachsample,qRT-PCRreactionwas
performedwiththreeindependentbiologicalreplicatesandtherelativemRNAexpression
levelwascalculatedby2-ΔΔCtasbefore[42].StatisticaldifferencewasexaminedbyDuncan’s
PLOSONE|DOI:10.1371/journal.pone.0173043 March1,2017 4/20
HD-Zipgenefamilyincassava
multiplerangetest(n=3)andexpressionvalueswithdifferentlettersweresignificantat
P<0.05.
Results
IdentificationofHD-Zipgenesincassava
Intotal,125HD-Zipcandidatesequenceswerefoundthroughsearchingagainstcassava
genomeusingHMMandBLASTPmethods.Afterremovingredundantsequencesandcon-
firmingthepresenceofbothHDandLZdomainsbyCDDandPFamdatabases,57non-
redundantHD-Zipgenes(designatedasMeHDZ01-57)werefinallyretainedandusedforfur-
theranalysis.Theseproteinsvariedfrom79(MeHDZ57)to853(MeHDZ01)aminoacidin
length,andtheirmolecularweightrangedfrom9.12(MeHDZ05)to93.83(MeHDZ13)kDa
(S3Table).
PhylogeneticandevolutionaryanalysisofHD-Zipgenes
TostudytheevolutionaryrelationshipofHD-Zipgenesbetweencassavaandotherspecies,
proteinsequencesof57MeHDZgenes,togetherwith48fromArabidopsisand44fromrice
[36],werealignedandusedforphylogeneticanalysis(S4TableandS1File).AsshowninFig1,
phylogenetictreewellclusteredtheseMeHDZproteinsintofourgroups(subfamilyI-IV),
Fig1.PhylogeneticclusteringofHD-ZipproteinsfromArabidopsis,riceandcassava.Theunrooted
phylogenetictreewasconstructedbyMaximum-Likelihoodmethodwith1000bootstrapreplicates.The
HD-Zipproteinsequencescorrespondingto44fromrice(prefixedwith‘Os’),48fromArabidopsis(prefixed
with‘At’),and57fromcassava(prefixedwith‘Me’)werewellseparatedintofourmajorgroups(I-IV).
doi:10.1371/journal.pone.0173043.g001
PLOSONE|DOI:10.1371/journal.pone.0173043 March1,2017 5/20
HD-Zipgenefamilyincassava
togetherwithknownHD-ZipItoIVmembersofriceandArabidopsisrespectivelyineach
clade,indicatingthattherewerefourmajortypesofMeHDZgenesincassava.
Inaddition,HD-Zipproteinsfromsevenspeciesincludingfiveeudicots(Arabidopsis,cas-
sava,Vitisvinifera,Medicagotruncatula,andPopulustrichocarpa)andtwomonocots(riceand
maize)werealsoalignedandfurtherdividedintosubclasses(S1FigandS2File),basedon
whichfourdifferentevolutionaryscenarioswererevealed:(1)HD-Zipgenesthatweregener-
atedbeforethedivergenceofeudicotsandmonocotsandtheywerehighlyconserved,asthey
werepresentedinallofthesevenspeciesalthoughthenumberofHD-Zipgenesmightbevar-
ied.Moreover,thisphenomenonwasobservedineachsubfamily,e.g.,subclassa,f,iinsub-
familyI,subclassm,pinsubfamilyII,subclasstinsubfamilyIII,andsubclassDinsubfamily
IV;(2)HD-Zipgenesthatweregeneratedbeforethedivergenceofeudicotsandmonocotsbut
theirhomologsinsomespeicesmightbelostduringtheevolution,becausetheywerepre-
sentedinbotheudicotsandmonocotsbutnotinallofexaminedspecies,e.g.,subclassj,kin
subfamilyII,subclassq,rinsubfamilyIII,andsubclassu,winsubfamilyIV;(3)HD-Zipgenes
thatweregeneratedafterthedivergenceofeudicotsandmonocots,astheywereexclusively
presentedineudicots.Theywereeitherhighlyconserved(e.g,subclassc,h,n,z,A,andCcon-
tainingHD-Zipgenespresentedinalloffiveexaminedeudicots)ornotveryconserved(e.g,
subclassb,d,s,andycontainingHD-Zipgenesonlypresentedinpartialofexaminedeudicots);
(4)HD-Zipgenesthatweregeneratedafterthedivergenceofeudicotsandmonocots,butthey
wereexclusivelypresentedinmonocots.Forexamples,subclassg,l,o,x,B,andFonlycon-
tainedHD-Zipgenesfromriceandmaizeinthisstudy.Takentogether,theseresultsindicated
acomplicatedevolutionaryprocessofHD-Zipgenesinplants.
Genestructureandconservedmotifanalysis
TogaininformationregardingthegenestructureofMeHDZgenes,exon-intronstructure
analysiswasperformedforeachMeHDZgeneindividually(Fig2).Overall,mostclosely
relatedmemberswithinthesamesubfamilysharedsimilarexon-intronstructureandintron
numbers.Therewere23membersclusteredinMeHDZIsubfamilyandmostofthemcon-
tainedonetothreeexons.Therewere14membersinMeHDZIIsubfamily,andtheycon-
tainedthreetofourexons.Therewere9and11membersinMeHDZIIIandIVsubfamilies,
respectively.ComparedwithMeHDZIandIIsubfamilies,muchmoreexonswerefoundin
MeHDZIII(18ineach)andIV(9.6onaverage)subfamilies,respectively.
Intotal,15conservedmotifswerecapturedbyMEMEsoftwareandsubsequentlyannotated
withInterProScandatabase(Fig3,S2FigandS5Table).Asexpected,HDdomain(motif1
and2)andleucinezipperdomain(motif4)werefoundinallidentifiedMeHDZmembers.
SomemotifsthatpresentedonlyincertainMeHDZsubfamilieswereobserved.Forexamples,
STARTdomain,whichwasrepresentedbymotif3,5and7,wasnotfoundinmembersof
MeHDZIandIIsubfamiliesbutpresentedinallmembersofMeHDZIIIandIVsubfamilies;
MEKHLAdomain,representedbymotif6,waspresentedonlyinmembersofMeHDZIII
subfamily.Besidesthewellannotated,someunknownmotifs(e.g.,motif9)werealsofound
specificallyinsomesubfamilies,indicatingthatthesemotifsmightbekeyelementsforthesub-
family-specificfunctions.
ChromosomallocationandgeneduplicationofMeHDZgenes
GenomicdistributionofMeHDZgeneswasdeterminedbytheirchromosomalpositions.To
thebestofourknowledge,thisisthefirstreportconcerningthegenomicdistributionofgene
familyincassava.ExceptforthreemembersofMeHDZ31,-47and-53thatwerelocatedon
scaffolds,theremaining54MeHDZgenesweremappedon18cassavachromosomes(Fig4;S6
PLOSONE|DOI:10.1371/journal.pone.0173043 March1,2017 6/20
HD-Zipgenefamilyincassava
Fig2.Exon-intronstructureof57MeHDZgenesincassava.Thephylogenetictreewasconstructedby
Maximum-Likelihoodmethodwith1000bootstrapreplicates.TheMeHDZgeneswereclusteredintofour
groups,ItoIV.Exon-intronanalysiswasperformedusingGSDS(http://gsds.cbi.pku.edu.cn/).Lengthofexons
andintronsofeachMeHDZgenewasdisplayedproportionally.Upstreamwasdefinedasthesequencesbefore
thetranscriptionstartsite‘ATG’,whiledownstreamwasdefinedasthesequencesafterthestopcodon(e.g.,
‘TAA’).
doi:10.1371/journal.pone.0173043.g002
Table).MeHDZ01andMeHDZ06,togetherwithMeHDZ37andMeHDZ41,whichmayresult
fromdifferentalternativesplicing,weremappedtothesamelocionchromosome3and1,
respectively.Overall,anunequaldistributionofMeHDZgeneswasrevealedincassava.For
example,amaximumofeightMeHDZgeneswerepresentedonchromosome1,incontrast,
onlyoneMeHDZgenewasrespectivelylocatedonchromosome7,10and18,whichcontained
theminimumMeHDZgeneincassava(Fig4).
Segmentalduplicationandtandemduplication,whichplayedimportantrolesinexpanding
newmembersduringtheevolutionofagenefamily,wereinvestigatedtoelucidatethepoten-
tialevolutionmechanismofMeHDZgenesincassava.Totally,tenMeHDZgenes,formingfive
pairs,wereidentifiedassegmentalduplications(Fig4;S3Fig).Ofwhich,twopairswerefrom
MeHDZIIsubfamily,whiletheotherthreepairswerefoundintheremainingthreesubfami-
lies,respectively.However,noneoftheMeHDZgenesseemedtobegeneratedfromtandem
duplicationsinouranalysis.Theseresultsimplicatedthatsegmentalduplicationeventswere
themaincausefortheexpansionofMeHDZgenesincassava.
PLOSONE|DOI:10.1371/journal.pone.0173043 March1,2017 7/20
HD-Zipgenefamilyincassava
Fig3.ConservedmotifsofMeHDZproteinscorrespondingtotheirphylogeneticrelationships.The
conservedmotifswereidentifiedbyMEMEsoftware.Motifswereindicatedbydifferentcoloredboxeswiththe
motifnumber,whilenon-conservedsequenceswererepresentedbygreylines.Lengthofmotifswas
exhibitedproportionally.
doi:10.1371/journal.pone.0173043.g003
Cis-elementidentificationofMeHDZgenes
TodemonstratethepossibleregulatorymechanismofMeHDZgenes,1.5kbpromoter
sequencesupstreamofthetranscriptionstartsite(ATG)wereextractedandusedforcis-ele-
mentprediction.Asexpected,manystress-relatedelements,suchasMBSfordrought,HSEfor
heat,LTRforlow-temperature,andWUN-motifforwound-response,wereidentifiedinthe
promoterregionofMeHDZgenes(S7Table).Additionally,manyhormone-relatedcis-ele-
ments,includingABREforABA,TGA-elementforauxin,EREforethylene,P-boxand
GARE-motifforgibberellicacid(GA),CGTCA-motifandTGACG-motifformethyljasmo-
nate(MeJA),andTCA-elementforsalicylicacid,werealsofound(S7Table).Moreover,other
cis-elementssuchascircadianforcircadiancontrol,AREforanaerobicinduction,andas
muchas20cis-elementsforlightresponsewerealsoobserved.
ToexploretheparticularrolesofMeHDZgenesineachsubfamily,overrepresentationanal-
ysiswasconductedrespectivelyforthesecis-elements,andsubfamily-specificoverrepresented
mannerswererevealed(S7Table).Severalcis-elementsrelatedtoabioticstress,e.g.,ABRE,
PLOSONE|DOI:10.1371/journal.pone.0173043 March1,2017 8/20
HD-Zipgenefamilyincassava
HSE,andTC-richrepeats,weresignificantlyoverrepresentedinMeHDZIandIIsubfamilies,
whileGARE-motifandMBSwereoverrepresentedinMeHDZIIIandIVsubfamilies,respec-
tively.Further,lightresponsiveelements(e.g.,Sp1andG-box)alsoexhibitedsubfamily-spe-
cificoverrepresentedmanners.Together,theseresultssupportedthatsubfamily-specificroles
ofHD-Zipsubfamiliesmightberelatedtotheregulatoryelementsinthepromoterregion.
ExpressionprofilesofMeHDZgenesindifferenttissuesofthree
genotypes
TorevealexpressionprofilesofHD-Zipgenesindifferenttissues,previousgeneratedRNA-seq
forthesamplesofstem,leaf,earlystorageroot(ESR),middlestorageroot(MSR),andlaststor-
ageroot(LSR)inonewildsubspecies(W14)andtwocultivars(Arg7andKu50)wasused.
ExpressedgeneswerearbitrarilydefinedasthosewithFPKM(FragmentsPerKilobaseofexon
perMillionfragmentsmapped)>1.Overall,differentexpressionprofileswererevealedwithin
eachsubfamilyofMeHDZgenes(Fig5A;S8Table).InMeHDZIsubfamily,mostgeneswere
expressedinalltestedorgans,ofwhichMeHDZ22andMeHDZ23showedthehighestexpres-
sionlevel.GenessuchasMeHDZ21,-47,-49and-54wereexpressedonlyinafewsamples,
whiletheotherssuchasMeHDZ26,-41and-51werenotexpressedatall.Similargeneexpres-
sionpatternswereobservedinMeHDZIIsubfamily.Forexample,MeHDZ27,-28,-31,-32,
-40and-44wereexpressedinalltestedsamples,whileMeHDZ56werenotexpressedand
theremaininggeneswereexpressedinpartialsamples.ComparedwithMeHDZIandIIsub-
families,mostgenesinMeHDZIIIandIVsubfamiliesexhibiteddistincttissue-specificexpres-
sionpatterns.Forexamples,eightoutofninegenes(exceptMeHDZ06thatwasnotexpressed
inanysamples)fromMeHDZIIIsubfamilywerelowerexpressedinleavesthaninothertis-
sues;andsevenoutofelevengenes(exceptMeHDZ10thatshowedconstitutiveexpression,
MeHDZ17andMeHDZ46thatwerenotexpressed,andMeHDZ43thatexpressedhigherin
rootthaninothersamples)fromMeHDZIVsubfamilywerelowerexpressedinrootsthanin
Fig4.PhysicallocationsofMeHDZgenesoncassavachromosomes.Chromosomenumberswere
indicatedatthetopofeachchromosome.MeHDZgenesfromdifferentsubfamilieswereindicatedbydifferent
colors.FivepairsofMeHDZgenesconnectedbydashlineswereresultedfromsegmentalduplication.
doi:10.1371/journal.pone.0173043.g004
PLOSONE|DOI:10.1371/journal.pone.0173043 March1,2017 9/20
HD-Zipgenefamilyincassava
Fig5.ExpressionprofilesofMeHDZgenesincassava.(A)ExpressionofMeHDZgenesindifferent
tissuesofthreegenotypes.Log2-transformedFPKMvaluewasusedtoplottheheatmap.ESR:earlystorage
root;MSR:middlestorageroot;LSR:laststorageroot.(B)ExpressionofMeHDZgenesinresponseto
droughtandPEGtreatments.Log2basedfoldchangeoftreatment/controlwasusedtoplottheheatmap.FL:
foldedleaf;FEL:fullexpandedleaf;BL:bottomleaf;RT:root.Thenumbersattachedbehindsamples
representedthetimepointsatwhichsampleswerecollected:e.g.,03and24represented3and24h,
respectively.
doi:10.1371/journal.pone.0173043.g005
others(Fig5A).Takentogether,theseresultssuggestedthatthefunctionsofMeHDZgenes
weredivergedthroughconstitutiveexpression,non-expression,ortissue-specificexpression
duringcassavaevolutionprocesses.
Differentexpressionpatternswerealsorevealedinthesametissueamongdifferentgeno-
types,especiallybetweenthecultivatedandthewildsubspecies.ComparedtoW14,theexpres-
sionof22,25and15HD-Zipgeneswasgreatlychanged(absolutevalueoffold-change>2)in
stem,leafandMSRinArg7,while15and17HD-Zipgenesweredifferentiallyexpressedinleaf
andMSRinKu50,respectively(S8Table).
ExpressionanalysisofMeHDZgenesrespondingtodroughtindifferent
genotypes
TodemonstratetheresponsesofMeHDZgenestodroughtstress,seedlingsofonewildsubspe-
cies(W14)andtwocultivatedvarieties(Arg7andSC124)weresubjectedto12-daywaterwith-
holdingtreatment,andtheirleafandrootsamplesofcontrolanddrought-treatedplantswere
PLOSONE|DOI:10.1371/journal.pone.0173043 March1,2017 10/20