Identification of Gene Networks for Residual Feed Intake in Angus Cattle Using Genomic Prediction

19 Pages·2016·1.35 MB·English

RESEARCHARTICLE
Identification of Gene Networks for Residual
Feed Intake in Angus Cattle Using Genomic
Prediction and RNA-seq
KristinaL.Weber1*,BryanT.Welly2,AlisonL.VanEenennaam2,AmyE.Young2,Laercio
R.Porto-Neto3,AntonioReverter3,GonzaloRincon1
1 VMRDGeneticsR&D,ZoetisInc.,Kalamazoo,MI,UnitedStatesofAmerica,2 DepartmentofAnimal
Science,UniversityofCaliforniaDavis,Davis,CA,UnitedStatesofAmerica,3 CSIROAgriculture,
QueenslandBiosciencePrecinct,St.Lucia,QLD,Australia
*[e-mail]
Abstract
Improvementinfeedconversionefficiencycanimprovethesustainabilityofbeefcattlepro-
duction,butgenomicselectionforfeedefficiencyaffectsmanyunderlyingmolecularnet-
OPENACCESS
worksandphysiologicaltraits.Thisstudydescribesthedifferencesbetweensteerprogeny
Citation:WeberKL,WellyBT,VanEenennaamAL,
oftwoinfluentialAngusbullswithdivergentgenomicpredictionsforresidualfeedintake
YoungAE,Porto-NetoLR,ReverterA,etal.(2016)
IdentificationofGeneNetworksforResidualFeed (RFI).Eightsteerprogenyofeachsirewerephenotypedforgrowthandfeedintakefrom
IntakeinAngusCattleUsingGenomicPredictionand 8mo.ofage(averageBW254kg,withameandifferencebetweensiregroupsof4.8kg)
RNA-seq.PLoSONE11(3):e0152274.doi:10.1371/ untilslaughterat14–16mo.ofage(averageBW534kg,siregroupdifferenceof28.8kg).
journal.pone.0152274
Terminalsamplesfrompituitarygland,skeletalmuscle,liver,adipose,andduodenumwere
Editor:RamonaNatachaPENAiSUBIRÀ,University
collectedfromeachsteerfortranscriptomesequencing.Geneexpressionnetworkswere
ofLleida,SPAIN
derivedusingpartialcorrelationandinformationtheory(PCIT),includingdifferentially
Received:November9,2015
expressed(DE)genes,tissuespecific(TS)genes,transcriptionfactors(TF),andgenes
Accepted:March12,2016 associatedwithRFIfromagenome-wideassociationstudy(GWAS).Relativetoprogenyof
Published:March28,2016 thehighRFIsire,progenyofthelowRFIsirehad-0.56kg/dfinishingperiodRFI(P=0.05),
-1.08finishingperiodfeedconversionratio(P=0.01),+3.3kg^0.75finishingperiodmeta-
Copyright:©2016Weberetal.Thisisanopen
accessarticledistributedunderthetermsofthe bolicmid-weight(MMW;P=0.04),+28.8kgfinalbodyweight(P=0.01),-12.9feedbunk
CreativeCommonsAttributionLicense,whichpermits visitsperday(P=0.02)with+0.60min/visitduration(P=0.01),and+0.0045carcassspe-
unrestricteduse,distribution,andreproductioninany cificgravity(weightinair/weightinair—weightinwater,apredictorofcarcassfatcontent;P
medium,providedtheoriginalauthorandsourceare
=0.03).RNA-seqidentified633DEgenesbetweensiregroupsamong17,016expressed
credited.
genes.PCITanalysisidentified>115,000significantco-expressioncorrelationsbetween
DataAvailabilityStatement:Sequencedatafrom
genesand25TFhubs,i.e.controllersofclustersofDE,TS,andGWASSNPgenes.Path-
theRNA-seqexperimentisavailablefromtheNCBI
SRAunderbioprojectPRJNA311009. wayanalysissuggestslowRFIbullprogenypossessheightenedgutinflammationand
reducedfatdeposition.Thismulti-omicsanalysisshowshowdifferencesinRFIgenomic
Funding:Fundingforthisstudywasprovided
throughagrantfromZoetisInc.,#201224918.The breedingvaluescanimpactothertraitsandgeneco-expressionnetworks.
funderprovidedsupportintheformofsalariesfor
authors[GR,KW],butdidnothaveanyadditional
roleinthestudydesign,datacollectionandanalysis,
decisiontopublish,orpreparationofthemanuscript.
Thespecificrolesoftheseauthorsarearticulatedin
the‘authorcontributions’section.
PLOSONE|DOI:10.1371/journal.pone.0152274 March28,2016 1/19
GeneNetworksforRFIinAngusCattle
CompetingInterests:KWandGRarecurrently Introduction
employedbythefundingagency,ZoetisInc.This
Thelargestvariablecostsinbeefcattleproductionarefeedandland[1].Feedcostscanbe
doesnotaltertheauthors'adherencetoPLOSONE
policiesonsharingdataandmaterials. reducedbyimprovingfeedconversionefficiencywithoutsacrificingproductiontraits;one
optionforthisistoconsiderresidualfeedintake(RFI;[2]).RFIisthedifferencebetween
observedfeedintakeandexpectedintakebasedonbodysizeandrateofproduction(growth,
milkproduction,etc.)overaperiodoftime.AnanimalwithalowRFIhasimprovedfeedcon-
versionefficiencysinceitconsumeslessfeedthanexpectedforitsmaintenanceandgrowth
requirements.TheheritabilityofRFIhasbeenestimatedtobeashighas42%ingrowingbeef
cattle[3],suggestingithasastronggeneticcomponentandwouldrespondtoselectivebreed-
ing.Whileitistooexpensiveforallbulltestingcentersandfeedlotstomeasureindividualfeed
intakeinordertocalculateRFI,referencepopulationsofanimalspossessingbothRFIpheno-
typesandhighdensitySNPgenotypeshavebeenusedtogeneratepredictionsoftotalgenetic
merit[4]forRFIinmajorcattlebreeds[5–7].Thesegenomicpredictionscanbeaselection
toolintheabsenceofdirectphenotypicdata.
Theaimofthisstudywastodeterminewhetherprogenyofbullswithdivergentgenomic
predictionsforRFIdisplayedphysiologicalandtranscriptomicdifferences.Weusedtwoinflu-
entialAngusbullswithhighorlowgenomicallypredictedbreedingvaluesforRFItoproduce
steerprogeny,whichwerephenotypedforRFIandalargevarietyoftraitsfrompost-weaning
throughslaughterandcarcassevaluation.Keytissuesforgrowthandmetabolism,including
thepituitary,skeletalmuscle,liver,visceraladipose,andduodenum,werecollectedforRNA-
seqanalysis.Partialcorrelationandinformationtheory(PCIT)wasusedtogenerategenenet-
works,linkingtissuespecific(TS)anddifferentiallyexpressed(DE)geneswithknowntran-
scriptionfactors(TF)andgenesharboringSNPfromalargegenome-wideassociationstudy
(GWAS)study.TheDEgenesandgenenetworksdescribedhereprovidenewinsightsintoreg-
ulatoryandgeneexpressionpathwaysofRFI.
Materials&Methods
AnimalManagementandPhenotypeCollection
UsinggenomicpredictionsforRFIinAnguscattle[8]producedbyZoetisInc.(FlorhamPark,
NJ),twoinfluentialpurebredAngusbullswereselectedasbreedingsires.Onebullwaspredicted
tobeinthetop1%intheAngusbreedin2010forRFI,andtheotherinthebottom10%,witha
differenceinRFIbreedingvalueof0.32kgDM/dorapproximately3.5%ofdailyintake.Semen
fromthesebullswasusedtoinseminateagroupofcommercialcowsofpredominantlyAngus
backgroundattheUniversityofCaliforniaSierraFoothillResearchandExtensionCenter
(BrownsValley,CA).Eightsteeroffspringpersirewereselectedforparticipationinthestudy,
matchedbetweensiregroupsbyweaningweight.Eightanimalspergroupprovides80%power
todetectdifferentialgeneexpressionwithfoldchangegreaterthanorequalto1.8forthetop
80%ofgenesbycoverageusingabiologicalcoefficientofvariationof0.4andalphaof0.05([9]).
Theuseofhalf-siblingprogenyoftwobullsisexpectedtoreducethevariabilitywithinlowor
highRFIpopulations,butdoesconfoundRFIwithothertraitsforwhichthebullsdiffer.Siredif-
ferencesforotherbreedingvaluesareprovidedinS1Table.Basedonmolecularandtraditional
breedingvalues,thelowRFIsireispredictedtohavelowerintake,higherweightsfrombirth
throughtocarcass,largerheightandrib-eyearea,andreducedbackfatthickness.
AlllivestockweremanagedunderprotocolsapprovedbytheInstitutionalAnimalCareand
UseCommittee(IACUC)oftheUniversityofCalifornia.Postweaning,at~8moofage,the
animalsweretransferredtoacommercialfeedlot(SnyderLivestockCo.,Yerington,NV;http://
www.slcnv.com/),wheretheywerehousedinasinglepenequippedwithGrowSafeunits
PLOSONE|DOI:10.1371/journal.pone.0152274 March28,2016 2/19
GeneNetworksforRFIinAngusCattle
(GrowSafeSystemLtd.,Airdrie,AB,Canada)fordailyfeedintakemeasurementoveraperiod
of70d(henceforthreferredtoasthegrowingperiod)aftera14dadaptationperiod.Animals
werefedastandardgrowerdiet(Table1).Feedwasprovideddailyinthemorningandpushed
towardsthecattletwoadditionaltimesduringthedaytoencouragefeeding.Feedintakewas
measuredcontinuouslytothe0.01kgbytheGrowSafesystem.Datafrom2dwereomitteddue
tosystemfailure.Feedingbehaviorwasassessedbasedonfrequencyanddurationofbunkvis-
its.Bodyweightsweretakenevery14d.
Followingthegrowingperiod,theanimalsweretransferredtotheUniversityofCalifornia,
DavisAnimalScienceFeedlotfacility,wheretheywereindividuallyhousedin7.5m2(1.5x5
m)pensforaminimumof70daftera14dadaptationperiodoruntiltheyreachedabodycon-
ditionappropriateforslaughter(91donaverage),hereafterreferredtoasthefinishingperiod.
Totalfeedofferedandamountrefusedweremeasuredto0.05kgprecision,anddailyintake
wascalculatedasthedifferencebetweenofferedandrefusedfeedperanimalperday.Body
weightsweremeasuredevery14dbeforeAMfeeding.
Forfiverandomlyselecteddaysduringthefinishingperiod,thedailyrationwasprovided
viaatie-stallGreenFeedunit(C-LockInc.,RapidCity,SD)inordertocollectmethane
Table1. Compositionofthetotalmixedration(TMR)usedinthegrowingandfinishingperiods. Fin-
ishingdietpossessedahigherproportionofgrainandhigherenergycontent(Mcal/kg).
Item Growing Finishing,Mean(SE)
Ingredient,%
Cornsilage 32.56
Cornsteepwater 2.75
Rolledcorn 62.65
Distillersgrain 17.23
Almondhulls 20.00
Ricebran 18.28
Wheat 10.00
Alfalfahay 8.00 7.83
Oathay 2.00 3.91
Molasses 4.74
Condensedmolassessolubles 2.89
Whey 2.30
Fat 1.96
Limestone 1.28
Calciumcarbonate 0.50
Whitesalt 0.46 0.26
Suspensionagent 0.18
Magnesiumoxide 0.13
Traceminerals 0.06
Rumensin 0.02 0.01
Analyzed,%DMbasis
Crudeprotein 11.36 13.30(0.86)
Crudefiber 16.35
ADF 9.00(1.43)
NDF 16.80(1.34)
Calculated
NEm[10],Mcal/kg 1.263 2.112(0.047)
NEg[10],Mcal/kg 0.723 1.437(0.039)
doi:10.1371/journal.pone.0152274.t001
PLOSONE|DOI:10.1371/journal.pone.0152274 March28,2016 3/19
GeneNetworksforRFIinAngusCattle
emissionmeasurementsduringfeeding.Foroneofthosedays,thetimingandquantityoffeed
consumedfromtheGreenFeedunitwascontrolledtominimizevariationduetoanimalfeed-
ingbehavior.Duringthatday,approximately1kgoftotalmixedration(TMR)wasprovidedat
alternating1.5and3hintervalsinordertomeasurethepeaksandtroughsofmethaneproduc-
tion[11];ontheotherdays,feedprovidedfromtheGreenFeedunitwasapelletedformofthe
finishingperiodTMRcomponents,providedadlibitum.Fornon-GreenFeeddays,dailyration
wasdividedintofourpartsofferedthroughouttheday(earlyandlatemorningandafternoon)
inorderfortheintermittentfeedingallowedbytheGreenFeedtonotpresentahugedeparture
fromhabituatedfeedingpatterns.
Forboththegrowingandfinishingperiods,RFIisdefinedastheresidualofdrymatter
intakeadjustedforrateofgainandbodysize[12],
DMI¼b þb ADGþb MMWþRFI
0 1 2
whereDMIisaveragedailydrymatterintake;ADGisaveragedailygain,estimatedfromthe
regressionof14d-bodyweightsontime;MMWismetabolicmid-weight,or(meanperiod
bodyweight)0.75;b ,b andb arethepartialregressioncoefficientsassociatedwiththeinter-
0 1 2
cept,ADGandMMWT,respectively.
Atslaughter,standardbeefcarcasstraitphenotypeswerecollected,withtheadditionof
organweightsandcarcassspecificgravity,apredictorofthefatcontentofthecarcass[13],esti-
matedfromCW/(CW-CW )whereCWisthedrycarcassweightandCW istheweight
H2O H2O
ofthecarcasssuspendedinwater(measuredtonearest1g).
RNAIsolation,cDNALibraryConstruction,andTranscriptome
Sequencing
Samplesofpituitary,skeletalmuscle,liver,visceraladipose(KPHfat),andduodenumwerecol-
lectedwithin25minutesofslaughter,snapfrozeninliquidnitrogentopreserveRNAintegrity,
andfrozenat-80°C.Tocollectthepituitary,theskullwasbisectedlongitudinallyandthepitui-
taryidentifiedbyabovinephysiologist.ToextractRNA,approximately200mgoffrozentissue
wasimmersedinliquidnitrogen,groundwithamortarandpestle,andhomogenizedin2mL
TRIzol1(ThermoFisherScientific,Waltham,MA)usingaMiniBeadBeater-8(BiospecProd-
ucts,Bartlesville,OK)forupto10s,exceptfattissuewhichwasshakenratherthanbeadbeaten
toavoidRNAdegradation.TotalRNAwaspurifiedusingtheTRIzolstandardprotocol
(ThermoFisherScientific,Waltham,MA),andresuspendedin25μLofRNase-freewater.
RNAquantityandqualitywereassessedusingNanoDrop1spectrophotometer(ThermoSci-
entific,Wilmington,DE)andAgilent2100Bioanalyzer(SantaClara,CA).Foreachtissuetype,
mean260/280ratiorangedfrom1.84to1.93(0.10SD),andmeanRNAintegritynumberfrom
7.0to8.1(2.0SD).TheIlluminaTruSeqstrandedmRNAHTsamplepreparationprotocol
(SanDiego,CA)wasusedtogeneratedual-indexedcDNAlibrariesforeachsampleusing1μg
oftotalRNAasinput.SuccessoflibrarypreparationwasdeterminedbyAgilent2100Bioanaly-
zer(SantaClara,CA),afterwhichall80samplesweremultiplexedintoasingleequimolarpool
toavoidsequencingrunbias.Singleend100bpsequencingwasconductedonanIllumina
HiSeq2000followingstandardprotocols.
SequenceAnalysisandMixedModelAnalysis
Eightysamples,fivetissuesfromsixteenAngussteers,wereanalyzedwithRNA-seq.Failed
readswerediscardedbeforeanalysis.Sampleswithlessthan10millionreads(oneadiposeand
oneduodenalsample)werediscarded.Fortheremainingsamples,weobtained45.2Mreads
persampleonaverage,withmeanPhredscore35.8(0.13SD),medianPhred38,andlessthan
PLOSONE|DOI:10.1371/journal.pone.0152274 March28,2016 4/19
GeneNetworksforRFIinAngusCattle
7%ofreadsbelowPhredscore30.Noadditionalsequencefilteringwasperformed.Sequence
readsweremappedtothebovinereferencegenome(UMD3.1forBTA1-29andXandBtau
4.6.1forBTAY)usingCLCGenomicsWorkbench7RNA-seqtool(RedwoodCity,CA).Reads
wererequiredtomapwith80%similarityfor90%oftheirlength,which93%ofreadsachieved.
Weusedreadsperkilobaseofgenepermillionmappedreads,RPKM[14],astheunitof
expression.Agenewasconsideredtissuespecific(TS)ifonetissuerepresentedatleasttwo-
thirdsofthetotalexpressionacrosstissues.HighlyexpressedTSgeneswerecomparedwithtis-
sueexpressionprofilesinBostaurusandothervertebratespeciesbasedontheEMBL-EBI
ExpressionAtlas(https://www.ebi.ac.uk/gxa/home)toverifythattheywereconsistentwith
previouslyobservedexpressionpatternsinthesetissues.
Geneexpressionbysiregroupwasestimatedusingamixedmodelapproach[15,16]in
ordertodeterminewhetheranygenesweredifferentiallyexpressedacrosstissuesduetosire
effects.
Y ¼μþGþGT þGA þGS þe
ijkl i ij ik il ijkl
wherelog2-transformedRPKM(Y )wasmodeledasafunctionoftherandomeffectsofgene
ijkl
(G),genebytissue(GT ),genebyanimal(GA ),andgenebysiregroup(GS )forigenes,jtis-
i ij ik il
sues,kanimals,andlsires.Randomresidual(e )wasassumedtobeindependentandidenti-
ijkl
callydistributed.VariancecomponentanalysiswasperformedusingVCE6software(Eildert
Groeneveld,Friedrich-Loeffler-Institut,ftp://ftp.tzv.fal.de/pub/vce6/).Agenewasconsidered
differentiallyexpressedbetweensiregroups(DE)ifthegenebysireeffectwasatleasttwostan-
darddeviationsfromthemean,correspondingtoP<0.01.TocompareTSandDEgeneswith
knowntranscriptionfactors(TF),allknownTFfromAnimalTFDB[17](http://www.bioguo.
org/AnimalTFDB/)werefilteredforthosebothabundantandmostconsistentlyassociated
withTSandDEgenesbetweensiregroupsbasedonregulatoryimpactfactormetrics[16,18]
whichweightsaveragegeneexpressionacrosssamples(G),thedifferentialexpressiondueto
i
siregroup(GS ),andtheco-expressioncorrelationbetweenTFandDEorTSgenes.Genes
il
harboringSNPidentifiedbyBolormaaetal.[19]asassociatedwithRFI(hereafterreferredto
asSNPgenes)wereincludedtodeterminehowcloselythisexperimentalignswithaGWAS
studyfromabroaderbeefcattlepopulation.
NetworkandPathwayAnalysis
AgenenetworkforRFIwithinthesefivetissueswasderivedusingtheTS,DE,TFandSNP
genesasnodes,andsignificantconnectionsbetweenthemweredeterminedviapartialcorrela-
tionandinformationtheory(PCIT)algorithmsandsoftwareasdescribedbyReverterand
Chan[20].PCITreliesoncalculatingthecorrelationbetweenapairofgenesafteraccounting
forallothergenes(thus,thepartialcorrelation).Connectionsbetweengenenodeswere
acceptedwhenthepartialcorrelationwasgreaterthantwostandarddeviationsfromthemean
(P<0.01).Pairwiseco-expressedgeneswereimportedintoCytoscapesoftware[21](http://
www.cytoscape.org/),tobevisualizedandclusteredintonetworksusingaforce-directededge-
weightedspringembeddedlayout[22].FurtherfunctionalanalysiswasperformedusingInge-
nuityPathwayAnalysis1(RedwoodCity,CA;www.qiagen.com/ingenuity)andBlast2GO
PRO[23].ForIPAanalyses,thehuman/mouse/ratgeneorthologueIDsandthefoldchanges
forallDEbovinegenesexpressedineachtissuewerecomparedtotheIPAdatabaseofpath-
ways,networks,diseases,functions,andregulatorsforgenesetenrichment(P<0.05)andacti-
vationZ-score(|Z|>2)analysis.ForBlast2GO,allannotatedbovinegeneswereimported
fromEnsemblBiomart(http://www.ensembl.org/biomart/martview/),andsignificantly
enrichedGOterms(P<0.05)foragivengenesubsetwereidentifiedusingFisher’sexacttest.
PLOSONE|DOI:10.1371/journal.pone.0152274 March28,2016 5/19
GeneNetworksforRFIinAngusCattle
Results
SireGenomicBreedingValuesPredictedDifferenceinProgeny
Phenotypes
ThehighandlowRFIbulls’progenywerephenotypedforRFIandavarietyofgrowth,intake,
metabolism,andbodycompositiontraitstodeterminewhichtraitsdifferentiatedthem.Aspre-
dictedfromtheirsires’breedingvalues,theprogenyofthelowRFIsirehadlowerRFIthanthe
progenyofthehighRFIsireinthefinishingphase(-0.57kgDM/d),withasimilartrendinthe
growingphase(-0.62kgDM/d;Table2).AswithRFI,feedconversionratio(DMI/ADG)was
decreasedinprogenyofthelowRFIsirerelativetothehighRFIsireinthefinishingphase.Sire
breedingvaluespredictprogenyofthelowRFIsiretoexhibitreducedintakewithsimilarADG
relativetoprogenyofthehighRFIsire,whichwasobservedinthegrowingphase(siregroup
differencesof-0.48kg/dDMIand0.00kg/dADG)butnotinthefinishingphase(siregroup
differencesof0.01kg/dDMIand0.22kg/dADG).Largerfinishingandcarcassweightswere
observedintheprogenyofthelowRFIsire,asmightbepredictedfromthehigherpedigree-
basedbreedingvaluesofthatsireforweaningweight(WW),yearlingweight(YW),andmature
weight(MW).ThelowRFIsiregrouphadhighercarcassspecificgravities,indicatingthat
theseanimalswerelessfatorleanerforthesamecarcassweight.Nosignificantdifferences
wereobservedinothercarcasstraitsorinheart,liverorkidneyweights.
Table2. Phenotypicdifferencesbetweensiregroups(lowRFIsire–highRFIsire)comparedwithexpectedprogenydifferencebasedonbullbreed-
ingvalues. Differencesbetweensireprogenygroupsforgrowth,intake,feedefficiency,feedingbehavior,methaneproduction,andcarcasstraitsarepro-
vided,includingtrait,traitmean,meandifferencebetweensiregroups,P-valueforthatdifference(P<0.05bolded),andexpectedprogenydifferencebased
onsirebreedingvalueswhereavailable.
TraitType Trait Period Mean(SE) SireGroupDifference P-value ExpectedProgenyDifference
Growth, MMW,kg0.75 Growing 74.5(1.1) 0.8 7.3E-01
Intake,& Finishing 100.4(0.9) 3.3 4.0E-02
Feed ADG,kg/d Growing 1.71(0.05) 0.00 9.8E-01 0.00
Efficiency Finishing 1.40(0.06) 0.22 8.0E-02
DMI,kgDM/d Growing 9.55(0.25) -0.48 3.4E-01 -0.18
Finishing 8.52(0.20) 0.01 9.8E-01
Feedconversion Growing 5.63(0.17) -0.28 4.1E-01
ratio(DMI/ADG) Finishing 6.21(0.23) -1.08 1.0E-02
RFI,kgDM/d Growing 0.00(0.19) -0.62 9.0E-02 -0.16
Finishing 0.00(0.15) -0.57 5.0E-02
Feeding BunkVisits,visits/d 57.7(2.9) -12.9 2.0E-02
Behavior BunkVisitDuration,min/visit 1.87(0.13) 0.60 1.0E-02
(GrowingPeriod) FeedingDuration,min/d 96.7(3.1) 7.04 2.6E-01
MethaneProduction Controlledfeedamountand 165.3(13.1) 6.6 8.1E-01
interval,gCH /d
4
(FinishingPeriod) Adlibfeeding,gCH /d 161.8(9.5) 8.7 6.6E-01
4
Carcass Finalweight,kg 533.9(6.1) 28.8 1.0E-02 MW33.6
Hotcarcassweight,kg 323.4(4.1) 12.6 1.1E-01 10.5
Rib-eyearea,cm2 80.2(1.9) 6.8 7.0E-02 5.7
Backfatthickness,mm 12.9(0.4) -0.5 5.4E-01 -0.5
Yieldgrade 3.0(0.1) -0.3 2.6E-01
Marblingscore 6.8(0.3) 0.0 1.0E+00 0.0
Carcassspecificgravity 1.08(0.002) 0.0045 3.0E-02
Heartweight,kg 2.01(0.10) 0.25 2.0E-01
Liverweight,kg 6.00(0.17) -0.03 9.2E-01
Kidneyweight,kg 1.00(0.03) -0.03 6.4E-01
doi:10.1371/journal.pone.0152274.t002
PLOSONE|DOI:10.1371/journal.pone.0152274 March28,2016 6/19
GeneNetworksforRFIinAngusCattle
Whilefeedingbehaviorwasonlymeasuredinthegrowingperiod,progenyofthelowRFI
sirevisitedthefeedbunklessoften(-12.9visits/d)andstayedlonger(0.6min)ateachvisit.No
significantdifferenceswereobservedformethaneproduction.
DifferentialExpressionandPathwayAnalysis
Of24,737annotatedbovinegenes,7,721werenotexpressed(RPKM<0.2)inanyofthetissues
surveyed.Weidentified1,026TSgenes,including285inthepituitary,220inskeletalmuscle,
275intheliver,33inadipose,and213intheduodenum,and633DEgenesacrossalltissuesin
thelowversushighRFIsiregroups(S2Table).Pathwayanalysisshowedmostdifferentially
expresseddiseaseandbiofunctions(P<0.05and|Z|>2)correspondedtoactivationofthe
immuneresponseintheduodenumanddownregulationoffatdepositioninadiposeandmus-
cletissue(Table3).Thecombinedeffectofdownregulationofphosphoenolpyruvatecarboxy-
kinase2(PCK2),thyroidhormoneresponsive(THRSP),fattyacidsynthase(FASN),stearoyl-
CoAdesaturase(SCD),acetyl-CoAcarboxylasealpha(ACACA),acyl-CoAsynthetaselong-
chainfamilymember1(ACSL1),glycerol-3-phosphateacyltransferase(GPAM),and1-acylgly-
cerol-3-phosphateO-acyltransferase9(AGPAT9)supportadeactivationofaregulatorynet-
workcontrollingfattyacidmetabolismintheadiposetissueoftheprogenyofthelowRFIsire
(Fig1).Thisisconsistentwiththeincreaseincarcassspecificgravityoftheprogenyofthelow
RFIsire,whichsuggestsareducedbodyfatcontent.
Ofthe633DEgenes,122wereTS,6wereTF,and12harborSNPidentifiedbyBolormaa
etal.[19]asassociatedwithRFIinalargecattlepopulation(Fig2).The6DETFareETSvari-
ant4(ETV4),group-specificcomponent/vitaminDbindingprotein(GC,alsoliver-specific),
highmobilitygroupbox1protein(HMGB1),sexdeterminingregionY-box6(SOX6),transdu-
cinbeta-like1Y-linked(TBL1Y),andanuncharacterizedprotein(ENSBTAG00000036343).
ETV4ismosthighlyexpressedinthepituitaryandduodenumandisupregulatedintheprog-
enyofthelowRFIsire.HMGB1,SOX6,andTBL1Yareexpressedinallassayedtissuesbut
mosthighlyintheduodenum.HMGB1andTBL1Yareupregulatedinalltissues,andSOX6is
downregulatedintheduodenum,pituitary,andskeletalmuscleandupregulatedinliverand
adipose.ENSBTAG00000036343isanadipose-specificheatshockproteinandstronglydown-
regulated(-5.2foldchange).The12DEgenesharboringGWASSNPareacyl-CoAsynthetase
long-chainfamilymember1(ACSL1),carcinoembryonicantigen-relatedcelladhesionmole-
cule20(CEACAM20),chordin-like2(CHRDL2),claudin2(CLDN2),glypican3(GPC3),glu-
tamatereceptorionotropicAMPA2(GRIA2),interleukin2receptorgamma(IL2RG),lipoma
HMGICfusionpartner-like1(LHFPL1),lin-52DREAMMuvBcorecomplexcomponent
(LIN52),ryanodinereceptor3(RYR3),transmembraneprotein8C(TMEM8C),andZW10
interactingkinetochoreprotein(ZWINT).WhilerelativelyfewDEgenesareTForSNPgenes,
somehavekeybiologicalfunctionsintheimmunesystem,lipidmetabolism,andmusclepro-
cesses.GCistheprimarycarrierofvitaminDinthebloodandamacrophageactivatingfactor
whichisusedinsomecancertherapies[24].IL2RGisasignalingcomponentofmanyinterleu-
kinreceptorsandloss-of-functionofthisgenehasbeenassociatedwithX-linkedimmunodefi-
ciencyinhumans[25].ACSL1isanintermediateoffattyacidmetabolism.RYR3isinvolvedin
releasingcalciumfromintracellularstores.Overall,manyDEgeneswereTSandonlyafew
werealsoTForGWASSNPgenes,manysupportingimmunesystemandfatmetabolismroles.
PCITNetworkAnalysis
Partialcorrelationinformationtheorywasusedtodeterminesignificantpartialcorrelations
betweenTS,DE,TF,andSNPgenes(N=115,058edges;Fig3).
PLOSONE|DOI:10.1371/journal.pone.0152274 March28,2016 7/19
GeneNetworksforRFIinAngusCattle
Table3. PathwayanalysisofdifferentiallyexpressedgenesbetweenlowandhighRFIsiregroupsbytissue,forpituitary(P),skeletalmuscle(M),
liver(L),adipose(A),andduodenum(D). Topbiofunctionsby|Z-score|fromIPAareprovided,expressedaslowrelativetohighRFIsireprogenygroupdif-
ferences;allaresignificant(P<0.05).Functionsassociatedwiththeimmunesystemorfatdepositionaregroupedtogether,allothersignificantfunctionsat
thebottom.
DiseasesorFunctionsAnnotation P M L A D DiseasesorFunctionsAnnotation P M L A D
Immunesystem Fatdeposition
activationofantigenpresentingcells 2.6 differentiationofadipocytes 2.4
immuneresponseofantigenpresenting 2.0 massofadiposetissue 2.2
cells
inflammationofbodyregion 2.1 concentrationoftriacylglycerol -2.5
chemotaxisofgranulocytes 2.0 synthesisoftriacylglycerol -2.4 -2.6
recruitmentofgranulocytes 2.2 quantityofdiacylglycerol -2.0
inflammationoflargeintestine -2.7 concentrationofacylglycerol -2.1
colitis -2.7 synthesisofacylglycerol -2.2
activationofleukocytes 2.3 2.7 metabolismofmembranelipidderivative -2.2
immuneresponseofleukocytes 2.2 concentrationof -2.3 -2.3 -3.1
1,2-dipalmitoylphosphatidylcholine
recruitmentofleukocytes 2.7 2.8 2.8 concentrationofphosphatidylcholine -2.2 -2.1
activationofmacrophages 2.5 beta-oxidationoffattyacid -2.2
activationofmononuclearleukocytes 2.3 accumulationoflipid -3.0 -2.4
cellmovementofmononuclear 2.0 bindingoflipid 2.0
leukocytes
responseofmononuclearleukocytes 2.6 fluxoflipid 2.6 2.5
activationofmyeloidcells 2.3 transportoflipid 2.8
responseofmyeloidcells 2.4 transportofcholesterol 2.1
chemotaxisofneutrophils 2.3 fattyacidmetabolism 2.1
activationofphagocytes 2.4 transportofsteroid 2.2
immuneresponseofphagocytes 2.4 steroidmetabolism -2.2
recruitmentofphagocytes 2.5 2.1 2.4 concentrationofhormone -3.0
responseofphagocytes 3.0 2.2 bindingofcarbohydrate 2.0 2.2
viralinfection -2.3 uptakeofcarbohydrate 2.1
inflammatoryresponse 2.3 3.2 insulinresistance 2.1
Othertraits
morbidityormortality -2.2 relaxationofartery 2.7
necrosis -2.5 recruitmentofbloodcells 2.4 2.6 2.6
depositionofextracellularmatrix -2.2 -2.2 vasoconstriction 2.0
contractilityofmuscle 2.2 bleeding -2.3
proliferationofmusclecells 2.6 functionofcardiacmuscle 2.1
formationoffibrils -2.2 activationofcells 2.3 2.0 2.2
fibrogenesis -2.0 2.2 cellmovement 2.1 2.6
injuryofliver 2.1 engulfmentofcells 2.6
exportofmolecule 2.3 migrationofcells 2.1 2.8
neurologicalsigns 2.2 damageofendothelialcells 2.0 -2.2
secretionofprotein 2.2 damageofepithelialtissue 2.5
softtissueneoplasm -2.1 stimulationofepithelialcells -2.1
doi:10.1371/journal.pone.0152274.t003
Therewere6majorclusterstothePCITnetwork,whichalignwiththetissueofhighestrela-
tiveexpression(withmusclehavingtwoclustersandallothertissueshavingoneclusterpertis-
sue),asmorethan90%ofedgesarebetweengeneswithintissue.Thenumberofedgesfrom
onetissuetoanother(e.g.fromlivertoanyothertissue)wereproportionaltothenumberof
nodesforthattissue(R2=0.94).Commensuratewithitsroleasahormonalregulator,ofthe
PLOSONE|DOI:10.1371/journal.pone.0152274 March28,2016 8/19
GeneNetworksforRFIinAngusCattle
Fig1.DownregulationofregulatorynetworkcontrollingfatdepositionamonglowRFIsireprogeny.
Thisfiguredisplaystheregulatorynetworkconnectingconcentrationoftriacylglycerol,accumulationoflipid,
quantityofinsulinintheblood,quantityofconnectivetissue,andmigrationofphagocytes,overlaidwithfold
changesobservedinadiposetissue.Red=up-regulated,green=down-regulated,orange=predicted
activation,blue=predictedinhibition,andyellow=inconsistent.
doi:10.1371/journal.pone.0152274.g001
top50genesrankedbynumberofedgestoothergenes,38aremosthighlyexpressedortissue
specifictothepituitary,suchasluteinizinghormone/lutropinsubunitbeta(LHB;280edges),
Gprotein-coupledreceptor173(GPR173;273edges),EF-handdomainC-terminalcontaining
1(EFHC1;268edges),prolactin(PRL;267edges),andgrowthhormone/somatotropin(GH1;
266edges).ThedistributionofnumberofedgespergeneisprovidedinFig4.
ThemajorityofedgesarebetweenTSgenes(54%)orbetweenTSandDE,TF,orSNPgenes
(32%).ThereareproportionallymoreedgestoDEgenesthantoSNPgenes(27%versus20%).
TwogenesidentifiedinBolormaaetal.[19]aspotentialQTLforRFIwereSrchomology2
domaincontainingtransformingprotein3(SHC3)andinsulin-likegrowthfactorbindingpro-
tein2(IGFBP2).SHC3,whilenotDEoraTF,wassignificantlycorrelatedwith18DEgenes
(predominantlyexpressedinthepituitaryandduodenum),17SNPgenes,5TF[Mshhomeo-
box1(MSX1),pairedbox8(PAX8),zincfingerprotein300(ZNF300),v-etsavianerythroblas-
tosisvirusE26oncogenehomolog1(ETS1),andsignaltransducerandactivatorof
transcription6,interleukin-4induced(STAT6)],and1TSgene.IGFBP2wasadjacenttobut
didnotharborGWASSNP,wasnotusedasaSNPgeneinthisstudynorwasitidentifiedas
DEorTS,butrelatedproteinsIGFBP1and4wereTSandDEandcorrelatedwith219and18
genesintheliver,respectively.FivepercentofedgesconnecttoTF,ofwhichthemosthighly
connectedTFareMSX1(235edges),GC(228),v-mybavianmyeloblastosisviraloncogene
homolog(MYB;203),ETS1(200),SIXhomeobox3(SIX3;194),SAMpointeddomaincontain-
ingETStranscriptionfactor(SPDEF;192),zincfingerE-box-bindinghomeobox1(ZEB1;
187),LIMhomeobox3(LHX3;182),andmyogenicdifferentiation1(MYOD1;165).TFand
genesconnectedtothembyanedgewereextractedandvisualizedasanewnetwork(Fig5).
IntheTFnetworkdiagram,theclustersbytissuearepresent,withsmallclustersaroundTF
thatregulategeneexpressionwithinandacrosstissues.TheTFwiththegreatestnumberof
PLOSONE|DOI:10.1371/journal.pone.0152274 March28,2016 9/19
GeneNetworksforRFIinAngusCattle
Fig2.VenndiagramofDE,TS,TF,andSNPgenes.
doi:10.1371/journal.pone.0152274.g002
edgesinthepituitaryare,indescendingorder,MSX1,ETS1,SIX3,ZEB1,LHX3,myelintran-
scriptionfactor1(MYT1),zincfingerE-boxbindinghomeobox2(ZEB2),friendleukemia
integration1transcriptionfactor(FLI1),andGATAbindingprotein3(GATA3).GATA3is
linkedthroughDEandSNPgenestoE2Ftranscriptionfactor1(E2F1)intheduodenum,in
whichtheothermajorTFareMYB,SPDEF,andcaudaltypehomeobox1(CDX1).Similarly,
ZEB1andFLI1arelinkedwithhematopoieticallyexpressedhomeobox(HHEX)intheliver,in
whichGC,uncharacterizedproteinENSBTAG00000030470,signaltransducerandactivatorof
transcription3(STAT3),andE2Ftranscriptionfactor4(E2F4)areothermajorTF.ZEB2is
linkedwithT-cellacutelymphocyticleukemia1(TAL1),theonlymajoradiposeTF,which
clusterswithmuscleTFsMYOD1,transcriptionterminationfactor,RNApolymeraseI(TTF1),
ZincfingerX-chromosomalprotein(ZFX),iroquoishomeobox3(IRX3),zincfingerand
SCANdomaincontaining21(ZSCAN21),andzincfingerprotein35(ZNF35).Mostclusters
areanapproximatelyevenmixtureofDEandSNPgenes,withtheexceptionoftheGC/
ZNF160-centricclusterintheliverwhichisricherinDEgenes,andaclusterofSNPgenes
linkedtoZFX,IRX3,andzincfingerprotein35(ZNF35)inskeletalmuscle.Comparingthese
PLOSONE|DOI:10.1371/journal.pone.0152274 March28,2016 10/19