Journal of Biogeography. 2019;46:1163–1177. | 1163 wileyonlinelibrary.com/journal/jbi Received:17July2018 | Revised:3February2019 | Accepted:17March2019 DOI: 10.1111/jbi.13576 RESEARCH PAPER Multi‐variable approach pinpoints origin of oak wood with higher precision Linar Akhmetzyanov1 | Allan Buras1 | Ute Sass‐Klaassen1 | Jan den Ouden1 | Frits Mohren1 | Peter Groenendijk2,3 | Ignacio García‐González2 ThisisanopenaccessarticleunderthetermsoftheCreativeCommonsAttribution-NonCommercialLicense,whichpermitsuse,distributionandreproduction inanymedium,providedtheoriginalworkisproperlycitedandisnotusedforcommercialpurposes. ©2019TheAuthors.Journal of BiogeographyPublishedbyJohnWiley&SonsLtd Editor:Dr.CarinaHoorn 1ForestEcologyandForestManagement Group,WageningenUniversity&Research, Wageningen,TheNetherlands 2DepartamentodeBotánica,Escola PolitécnicaSuperiorde Enxeñaría,UniversidadedeSantiagode Compostela,Lugo,Spain 3DepartmentofPlantBiology,Institute ofBiology,UniversityofCampinas– UNICAMP,Campinas,SP,Brazil Correspondence LinarAkhmetzyanov,ForestEcologyand ForestManagementGroup,Wageningen University&Research,POBox47,6700AA Wageningen,TheNetherlands. Email:
[email protected] Funding information FP7People:Marie-CurieActions,Grant/ AwardNumber:607545 Abstract Aim: Spatialvariationsofenvironmental conditions translate intobiogeographical patternsoftreegrowth.Thisfactisusedtoidentifytheoriginoftimberbymeansof dendroprovenancing.Yet,dendroprovenancingattemptsarecommonlyonlybased onring-widthmeasurements,andlargelyneglectadditionaltree–ringvariables.We explorethepotentialofusingwoodanatomyasadendroprovenancingtool,andinvestigatewhetheritincreasestheprecisionofidentifyingtheoriginofoakwood. Sincedifferenttree–ringvariablesholddifferentinformationonenvironmentalconditions prevailing at specific times of the growing season—which vary between source regions—we hypothesize that their inclusion allows more precise dendroprovenancing. Location: Europe,Spain. Taxon: Quercus roburL.,Quercus petraea(Matt.)Liebl.,Quercus fagineaLam.,Quercus pyrenaicaWilld. Methods: WesampledfouroakspeciesacrossNorthernSpain,i.e.fromtheBasque countryandCantabriaand—intheBasquecountry—fromlowtohighelevation(topographic/latitudinalgradient).Wemeasuredmultipletree–ringvariablesto(a)extract complementaryvariables;(b)presentstatisticalrelationsamongthem;(c)analyseregion-specificvariationintheirpatternsbasedontime–seriesofindividualtrees;and (d)determineunderlyingclimate–growthrelationships.Leave-one-outanalysiswas usedtotestwhetheracombinationofselectedvariablesalloweddendroprovenancingofarandomlyselectedtreewithinthearea. Results: Acombinationoflatewoodwidth(LW)andearlywoodvesselsizewasused topinpointtheoriginofoakwoodwithhigherprecisionthanringwidthorLWonly. VariationinLWpinpointedthewoodtoeastandwestareas,whereasvariationin vesselsassignedwoodtolocationsalongalatitudinal/topographicgradient.Theclimatictriggersbehindthesegradientsarerespectivelyaneast–westgradientinJune– July temperature and a north–south gradient in winter/spring temperatures. The leave-one-outanalysessupportedtherobustnessoftheseresults.
1164 | AKHMETZYANOV ET Al. 1 | INTRODUCTION Trees are affected by a set of spatially and/or temporarily varying environmentalfactors,suchasprecipitation,temperature,competition,soiltype,andmanyothers.Thecombinedeffectoftheselocal environmentalfactorsontreegrowthis reflectedin thewidthand structureofthetreerings,andcanbeusedtocharacterizeagiven species growing under specific site conditions (Babst etal., 2013; Holz, Hart, Williamson, Veblen, & Aravena, 2018; Martin-Benito & Pederson,2015;Vicente-Serrano,Camarero,&Azorin-Molina,2014). IndryMediterraneanareas,wheretreeshavelimitedaccesstoground water, drought-induced water shortage was observed as the main growth-limitingfactor(deLuisetal.,2011;Martin-Benito,Beeckman, &Cañellas,2013;Olanoetal.,2014).Underincreasedvapour–pressuredeficit treestranspiremore,whichleadstoafasterdepletion of plant-available soil moisture. A decreasing soil–water potential initiates stomatal closure in the leaves to prevent drought-induced mortality(McDowelletal.,2008).Theassociatedreductionofphotosynthesislimitstreegrowthandleadstotheformationofnarrowtree rings(Bréda,Huc,Granier,&Dreyer,2006).Incontrast,athigherlatitudesandelevations,temperaturehasbeenfoundasthemainfactor limitingtreegrowth(Grudd,2008;Körner&Paulsen,2004;McCarroll etal.,2013;Yang,He,Melvin,Zhao,&Briffa,2013).Lowtemperatures limitthetrees’photosyntheticrate,whichleadstoreducedgrowth rates(Petit,Anfodillo,Carraro,Grani,&Carrer,2011).Treesgrowing undermoderateconditions,e.g.inCentralEuropeanlowlandsorin Atlanticregions,exhibitmorecomplexrelationshipsbetweenenvironmentandgrowthwithoftenmultiplegrowth-limitingfactors(Gričar, de Luis, Hafner, & Levanič, 2013; Souto-Herrero, Rozas, & García- González,2018a)buttheirgrowthpatternsarestillindicativeforthe specificenvironment–growthinteractionsatagivensitethroughtime. Environmental information on growing conditions as stored in treeringsiswidelyusedtopredictgrowthresponsesoftreestoclimatechange(Holzetal.,2018;Williams,Michaelsen,Leavitt,&Still, 2010),tomodelpotentialfuturespeciesdistribution(Gutiérrez,Snell, &Bugmann,2016;vanderMaatenetal.,2017),toquantifyforest resistancetodroughtevents(Gazol,Camarero,Anderegg,&Vicente- Serrano, 2017) and to reconstruct growth-limiting climate factors during the pre-instrumental period (Čufar, De Luis, Zupančič, & Eckstein,2008;Rydvaletal.,2017;Schofield,Barker,Gelman,Cook, &Briffa,2016).Informationonspatialvariationofclimaticsignalsin treeringsisalsoimportantforarchaeologicalandhistoricalstudies, asitservesasabasisfordatinganddeterminingthegeographical originofhistoricaltimbers(“dendroprovenancing”)(Speer,2010).This canbeusedtoidentifysourceareasofwood,aswellasforreconstructing wood–trade networks (Domínguez-Delmás etal., 2015). Nowadays,dendroprovenancingisbecomingmoreandmoreimportantalsointermsofdetectionofillegallogging(Kagawa&Leavitt, 2010;Vlametal.,2018).Moreover,ancienttimbercontainsinformationaboutpastforest-managementpractices(Haneca,Boeren,Van Acker,&Beeckman,2006),andreflectstechnologicalchoices,such asselectionofspeciesforparticularelementsofshipsorbuildings, andalsotypesoftoolsusedforwoodprocessing(Billamboz,2003; Castro,2008;Tegel,Elburg,Hakelberg,Stäuble,&Büntgen,2012). Dendroprovenancing was successfully applied to identify the originoffoundationpilesunderhistoricalbuildingsintheNetherlands (Sass-Klaassen,Vernimmen,&Baittinger,2008)orwoodfromshipwrecks(Bridge,2011;Domínguez-Delmás,Nayling,Ważny,Loureiro, &Lavier,2013),andtounraveltheBalticwoodtrade(Ważny,2002). Dendroprovenancingofwoodfromshipwrecksisconsideredasthe mostchallenging,sinceboatsandshipsmaybefoundfarawayfrom theconstructionplace(Bridge,2012),whichisnotnecessarilyrelated tothetimbersourceareasand,singleshipwrecksmaycontainwood fromdifferentorigins(Haneca,Čufar,&Beeckman,2009). Until now,dendroprovenancingis almostexclusivelybasedon correlationanalysesbetweenring-widthseriesandasetofreferencering-widthchronologies.Astatisticallysignificantmatchleads to dating and provenancing of a given ring-width sequence (indicatedbytheareacoveredbythereferencechronology)(Fowler& Bridge,2017).Insomespecificareas,stableisotopicratioshavealso showngreatpotentialtobeusedinprovenancingstudies(Bridge, 2012). In combination with dendrochronological methods, carbon,oxygenandsulphurisotopeshavebeenusedtoimprovethe spatial resolutionof provenancing(e.g.Kagawaand Leavitt2010; Horacek2012),discriminatebetweenspecies(Horacek,Jakusch,& Krehan,2009),andtoidentifywoodorigin(Gori,Wehrens,LaPorta, &Camin,2015).Strontiumstableisotopeshavealsodemonstrated goodpotentialforprovenancingcedarshiptimber(Richetal.,2016). However,the isotopiccompositionof woodcanbeinfluencedby biologicalandphysicochemicalprocessesnotrelatedtogeographical/climaticfactors(Hajj,Poszwa,Bouchez,&Guerold,2017)and forarchaeologicalwoodfurthercareisneededasprocessessuchas diagenesisinsoilscanmodifywoodcomposition(Hajjetal.,2017). Extendingdendroprovenancingwithtechniquesthatusegeographicallyand climatically boundproxies has thepotentialto strongly improveitsaccuracy. Main conclusions: Integrationofmultiplewood–xylemanatomicalvariablesanalysed withmultivariatetechniquesleadstohigherprecisioninthedendroprovenancingof ring-porousoakspecies. KEYWORDS dendroprovenancing,earlywoodvessels,latewoodwidth,multi-variableapproach,Quercus spp.,region-specificgrowthpatterns,woodanatomy
| 1165 AKHMETZYANOV ET Al. Recently,thepotentialofusingadditionalwood-anatomicalvariablesbesidesringwidthhasalsobeeninvestigated.Esteban,Martin, dePalacios,andFernandez(2012)foundthatvariationofwoodbiometryofapinespecieswashigherandmorepronouncedamong provenancesthanwithinpopulations.Afirstattempttoincludeannualchangesinmaximumlatewooddensityofatreeringshowed thattheprobabilityofsuccessfullydatinghistoricalsamplescanbe greatlyincreased(Wilsonetal.,2017).Therecentdevelopmentof newtechniques(Gärtner&Nievergelt,2010),andimprovementsin imageanalyses(vonArx&Carrer,2014)openednewperspectives forusingmorequantitativewood-anatomicalvariablesfordendroprovenancingstudies(Hanecaetal.,2009). Inthisstudy,weexploredandtestedthepossibilitytouseamultivariableapproachonring-porousoak,aspeciesgroupthatiswidely usedinhistoricconstructions.Inthecontextofthemultidisciplinary projectForSEAdiscovery(htttp://forseadiscovery.eu),weinvestigated whetherusingoakearlywood(EW)vesselpropertiesandlatewood width(LW)canimprovetheprecisionofringwidthbaseddendroprovenancingand,ifso,whatanatomicalvariablesdrivetheimprovement insuchamulti-variableapproach.WehypothesizethatintegratingEW vesselsandLWleadstoamoreprecisedeterminationofthesource areaofagivenoaktimber,becausebothvariablescontaindifferent andcomplementaryclimatesignalswhichdesignatethewoodtospecificareaswherecombinationsoftheseclimatesignalsprevail. 2 | MATERIALS AND METHODS 2.1 | Study sites and species SampleswerecollectedfromnineoakforeststandsintheBasque Country (six sites) and Cantabria (three sites) in Northern Spain (Figure1). Theseforestswereselectedbecausetheyservedaswood-source areasforshipbuildinginthe16thto18thcenturies(MartínezRuiz, 1999).WesampledtreesintheBasquecountry(easternsites)andin Cantabriaprovince(westernsites)(Figure1)withanelevationalgradientintheBasquecountry(from80to1,181ma.s.l.).Theeastern sites(AZK,UZK,OIR,GOR/GOY,OTEandSVA)arelocatedalonga phytogeographicalgradient(Table1),withagradualchangeofvegetationfromEuro–SiberianintheNorthtomoreMediterraneantowardstheSouth,andagradientinelevationrangingfrom80ma.s.l closetothecoastto954ma.s.l.furtherinland(Table1).Thenortheasternsites(AZK,UZKandOIR)havetypicalEuro–SiberianvegetationdominatedbyQuercus robur.AZKandUZKwereintensively usedforcharcoalproductioninthepastbutarenotmanagednowadays.OIRisaprivate-ownedforestpatchinthevillageofOiartzun, withtreesscarcelydistributedonpastureland.Thecentraleastern site(GOR/GOY)islocatedinthetransitionalzonebetweenEuro– SiberianandMediterraneanvegetationandQuercus pyrenaicaisthe dominanttreespecieswithsomeinclusionsofQ. petraea.Thissite washeavilymanaged(coppiced)untilthemid-20thcenturybutwas notmanagedsince.WhereasthesoutheasternsitesOTEandSVA consistofMediterraneanforestsdominatedbyQ. faginea,withOTE onsteephillsandwerenotused/managedinthepast,SVAmostly consistsofoldtreesinflatpastureland(calleddehesainSpanish). Allwesternsites(BRR/BRY,ONQandORT;Figure1)featuretypical transitionalvegetationbetweenAtlanticandMediterraneanphytogeographicalregions,withthedominanceofQ. pyrenaica in ORT and Q. petraeainONQ.Thesesitesshareasimilarmildclimate,withouta pronouncedsummerdrought,andanaltitudinalrangefrom1,037to 1,181ma.s.l.Inouranalysesweseparateeasternfromwesternsites (ascategoricalvariables),andfurtherseparateeasternsitesinnorthern,centralandsoutherncategories(NE,CEandSErespectively). Selectedoakspeciesvaryintheirhabitatandecologyandhence showspecificdistributionalpatterns.Quercus robur and Q. petraea arewidelydistributedinEuropewhereasQ. faginea and Q. pyrenaica arelargelyrestrictedtotheIberianPeninsula(Table2). 2.2 | Acquisition of ring‐width and vessel‐size time series Sampling was performed in May/October 2015 and in October 2016.Betweensevenand25treesweresampledpersite(Table1). Thevariationofsamplesizeperstudysitewascausedbythedifferentnumberofsuitabletreesateachsite.Althoughoakspecies arepronetohybridizeandmanytaxaweredescribedfortheIberian Peninsula(AmaralFranco, 1990), weare confidentof ourspecies identificationaswecarefullyidentifiedeachtreesampledatstands withcoexistingoakspecies.Thelargestandapparentlyoldestoak treeswithnovisibledamagetocrownsorstemswereselectedand twotofour5-mmcoreswereextractedfromeachtree. Sampleswereairdried,gluedontowoodenholders,cutusingaslidingmicrotome(Gärtner&Nievergelt,2010),andsandedwithprogressivelyfinersandpaper(grainsizesfromP180toP1200,FEPAAbrasives). Tyloses (parenchyma outgrowths into vessels) and wood dust were removedusingahigh-pressurewaterblast,afterwhicheachcorewas stainedusingblackink.ChalkwasrubbedonthesamplestofillEWvessels,thusachievingbettercontrastofthevesselsforimageanalyses. Total ring width (TRW), EW and LW were measured using a Velmex measuring table (Velmex TA Unislide; Velmex Inc., Bloomfield,NY)coupledwithastereo-microscope(OlympusSZ60). Crossdatingwasdoneintwosteps:firstvisuallyandthenstatistically.Visualcrossdatingwasdonebyfirstmatchingthering-width patternsofindividualtrees,andthendifferenttreesfromthesame site.Statisticalcrossdatingwasdoneusingcofecha6.02P(Grissino- Mayer,2001).Crossdatingisconsideredasthemostimportantprincipleofdendrochronology(Fritts,1976),asitallowslinkingeachtree ringtotheexactcalendaryearofitsformation. For measurement oftheEW vesselarea,a sequence ofoverlappingimageswitharesolutionof13,780dpiwastakenfromthe woodsurfaceusingadigitalcamera(CanonEOS600D)coupledto the stereo-microscope (Figure2). The images were stitched using PTGui 9.2.0 (New House Internet Services B.V., Rotterdam, The Netherlands).Wemeasuredvesselsonasubsetoftreesandcores showing thehighestinter-correlation (Table1)asvessel measurementsarelabourintensive,andthuscanbeconductedonalimited
1166 | AKHMETZYANOV ET Al. numberof trees only.Moreover, some coreswere brokenor had cracks,whichmadethemlesssuitableforwood-anatomicalanalyses. Earlywoodvesselsweremeasuredsemi-automaticallyinimaGeJ 1.8.0(Schneider,Rasband,&Eliceiri,2012)usingtheVesselJplugin (I.García-González,unpubl.data)andfollowingthestepsdescribed in Souto-Herrero etal. (2018a). This plugin allows for efficient discriminationofthebrightchalk-filledvessel-lumenareafromthe dark-stained background (the remaining wood structures). After scanningandstitching,imageswereconvertedto8-bitblack-and- white and vessels were distinguished from the background using gray-scale thresholds. To ensure that undesired objects were excluded,additionalfilterswereappliedremovingellipticalobjects(less FIGURE 1 Locationofstudysitesanddistributionofsampledtreesateachsite(a–i).Dotsrepresentselectedforeststands(seeTable1 foranexplanationoftheabbreviationsofsitenames),trianglesrepresentselectedtreeindividuals,Quercus robur—inyellow,Quercus pyrenaica—inblue,Quercus petraea—ingreenandQuercus faginea—inred,respectively.Siteswithdoublenamesindicatethecollectionoftwo differentspeciesfromthesameforeststand.Transparentred,transparentgreenandtransparentpurplerefertotheMediterranean,Atlantic andAlpinebiogeographicalregions,respectively(https://www.eea.europa.eu/wileyonlinelibrary.com]).LayerCredits:Source:Esri[Colour figurecanbeviewedatwileyonlinelibrary.com] CANTABRIABASQUE COUNTRY (a) (b) (c) (d) (f)(g)(h)(i) (e)
| 1167 AKHMETZYANOV ET Al. TABLE 1 Overviewonstudy-sitecharacteristicsandsamplingdesign Forest Code Coord. Phyto‐geographic region Samp. region Bed rock Tavg (°C) Prec (mm) Elevation (m a.s.l.) Species TreesTR TreesVESS Azkorte AZK 43°14′N 2°09′W ENE SL 12 1,101 523 Quercus robur 7 6 Oiartzun OIR 43°18′N 1°51′W ENE SS 11.8 1,189 115 Q. robur 10 9 Uzkanga UZK 43°17′N 2°18′W ENE SL 12.8 1,199 80 Q. robur 9 9 Gordoa GOR 42°55′N 2°22′W TCE SS 9.7 878 1,129 Q. petraea 6 6 Gordoa GOY 42°55′N 2°22′W TCE SS 9.7 878 1,129 Q. pyrenaica 16 7 Oteo OTE 42°43′N 2°23′W MSE SL 10.4 653 829 Q. faginea 24 6 SanVicentede Arana SVA 42°43′N 2°21′W MSE SL 10.4 653 954 Q. faginea 25 7 Barrio BRR 43°04′1N 4°39″W TWSS 11.2 963 1,037 Q. petraea 11 10 Barrio BRY 43°04′N 4°39′W TWSS 11.2 963 1,037 Q. pyrenaica 12 9 OnQuemada ONQ 43°05′N 4°42′W TWSS 11.2 963 1,118 Q. petraea 16 10 Orticeo ORT 43°10′N 4°32′W TWSS 11.2 963 1,181 Q. pyrenaica 18 10 Note.Bedrock:SS,sedimentarysilicate;SL,sedimentarylimestone;Coord.,coordinatesofthestudysite;Phytogeographicalregions:E,Euro–Siberian;M,Mediterranean;T,thetransitionfromEtoM; Samp.region—samplingregion:NE,CE,SEandW,northeastern,centraleastern,southeasternandwesternsamplingregions;Tavg,yearlyaveragetemperature(°C)ofthecorrespondingstudysite;Prec, annualprecipitationofthecorrespondingstudysite(mm),(climatedatatakenfromCRUTS4.01(Harrisetal.,2014));TreesTR,numberoftreesusedforringwidth,earlywoodwidthandlatewoodwidth measurements;TreesVESS,numberoftreesusedforearlywoodvesselmeasurements.
1168 | AKHMETZYANOV ET Al. thantwiceaslongaswide)andobjectssmallerthan10,000μm².This particularsizethresholdwastakenfromFontiandGarcía-González (2004)andFonti,Solomonoff, andGarcía-González(2007)toexcludethesmallervesselsformedlaterintheseason,whichmayblur theclimaticsignal(García-González&Fonti,2006),andalsobecause mostofthesignalisrecordedbythelargestvessels,usuallylocated inthefirstrow(García-González,Souto-Herrero,&Campelo,2016; González-González, Rozas, & García-González, 2014). The images withthedetectedvesselswerethencheckedfordetectionerrors andmanualimagecorrectionswereperformedwhereneeded. Toassigneachvesseltothecorrespondingdatedtreering,and toidentifyvesselsbelongingtothefirst-formedrowatthebeginning of the tree ring, we used the software package ”Autovasos” (García-González,unpubl.data),againfollowingSouto-Herreroetal. (2018a). From the EW-vessel measurements, 11 variables related tosize,position,andwater-conductingcapacitywerederivedfrom each tree ring: hydraulically weighted vessel diameter (Dh), mean vesselsarea,maximumvesselsize,numberofvessels(NV),totalvesselsarea,meanareaofthefirstrowvessels,hydraulicdiameter(Dh) ofvesselsnotbelongingtothefirstrow,meanareaofvesselsnot belongingtothefirstrow,maximumsizeofvesselsnotbelonging tothefirstrow,NVnotbelongingtothefirstrow,andtotalareaof vesselsnotbelongingtothefirstrow.Dhwascalculatedaccording totheequation: whereDn (in μm)isthediameter ofthevessel n(Sperry,Nichols, Sullivan,&Eastlack,1994).Dhhasbeenshowntoberelatedtoprevailingenvironmentalconditionsandtomaximizetheclimaticsignal (Pérez-de-Lis,Rossi,Vazquez-Ruiz,Rozas,&García-González,2016; Souto-Herrero,Rozas,&García-González,2017). Before further processing, the individual time series of each variableweredetrendedtoremovetrendsrelatedtoageandforest dynamics,usingacubicsmoothingspline(Cook,Briffa,Shiyatov,& Mazepa,1992)withafrequencyresponseof50%atawavelengthof 32years.Wesubsequentlyusedindividualtimeseriestocompute siteandregionalchronologies.Sitechronologiesincludedonlytrees ofeachsite×speciescombinations,whileregionalchronologiesincludedtreesfromthepredefinedgeographicalregions,namelywest, eastaswellasNE,CEandSE(intheeast)regionsofthestudyarea (seeAppendixS1inSupportingInformationforthedetailedinformationonsiteTRWandDhchronologies). 2.3 | Acquisition of a complementary set of variables Manyofthe14measuredorderivedvariables(3ringwidthand11 vessel-relatedvariables)arecorrelatedandcarryredundantinformation(García-Gonzálezetal.,2016).Toreducethenumberofvariablestoasubsetofdifferent,statisticallymeaningfulvariables,we usedprincipalcomponentanalysis(varimax-rotatedPCA)usingthe sitechronologiescalculatedfromthedetrendedindividualtimeseries.Theselectionwasbasedontheordinationofvariablesaccordingtothefirsttwoprincipalcomponents(PCs).Theresultingsetof complementaryvariables,thatcontaineddifferentpatternsofvariationlikelyrelatedtosite-specificandclimate-relatedsignals,were includedinsubsequentanalyses. 2.4 | Geographic patterns of time–series variability Discriminationbetweenpatternsofannualvariationintheselected variableswasdonebasedonthedetrendedindividualtimeseriesof thetreesfromtheninestudysites.Principalcomponentgradient analyses(PCGA)(Burasetal.,2016)wasperformedinatimeinterval coveredbymoststudytrees(1888–2014).PCGAisbasedonloadingsonthefirsttwoPCsobtainedfromanordinaryPCA,defining agradientamongtimeseriesbasedontheirpolarcoordinates—in particulartheangle—oftheloadings(Burasetal.,2016).PCGAwas usedtotestwhetherourpredefinedbiogeographicalgradientswere reflectedintheloadingswhichrepresentsingletreesonthefirst twoPCs.Toobtainstatisticalconfidenceintheobservedpatterns, wetestedwhetherthepolaranglesrepresentativeoftheloadings differedbetweensourceregionsusingWilcoxonranksumtest.This individualtreeapproach,i.e.usingtimeseriesofindividualtreesinsteadofsitechronologies,waschosensinceoneoftheaimswas to investigate whether individual series of the selected variables D h= ∑ N n=1D 5 n ∑ N n=1 D 4 n TABLE 2 Overviewonstudyspeciescharacteristics,distributionandecology Species Distribution Elevation Bedrock/Soil Quercus robur Europe(Eaton,Caudillo,Oliveira,&deRigo, 2016) 0–1,300ma.s.l.(Eatonetal., 2016) Fertile,moistsoils(Eatonetal., 2016) Quercus petraea Europe(Eatonetal.,2016) 0–2,000ma.s.l(Eatonetal., 2016) Fertile,moistsoils(Eatonetal., 2016) Quercus pyrenaica Morocco,IberianPeninsulaandSouth- WesternFrance(NietoQuintano,Caudillo,& deRigo,2016) 0–2,100ma.s.l.(NietoQuintano etal.,2016) Siliceousmountains,rarely limestoneanddolomites (NietoQuintanoetal.,2016) Quercus faginea(subsp. faginea) IberianPeninsula(BlancoCastroetal.,2005) 600–1,200ma.s.l.(BlancoCastro etal.,2005) Basicsoils(JiménezSancho, DíazFernández,Martín Albertos,&GilSánchez, 1998)
| 1169 AKHMETZYANOV ET Al. expresssite-and/orregion-specificpatterns.Ifso,thiswouldindicatethattheirorigincanbecorrectlyidentified.Forapplicationasa provenancingtool,thisisespeciallyrelevantasoftenonlytimeseries ofasingleorafewtreesareavailablefordendroprovenancing. 2.5 | Climate effect on PCGA gradient ThegradientobtainedfromthePCGAwasfurtherusedtoexplore mechanismspossiblydrivingit.Todothis,wefocusedonclimate correlations of single tree-ring series, instead of calculating conventionalclimate-growthrelationshipsbasedonsitechronologies (Chenetal.,2015;Dorado-Liñán,Akhmetzyanov,&Menzel,2017). WedownloadedgriddedCRUTS4.01temperatureandprecipitationdataovertheperiod1901–2016(Harris,Jones,Osborn,&Lister, 2014)fromtheKNMIclimateexplorer(https://climexp.knmi.nl/start. cgi)andcorrelatedtheseclimatedatawitheachsinglering-width andvesselseries.Toaccountfornon-normallydistributeddata,we usedSpearman'srankcorrelationintheseanalyses.Afterwards,we investigated whether these correlations varied significantly along thedeterminedPCGAgradientusingSpearman'srankcorrelation.In thisanalysisastrongcorrelationwouldindicateadirectionalvariationofclimatecorrelationsalongthePCGAgradient. Moreover,totesttheeffectofthegeographicaloriginonclimate correlations,wealsocomparedtheindividualclimate correlations betweenthegeographicalcategoriesusingWilcoxonranksumtest (for the comparison between east and west), as well as Kruskal– Wallistest(forthenorth,centralandsoutheasternsites).Finally,we computedso-calledregionalchronologies,i.e.theaverageoverall individualtimeseriesrepresentativeofeachgeographicalcategory and performed the corresponding climate correlations. This was donetosupporttheresultsfromthesingletreeclimatecorrelations. 2.6 | Species effect on PCGA gradient In addition to testing the resemblance of geographical categories in the PCGA-loadings, we also explored whether the variation of PCGA-loadingswasrelatedtospecies.Thatis,wetested—inasimilar mannerasforthegeographicalcategories—whetherthePCGA-load- ingsdifferedbetweenthespecies.Thiswasdone,ontheonehand, onalltime–seriesperparameter(i.e.latewoodwidthandvessels). Ontheotherhand,wealsodidthisanalysisonthebasisofdataconsideringonlythosesiteswheretwoofthestudyspecies(Q. petraea and Q. pyrenaica)wereco-existing(BRR/BRYandGOR/GOY). 2.7 | Identifying the area of origin TotestwhethertheoriginofarandomlyselectedtreecanbeidentifiedfromPCGA,weperformeda“leave-one-out”approach.For this,wesuccessivelyselectedeachtreefromthedatasettotest whetheritsorigincouldbecorrectlyidentifiedbymeansofPCGA. Inorder todoso,we excluded the respectivetree fromthedata pool,andperformedaPCGAontheremainingdatasetasexplained insection2.2.ThisanalysiswasdoneseparatelyforLWand vesseltime–series.Afterwards,wecorrelatedtheleft-outtime–series withthescoresofthefirsttwoPCsofthecorrespondingPCGAand transformedthesecorrelationvaluesintoPC-loadings.Fromthese loadingswecomputedEuclideandistancestotheloadingsofallthe othertrees(thosewhichwerekeptinthedataforthePCGA)and computed the categorical mean Euclidean distances referring to geographicalorigins.IfthemeanEuclideandistanceoftheleft-out treetotheothertreesfromthesamegeographicaloriginwaslower incomparisontotreesfromadifferentgeographicalorigin,wedefinedthisasacorrectassignment.Asameasureofsuccessofthis approach,thepercantageofcorrectassignmentswascalculated. AllcalculationsweredoneinR3.4(RDevelopmentCoreTeam, 2017)usingthe”dplR”(Bunn,2010),”psych”(Revelle,2017)and”treeclim”(Zang&Biondi,2015)packages. 3 | RESULTS 3.1 | Selection of complementary variables BasedonthesimilarityofloadingsonthefirsttwoPCsofVarimax rotatedPCA,thenumberofvariablesforfurtheranalyseswasreducedfrom14tothree,i.e.NV,Dh,andLW,sincetheyappeared to carry unique and complementary information for most of the studysiteslikelyrelatedtosite-specificenvironmentalsignals(see AppendixS2). FIGURE 2 Woodsurfaceimagesof: (a)Quercus robur,(b)Quercus petraea,(c) Quercus pyrenaicaand(d)Quercus faginea. EW,LWandVindicateearlywood, latewoodandvessels,respectively[Colour figurecanbeviewedatwileyonlinelibrary. com] (a) (b) (c) (d) EW LW V
1170 | AKHMETZYANOV ET Al. 3.2 | PCGA ThecommonoverlapperiodforthePCGAofLWdatawas127years (1888–2014)andincorporated64individualtreesrepresentingall studysitesandspecies.ThePCGAdemonstratedthatLWtimeseriesweregroupingaccordingtogeographicalcategories,i.e.there wasaclearseparationbetweentheloadingsrepresentingtreesof easternandwesternorigin(Figure3)(p<0.01,Wilcoxontest).Trees showedasimilargroupingaccordingtotheloadingsasderivedfrom NVtimeseries;however,thereweremoremismatchescomparedto LW(seeAppendixS3). Dhtimeseriesrevealednoclearseparationofloadings,neither across sites nor between the regions (see Appendix S4). Nevertheless,theapplicationofPCGAonlyconsideringtreesfrom theeasternsites,i.e.northeast—centraleast—southeast(NE—CE— SE),allowedgroupingthetreesaccordingtothegeographicalcategoriesbasedonDh(Figure4).Inthisanalysis,thecommonoverlap periodusedforthePCGAwas102years(1913–2014)involving35 trees.AKruskal–WallistestindicatedsignificantdifferencesofloadingsbetweentreesfromNE,CEandSE(p<0.01).Atthesametime, itwasnotpossibletodistinguishtheCEtreesfromtheothertwo sitesonthebasisofLWandNV(seeAppendixS5). 3.3 | Climate effect on PCGA gradient CorrelationanalysesbetweenthegradientderivedfromthePCGA ofLWtimeseriesandclimatecorrelationscoresindicateastrong effectforaverageJune–JulytemperatureonthePCGAgradient(see AppendixS6).TreesfromtheeasternsiteswerelimitedinLWformationbysummertemperatures,whereasthosefromthewestern sitesdidnotshowanysignificantrelationshipbetweenLWandearly summertemperature(Figures5a,b). The strength of the climate correlations also varied along the PCGA gradientas derivedfromDh(seeAppendixS7).Vesselsize (Dh)oftreesinNE,locatedatlowerelevationandclosertothesea wasnegativelycorrelatedwithFebruarytemperature,whereasfor CE and SE correlations with February temperatures were weaker (Figure6a,b).Atthesametime,DhoftreesfromSE,i.e.thehigher andmoresouthernsitesweremostlypositivelycorrelatedwithaverageMarch–Maytemperature,whilethisrelationshipwasnotfound fortheotherregions(Figure6c,d). Principal component gradient analyses indicated no speciesspecific effects on the observed gradient neither for LW nor for Dhwheninvolvingalltreeandoakspecies(seeAppendixS8).This observationwassupportedwhenapplyingPCGAonlytodatafrom sites(GOR/GOYandBRR/BRY)wheretwospeciesco-existed(see AppendixS9). 3.4 | “Leave‐one‐out” analyses The“leave-one-out”analysisindicatedthatmosttreeswerecorrectly assignedtotheirspecificsiteoforiginbasedonLWandDh(Table3). Thatis,forLW89%and100%ofeasternandwesterntreeswere assignedcorrectlyrespectively.However,whenanalysingLWofonly theeasternsites,nosampleswereassignedtothecentraleastsite.At thesametime,Dhchronologiessuccessfullydetectedsamplescomingfromnortheast(81%),centraleast(58%)orsoutheast(91%). FIGURE 3 Loadingsplotfortheprincipalcomponentgradient analysesbasedonlatewoodwidth.Eacharrowrepresentsthe loadingsonthefirsttwoprincipalcomponentsfromonetree.Value labelsonthex-andy-axes(0.19and0.07)representtheamountof explainedvariancebytherespectivecomponent[Colourfigurecan beviewedatwileyonlinelibrary.com] −0.2 −0.1 0.00.1 0.2 −0.2 −0.1 0.00.1 0.2 0.19 0.07 Eastern trees Western trees FIGURE 4 Loadingsplotfortheprincipalcomponentgradient analysesbasedonhydraulicdiameterbutonlyconsideringthe easternsites.Eacharrowrepresentstheloadingsonthefirsttwo principalcomponentsfromonetree.Valuelabelsonthex-and y-axes(0.19and0.09)representtheamountofexplainedvariance bytherespectivecomponent[Colourfigurecanbeviewedat wileyonlinelibrary.com] −0.3 −0.2 −0.1 0.00.1 0.20.3 −0.3 −0.1 0.00.1 0.20.3 0.19 0.09 North−Eastern trees Central−Eastern trees South−Eastern trees
| 1171 AKHMETZYANOV ET Al. 4 | DISCUSSION 4.1 | Added value of including complementary treering variables into dendroprovenancing of oak Our results indicate that a higher precision of provenancing can be achievedforring-porousoakwoodwhencombiningwood-anatomi- cal variables such as LW and size of EW-vessels, rather than using TRW(orLW)alone,asisusuallythecaseindendroprovenancing(e.g. Domínguez-Delmásetal.,2014;Haneca&Debonne,2012).Wefound thatannualvariationpatternsinLWandvesselDhofthestudyspecies wererelatedtothebiogeographicalorigin,especiallyforDh.Bothvariablesareprobablyreflectingdifferentclimaticsignals(winterandspring temperatureforDh,summertemperature forLW [Figures4and5]), andthereforeconditionsduringdistinctperiodsoftheyearthatlargely determinethegrowthofoakinaspecificarea.LWallowedustoassignoakstoBasque(eastern)andCantabrian(western)sites(Figure3), whereasDhdiscriminatedoakpopulationsacrossalatitudeandelevationgradient(NE-CE-SEinFigure4),whichalsorepresentedincreasingcontinentality.Thus,combiningLWandDhallowsforresembling region-specificclimateconditionsoverdifferentpartsofthegrowing season,whichhencecanbeusedtopinpointtheoriginofwood. Theproposedcombinedapproachisalsointerestingforotherareas, such as the maritime part of Central Europe, where classical dendroprovenancing of oak fails (Bridge, 2000), probably due to its complex signalinTRW(Bridge,2012).Moreover,TRWandLWofQ. pyrenaica innorthwesternSpainwerefoundtobeinfluencedbyforestdynamics, whichishardlythecaseforEWvessels(Souto-Herrero,Rozas,&García- González,2018b).Thisisespeciallyinterestingforarchaeologicalstudies, aswoodusedfor(ship-)constructionoftencomesfrommanagedforests. Wood-anatomical analyses can also be applied to waterlogged timber,whereasthere still isa knowledge gap about the change of woodcompositionunderwaterloggedorburiedconditions.Isotopic analysesdonotnecessarilyimprovewoodprovenancingprecisionof livingtreesaswell(Vlametal.,2018),thereforetheyneedtobefurther testedfordifferentspeciesandregions.However,isotopesmightcontainsite-specificsignalscomplementarytowood-anatomicalvariables (e.g. maximum summer temperature and precipitation in δ13C and δ18Ointemperateregions[Treydteetal.,2007], strontium—reflects thesoilsignature[Hajjetal.,2017]).Acombinedmulti-proxydendroprovenancingmethodthusisthemostpromisingapproachforfuture studies.Andintegrationwithothertechniquesrecentlytestedfordendroprovenancing,e.g.Fouriertransforminfraredspectroscopy(Traore, Kaal,&Cortizas,2018),seemstobroadenthispotential. Ourresultsalsodemonstratethatatree-centredapproach,i.e. using single trees as a main source of information to understand variabilityingrowthresponses,asproposedforecologicalstudies (Sass-Klaassen etal., 2016) allowspin-pointing tree originat both inter-andintra-regionalscales.Suchapproachalsoovercomesthe problemthatoftenonlyfewtreesareavailablefordatingandprovenancingpiecesofart(paintings,musicalinstruments,etc.). 4.2 | Spatial variation of climatic signal in different tree‐ring variables 4.2.1 | Climatic signal in latewood width Analyses of spatial variation of the climatic signal revealed that the main factor discriminating (based on LW time series) tree FIGURE 5 Relationshipsbetweenclimatecorrelationsandtheprincipalcomponentgradientanalyses(PCGA)-rankderivedfrom latewoodtimeseries.(a)Correlationsbetweenlatewoodwidthtime–seriesandaverageJune–Julytemperaturesforindividualtreesplotted againstindividualtreePCGA-ranks.PCGA-rankreferstotherankofatree-ringseriesinthegradientdefinedbyPCGA.(b)Correlations ofindividualtreelatewoodwidthserieswithaverageJune–Julytemperature,separatedaccordingtogeographicorigin.Blackdotsin(b) representthecorrelationcoefficientofregionalchronologieswithaverageJune–Julytemperature.Blueandredcolorsrefertoeasternand westerntreesrespectively.Notethatthey-axesspantheintervalfrom−0.5to0.5[Colourfigurecanbeviewedatwileyonlinelibrary.com] 0102030405060 −0.4 −0.2 0.00.2 0.4 PCGA−rank Jun−Jul temperature correlation score (a) West East −0.4 −0.2 0.00.2 0.4 (b)