Geography, climate and shifts in host plants distribution explain the genomic variation in the cactus moth

Geog aphy, clima e and shi s in hos plan s dis ibu ion explain he genomic a ia ion in he cac us mo h Daniel Po eda-Ma ı ´nez 1,2 *,Vı ´c o Nogue ales 3 , S ephen D. High 4 , Guille mo Loga zo 1 , B en C. Eme son 3 , Lau a Va one 1 † and Es eban Hasson 2 † 1 Fundacio ´n pa a el Es udio de Especies In asi as (FuEDEI), Hu lingham, A gen ina, 2 Ins i u o de Ecologı ´a Gene ´ ica y E olucio ´n de Buenos Ai es (IEGEBA-CONICET), Facul ad de Ciencias Exac as y Na u ales, Uni e sidad de Buenos Ai es, Buenos Ai es, A gen ina, 3 Island Ecology and E olu ion Resea ch G oup, Ins i u o de P oduc os Na u ales y Ag obiologı ´a (IPNA-CSIC), San C is o ´bal de La Laguna, Cana y Islands, Spain, 4 The U.S. Depa men o Ag icul u e –Ag icul u al Resea ch Se ice (USDA-ARS), Insec Beha io and Biocon ol Resea ch Uni (IBBRU), Tallahassee, FL, Uni ed S a es Landscape he e ogenei y and he hos plan use a e ac o s sugges ed o play de e minan oles in shaping he e olu iona y his o y o he bi o ous insec s. Howe e , he ole o he econfigu a ion o hos plan s dis ibu ions linked o Qua e na y clima e oscilla ions as d i e s o con empo a y popula ion gene ic s uc u e is s ill poo ly unde s ood. He e, we o mally examine he ela i e con ibu ion o such ac o s on in aspecific di e sifica ion using he Sou h Ame ican cac us mo h, Cac oblas is cac o um, an he bi o e insec specialized in he use o cac i as hos plan s. We assessed genomic a ia ion using genome- wide SNPs and mi ochond ial da a in popula ions sampled ac oss a b oad geog aphical g adien whe e mo hs eed on di e en cac us species. We in eg a ed demog aphic simula ions and ecological niche modeling in o a landscape genomics amewo k, o es al e na i e hypo heses o pas and cu en popula ion connec i i y o bo h C. cac o um and i s hos plan s. Regions exhibi ing highe genomic di e si y we e e alua ed o cong uence wi h a eas whe e sui able clima ic condi ions emained s able h ough ime. Ou esul s e ealed ha pas spa ial configu a ion o sui able habi a condi ions and shi s o hos plan s dis ibu ions a e he ac o s ha be e explain he in aspecific di e sifica ion. Genomic da a also suppo ed he hypo hesis ha a eas o long- e m habi a s abili y se ed as e ugia o C. cac o um, enabling he main enance o high le els o gene ic di e si y o e ime. O e all, ou s udy highligh s he impo ance o in eg a ing in e -specific in e ac ions and hei spa io- empo al dynamics o be e unde s and he ela i e impo ance o abio ic and bio ic ac o s d i ing he di e sifica ion p ocesses in he bi o ous insec s wi h b oad geog aphical and es ic ed hos anges. KEYWORDS cac us pes , cac oblas is cac o um, ddRAD, landscape genomics, paleoclima e, opun ia F on ie s in Ecology and E olu ion on ie sin.o g01 OPEN ACCESS EDITED BY Zi ong Li, Commonweal h Scien ific and Indus ial Resea ch O ganisa ion (CSIRO), Aus alia REVIEWED BY Cui Wang, Uni e si y o Helsinki, Finland Jo dan Golubo , Uni e sidad Au o ´noma Me opoli ana, Mexico *CORRESPONDENCE Daniel Po eda-Ma ı ´nez [email p o ec ed] † These au ho s sha e las au ho ship RECEIVED 18 July 2023 ACCEPTED 29 Augus 2023 PUBLISHED 19 Sep embe 2023 CITATION Po eda-Ma ı ´nez D, Nogue ales V, High SD, Loga zo G, Eme son BC, Va one L and Hasson E (2023) Geog aphy, clima e and shi s in hos plan s dis ibu ion explain he genomic a ia ion in he cac us mo h. F on . Ecol. E ol. 11:1260857. doi: 10.3389/ e o.2023.1260857 COPYRIGHT ©2023Po eda-Ma ı ´nez, Nogue ales, High , Loga zo, Eme son, Va one and Hasson. This is an open-access a icle dis ibu ed unde he e ms o he C ea i e Commons A ibu ion License (CC BY). The use, dis ibu ion o ep oduc ion in o he o ums is pe mi ed, p o ided he o iginal au ho (s) and he copy igh owne (s) a e c edi ed and ha he o iginal publica ion in his jou nal is ci ed, in acco dance wi h accep ed academic p ac ice. No use, dis ibu ion o ep oduc ion is pe mi ed which does no comply wi h hese e ms. TYPE O iginal Resea ch PUBLISHED 19 Sep embe 2023 DOI 10.3389/ e o.2023.1260857 1 In oduc ion The spa ial dis ibu ion o genomic a ia ion in na u al popula ions is influenced by a a ie y o p ocesses, including he dispe sal abili y o indi iduals, gene flow, bio ic in e ac ions, and landscape ea u es. Un eiling he ela i e oles o hese p ocesses is necessa y o unde s and how o ganisms in e ac wi h he en i onmen and how he la e shapes hei genomic makeup. In his ma e , he bi o ous insec s dis ibu ed along en i onmen al g adien s eeding on a se o hos plan s o e he oppo uni y o disen angling he ecological and geog aphical p ocesses d i ing popula ion di e en ia ion and gene flow dynamics in he e ogeneous landscapes (Bo e e al., 2012;Laukkanen e al., 2014). Among he ecological ac o s, he use o al e na i e hos plan s is hough o shape pa e ns o gene flow among popula ions in he bi o ous insec s (Funk e al., 2006;Fo bes e al., 2017). Mo eo e , changes in hos plan use oge he wi h adap a ion o di e en en i onmen al condi ions can lead o ecological specializa ion and di e sifica ion (isola ion by en i onmen , IBE; Wang and B adbu d, 2014). Geog aphy (isola ion by dis ance, IBD; Sla kin, 1993) and he spa ial dis ibu ion o sui able habi a s (isola ion by esis ance, IBR; McRae, 2006;McRae and Beie , 2007) may also be de e minan s o popula ion gene ic di e en ia ion ac oss he landscape (Pe e man e al., 2014;Vidal e al., 2019). Insec s wi h na ow die b ead hs and en i onmen al equi emen s a e expec ed o exhibi ma ked gene ic s uc u e as a consequence o gene flow dis up ion be ween popula ions exploi ing al e na i e hos s (Vidal and Mu phy, 2018) and/o inhabi ing a agmen ed landscape o sui able habi a s (O ego e al., 2021). Howe e , in widely dis ibu ed insec s eeding on a na ow a ay o hos plan s ac oss he landscape, con empo a y pa e ns o genomic a ia ion can also be influenced by a ange o ac o s, including dispe sal abili y, pa e ns o hos plan use and/o componen s o landscape he e ogenei y. Mo eo e , less s udied spa io-ecological ac o s, such as dis ibu ional shi s o hos plan s d i en by clima e oscilla ions h oughou he Qua e na y pe iod, may ha e also p omo ed cycles o geog aphical isola ion and seconda y con ac , shaping con empo a y pa e ns o gene ic a ia ion in he bi o ous insec s (Nogue ales e al., 2018). Ye , he ela i e con ibu ion o hese ecological ac o s and spa io- empo al landscape dynamics in d i ing popula ions di e en ia ion in axa wi h b oad clima ic niches and na ow eeding equi emen s a e s ill open ques ions in e olu iona y biology. To unde s and he ole o such spa ial-ecological p ocesses in species di e sifica ion and shaping gene flow in he e ogeneous landscapes, an in eg a i e app oach is necessa y, aking in o accoun his o ical and con empo a y en i onmen al da a, combined wi h a sui able insec model wi h na ow ophic equi emen s and dis ibu ed along en i onmen al g adien s. He e, we in es iga e pa e ns o genomic a ia ion in he cac us mo h, Cac oblas is cac o um (Be g) (Lepidop e a: Py alidae), a sou he n Sou h Ame ican he bi o e insec , ha inhabi s he e ogeneous landscapes, aces a a ie y o clima e condi ions, and eeds on se e al species o he genus Opun ia (Cac aceae, Opun ioideae). This mo h can be ound ac oss wide la i udinal and longi udinal g adien s, om he a id and semi-a id lands o no hwes e n and no heas e n Chaco, o he g asslands o he Pampa, and no h o he Pa agonian s eppe (McFadyen, 1985; Mo one, 2014;Va one e al., 2014). Along such b oad dis ibu ion, C. cac o um eeds on na i e species, such as O. anacan ha Speg., O. bonae ensis Speg., O. ela a Link & O o ex Salm-Dyck, O. megapo amica A echa ., O. penicillige a Speg., O. quimilo K. Schum., and O. iopla ense Fon ., and also he cul i a ed O. ficus-indica (L.) Mille , an exo ic species in he egion (Va one e al., 2012;Va one e al., 2014). Soon a e i s in oduc ion he p ickly pea O. ficus-indica was quickly adop ed as a c op in sou he n Sou h Ame ica (E in, 2012). The a ailabili y o his new hos plan ep esen s a ela i ely ecen expansion o hos use o C. cac o um. Cu en ly his mo h is a se ious h ea o he egional indus y o p ickly pea ui p oduc ion (Ochoa e al., 2007). The in oduc ion o he new hos led o he sugges ion ha he sp ead o O. ficus-indica p omo ed he dispe sion and geog aphical ange expansion o C. cac o um (Va one e al., 2014). P e ious s udies e ealed ha popula ions o C. cac o um a e gene ically di e en ia ed in a leas ou haplog oups co ela ed wi h geog aphy and a iable le els o gene ic di e si y among na i e and non-na i e mo h popula ions (Ma sico e al., 2011), and ha pa e ns o colo a ion o ma u e la ae and hos plan use a y ac oss egions (B ooks e al., 2012). Howe e , only limi ed conclusions can be d awn om hese s udies since (i) only a single mi ochond ial gene ic ma ke was used; (ii) sampling did no co e he en i e species ange; and (iii) he ela i e abundance o hos species was no conside ed (McFadyen, 1985;Ma sico e al., 2011;B ooks e al., 2012). Mo e ecen ly, an ex ensi e su ey e alua ing hos plan p e e ences ac oss he b oad ange o C. cac o um concluded ha pa e ns o hos use eflec a ia ion in hos a ailabili y a he han emale p e e ence (Va one e al., 2014). In his s udy, we use high- h oughpu nuclea sequence da a and Sange sequence da a o mi ochond ial DNA o in es iga e he p ocesses ha shaped pa e ns o di e sifica ion in he cac us mo h, C. cac o um. Ou sampling s a egy ex ended o e a b oad geog aphic ange aking ca e o including he di e si y o na i e and non-na i e hos s. We combine in e ences om phylogenomics and demog aphic modeling o unco e he empo and mode o di e sifica ion, pa e ns o gene flow and popula ion gene ic s uc u e. By in eg a ing ecological niche modeling in o a landscape genomic amewo k, we es a comp ehensi e sui e o compe ing scena ios o popula ion isola ion. Specifically we e alua e he hypo heses ha con empo a y spa ial pa e ns o popula ion gene ic di e en ia ion in C. cac o um a e p edominan ly a consequence o he spa ial econfigu a ion o sui able a eas and shi s in he dis ibu ion o hos plan s caused by Qua e na y clima ic oscilla ions whils he cu en pa e n o hos plan use and opoclima ic ac o s a e less influen ial in i s di e sifica ion. Since we es ima ed his o ical and con empo a y habi a sui abili y, and long- e m sui able en i onmen s end o suppo la ge e ec i e popula ion sizes (Ca na al e al., 2009; Soley-Gua dia e al., 2019), we also add ess whe he egions whe e sui able habi a condi ions emained s able since he las glacial maximum (LGM, ca. 21 ka) o he p esen exhibi g ea e gene ic di e si y han popula ions h i ing in less s able a eas. Finally, we e alua ed whe he he expansion o O. ficus-indica acili a ed Po eda-Ma ı ´nez e al. 10.3389/ e o.2023.1260857 F on ie s in Ecology and E olu ion on ie sin.o g02 con empo a y gene flow among geog aphically dis an popula ions. In his ma e , ou expec a ion is o find oo p in s o gene ic admix u e among dis an egions, especially among a eas highly impac ed by p ickly pea ag oindus y. 2 Me hods 2.1 Popula ion sampling and DNA ex ac ion Be ween 2018 and 2019, we collec ed 140 indi iduals o C. cac o um om 28 sampling si es spanning a b oad ange o i s na i e dis ibu ion in sou he n Sou h Ame ica (Figu e 1;Table S1). Indi iduals we e collec ed di ec ly om hos plan s, including se en na i e axa, O. anacan ha, O. bonae ensis, O. ela a, O. megapo amica, O. penicillige a, O. iopla ense, and O. quimilo, and he exo ic O. ficus-indica. To a oid he inclusion o siblings, only a single indi idual pe sampled plan and eggs ick was used o gene ic analysis. Be o e DNA ex ac ion, indi idual la ae we e dissec ed o emo e any ood esidue, and ~30 mg o clean issue was used o DNA ex ac ion wi h Qiagen DNeasy Blood & Tissue Ki ollowing manu ac u e ’s ins uc ions (Qiagen, Inc.). 2.2 Genomic lib a ies p epa a ion, fil e ing and assembling Indi iduals we e p ocessed in o wo genomic lib a ies ollowing he double diges es ic ion-si e Associa ed DNA (ddRAD) p ocedu e desc ibed in Pe e son e al. (2012). Genomic DNA was agmen ed using NspI and MboI es ic ion enzymes, and hen pu ified and liga ed o ba coded adap e s. Lib a ies we e sequenced in pai ed-end 125 bp mode on a HiSeq2500 Illumina ins umen . Demul iplexed aw eads we e checked, compiled and summa ized on Mul iQC 1.10.1 (Ewels e al., 2016). Reads we e assembled wi h IPYRAD .0.9.59 (Ea on and O e cas , 2020). A e e ence based assembly me hod was implemen ed o mapping he fil e ed eads agains he C. cac o um e e ence genome. DNA ex ac ed om an adul male (clean DNA, ee o ood con aminan s can be ob ained om males, as hey do no eed) was used o lib a y p epa a ion and sequenced in an Illumina No aSeq 6000 pla o m. The de no o assembly was accomplished wi h CLC Genomics Wo kbench 7.1 (see Po eda-Ma ınez e al., 2022 o de ails). Thus, he e e ence genome se ed also as a fil e o exogenous eads de i ed om la al ood. In subsequen IPYRAD assembly s eps, we allowed 20% o SNPs pe RAD locus, and sha ed he e ozygous si es occu ing ac oss a maximum o 25% o samples. A minimum o 90% was se o he minimum numbe o samples sco ed pe locus. We used c ools 1.15 (Danecek e al., 2011) o emo eSNPswi h minimum allele equency lowe han 3%, missing da a pe si e ac oss indi iduals exceeding 25%, and o keep SNPs wi h ead dep hs be ween 6× and 100×. Indi iduals wi h mo e han 20% missing da a we e excluded om u he analyses. To p une SNPs in linkage disequilib ium, PLINK 1.9 (Pu cell e al., 2007) was used wi h a window size o 50 bp, window shi o 5 and VIF h eshold o 2. SNPs unde selec ion we e iden ified using Bayescan .2.1 (Foll and Gaggio i, 2008) and excluded om he da ase . A agmen o he mi ochond ial gene encoding he Cy och ome oxidase subuni I AB C FIGURE 1 Geog aphic loca ion and hos species o he Cac oblas is cac o um popula ions sampled in his s udy (A). Regions ep esen biogeog aphic p o inces acco ding o Mo one e al. (2014). Panels (B, C) ep esen he maps o habi a sui abili y o C. cac o um du ing he p esen (B) and du ing he Las Glacial Maximum (LGM, ca. 21ky) (C) as in e ed by he mos -suppo ed ENM model in MAXENT. Po eda-Ma ı ´nez e al. 10.3389/ e o.2023.1260857 F on ie s in Ecology and E olu ion on ie sin.o g03 (COI) (m DNA) was amplified using p ime s C1-J-2183 and Pa II. Amplifica ions, PCR-p oduc pu ifica ion and sequencing we e pe o med as de ailed in Po eda-Ma ınez e al. (2022). 2.3 Assessing genomic popula ion s uc u e Popula ion gene ic s uc u e was in e ed using genome-wide SNP da a (nDNA he ea e ) and he spa se non-nega i e ma ix ac o iza ion app oach (sNMF) (F icho e al., 2014), as implemen ed in R using he LEA package (F icho and F ancois, 2015). The numbe o gene ic clus e s ha bes desc ibed ou da a was assessed wi h 100 epe i ions o each possible K alue ( om K = 1 o K = 10) and using he c oss-en opy c i e ion. In addi ion o sNMF, gene ic s uc u e was app oxima ed by using he majo axes o genomic a ia ion ob ained om a P incipal Componen Analysis (PCA), as implemen ed in he hie s a (Goude , 2005)R package. We used m DNA da a o cons uc a haplo ype ne wo k using he Median Joining algo i hm (Bandel e al., 1999)as implemen ed in PopA .1.7.1 (Leigh and B yan , 2015). 2.4 Phylogene ic analyses Phylogene ic ela ionships we e econs uc ed among he main gene ic g oups as in e ed in sNMF using nDNA da a and he coalescen -based me hod implemen ed in SNAPP (B yan e al., 2012). Due o compu a ional bu den, SNAPP analyses we e un including only fi e indi iduals pe gene ic clus e , selec ing hose wi h an ances y coe ficien highe han 0.85, acco ding o sNMF, and indi iduals ha had he lowes p opo ion o missing da a. One SNP pe locus was andomly selec ed esul ing in a new ma ix o 2,086 SNPs sha ed ac oss ips. This allows o ensu e independen biallelic ma ke s as well as ai ly spaced SNPs as equi ed in SNAPP (B yan e al., 2012). The de aul model pa ame e s we e used in SNAPP o U and V equal o one, and he analysis was un o 5,000,000 MCMC gene a ions, sampling e e y 1000 gene a ions in Beas .2.5.2 (Bouckae e al., 2014). The comple e se o ees was isualized in Densi ee .2.2.5 (Bouckae , 2010), emo ing he fi s 10% o he ees as bu n-in. A maximum c edibili y ee was gene a ed using T eeAnno a o .1.7.5 (D ummond e al., 2012). 2.5 Demog aphic his o y A coalescen -based simula ion app oach was implemen ed in Fas simcoal2 (Exco fie e al., 2013) based on he si e equency spec um (SFS) o nDNA da a o u he in es iga e he demog aphic his o y o he h ee main lineages iden ified in clus e ing analysis and species ee in e ence (Cen e , Eas and Sou h, see esul s sec ion). Ini ially, o u he e alua e he consis ency o phylogene ic ela ionships among C. cac o um popula ions in e ed in SNAPP, models ep esen ing he h ee possible opologies (A, B and C) we e conside ed in each one o he subsequen models o gene flow (Figu e S1). Fo each opology, di e gence imes, con empo a y and ancien e ec i e popula ion sizes we e es ima ed, as well as al e na i e scena ios o gene flow. Models 1–3 conside ed no gene flow among lineages; models 4–6 con empla ed scena ios o symme ic gene flow among lineages, while models 7–9 conside ed scena ios o asymme ic gene flow among lineages. Fo hese analyses, 20 indi iduals we e selec ed wi h he highes ances y coe ficien (>0.85) acco ding o sNMF (e.g., O ego e al., 2021). The olded join SFS was calcula ed conside ing a single SNP pe locus using he Py hon sc ip w i en by Isaac O e cas and a ailable a Gi Hub (h ps://gi hub.com/ isaaco e cas /easySFS). To emo e missing da a and minimize e o s wi h allele equency es ima es, each gene ic clus e was downsampled o 16 indi iduals yielding a o al o 2,565 a iable si es. Assuming ha in a iable si es we e no conside ed in he SFS calcula ion, we used he “ emo eZe oSFS”op ion in Fas simcoal2 and fixed he e ec i e popula ion size o one o he demes (Sou h lineage, NE SOUTH = 824,626) o enable he es ima ion o o he pa ame e s wi h Fas simcoal2. NE SOUTH was calcula ed acco ding o he equa ion NE = p/4m(Lynch and Cone y, 2003). Nucleo ide di e si y (p= 0.004) was es ima ed using a ian and in a ian si es wi h DNAsp 6 (Rozas e al., 2017). The mu a ion a e (m) was conside ed o be 2.9 × 10 −9 subs i u ions pe si e pe gene a ion, as p e iously es ima ed o he bu e flyHeliconius melpomene (Keigh ley e al., 2015) (Lepidop e a: Nymphalidae). Each model was un wi h 50 eplica es, conside ing 100,000–250,000 simula ions o he calcula ion o he composi e likelihood, 10–40 expec a ion–condi ional maximiza ion (ECM) cycles, and a s opping c i e ion o 0.001. Once maximum likelihood was es ima ed pe un, he bes demog aphic model was selec ed acco ding o Akaike’s in o ma ion c i e ion (AIC). AIC alues we e escaled in e ms o AIC di e ences (Di) acco ding o he o mula: Di = AICi −AICmin. A model wi h a DAIC alue o 0 and he highes AIC weigh (wi) se ed as he bes model. A pa ame ic boo s apping app oach was used o cons uc 95% confidence in e als o he es ima ed pa ame e s unning 100 boo s ap eplica es using ini ialized alues om he bes model (Exco fie e al., 2013). 2.6 En i onmen al niche modeling An en i onmen al niche modeling app oach was implemen ed o p edic con empo a y and his o ical geog aphic dis ibu ions o sui able a eas o he cac us mo h in i s na i e ange. En i onmen al niche models (ENMs) we e used o de e mine he impac o his o ical and con empo a y componen s o landscape he e ogenei y on C. cac o um di e sifica ion as is de ailed below in he landscape genomic analysis sec ion. ENMs we e buil using MaxEn .3.4.1 (Phillips e al., 2006) and he bioclima ic da a a ailable in CHELSA .1.2 da abase (h ps://chelsa-clima e.o g/ bioclim/). Model pa ame e s in MaxEn we e selec ed and op imized using he kuenm (Cobos e al., 2019) R package. Fo bo h empo al bioclima ic condi ions, p esen -day and Las Glacial Maximum (LGM, ca. 21 kya), we e alua ed 15 en i onmen al a iables e ie ed a 30 a c-sec (~1 km) o esolu ion. Fou bioclima ic a iables (Bio 8, 9, 18 and 19) we e disca ded o Po eda-Ma ı ´nez e al. 10.3389/ e o.2023.1260857 F on ie s in Ecology and E olu ion on ie sin.o g04 ha ing a ificial b eaks (Oli ei a e al., 2020). Highly co ela ed a iables (R > 0.9) acco ding o a iance infla ion ac o c i e ion we e excluded o downs eam analyses esul ing in a final da ase o six bioclima ic a iables (Bio 2, 3, 5, 13, 14 and 15). Sui abili y maps du ing he LGM we e ob ained by p ojec ing he p esen -day ENM on o LGM bioclima ic condi ions de i ed om he Communi y Clima e Sys em Model (CCSM4; B aconno e al., 2007) esul ing in wo sui abili y maps based on bioclima ic a iables da a (Se 1, Clima ic CURR and Clima ic LGM om he ea e ). Addi ionally, we cons uc ed ENMs o fi e o he se en na i e hos species o C. cac o um conside ing he a o emen ioned app oach o a iable selec ion and model building. ENM o O. penicillige a and O. bonae ensis could no be implemen ed due o limi ed occu ence da a. The esul ing ENMs ob ained o cu en and LGM condi ions om each hos species we e used as a iables o build an addi ional ENM o C. cac o um (Se 2, Hos CURR and Hos LGM om he ea e ). A hi d ENM was cons uc ed o C. cac o um conside ing he inpu o bo h he bioclima ic a iables (Se 1) and he hos plan ENMs (Se 2) (Se 3; Clima ic-Hos CURR and Clima ic-Hos LGM om he ea e ). Occu ence da a o he hos species and C. cac o um we e ob ained om field su eys made om 2007 o 2019 and complemen ed, a e de ailed cu a ion, wi h dis ibu ion in o ma ion a ailable a GBIF (www.gbi .o g). Redundan occu ences (e.g. poin s occu ing wi hin 1,500 km 2 )we e excluded using spThin (Aiello-Lammens e al., 2015) R package. A e hinning occu ence da a, 81 eco ds o C. cac o um loca ions and hei associa ed occu ance o 205 hos species emained [O. quimilo (50), O. megapo amica (47), O. iopla ense (39), O. ela a (38), and O. anacan ha (31)] and we e used o conduc ENMs. Model pe o mance o each scena io was e alua ed independen ly based on s a is ical significance (Pa ial ROC), omission a es (OR), and he Akaike in o ma ion c i e ion co ec ed o small sample sizes (AICc) using he kuenm (Cobos e al., 2019) R package. 2.7 Landscape genomic analyses A landscape genomic app oach was implemen ed o s udy po en ial ac o s ha could explain pa e ns o gene ic di e en ia ion wi hin C. cac o um.Asameasu eo gene ic di e en ia ion, pai wise F ST es ima es we e de i ed om genome- wide SNP da a using he Wei and Cocke ham (1984) me hod wi h he S AMPP (Pemble on e al., 2013) R package and 9,999 boo s ap eplica es. We e alua ed se e al plausible scena ios o popula ion connec i i y based on his o ical and con empo a y spa ial and ecological da a. Isola ion by esis ance scena ios (IBR) we e es ed by calcula ing esis ance su aces based on he sui abili y maps ob ained om ENMs o p esen -day and LGM condi ions conside ing di e en subse s o ac o s: (i) only clima ic a iables (Clima ic CURR and Clima ic LGM ); (ii) clima e-based habi a sui abili y maps conside ing he hos species dis ibu ion (Hos CURR and Hos LGM ); and (iii) combining clima ic a iables and clima e-based habi a sui abili y maps o hos species dis ibu ion (Clima ic-Hos CURR and Clima ic-Hos LGM ). Resis ance dis ances o all pai s o popula ions we e calcula ed using an eigh -neighbo cell connec ion scheme in Ci cui scape .5 (Hall e al., 2021) h ough Julia .1.5.2 (h ps://julialang.o g/). We also calcula ed esis ance dis ances based on a “fla landscape” whe e all cells ha e an equal esis ance alue (=1) ep esen ing a null model o isola ion by esis ance (IBR NULL ). An isola ion by dis ance scena io was also e alua ed using weigh ed opog aphic dis ances (IBD WTD ), which inco po a e an addi ional o e land dis ance co e ed by an o ganism due o changes in ele a ion imposed by he opog aphy. The IBD WTD dis ance ma ix was calcula ed on a digi al ele a ion model (DEM) a ~1 km o esolu ion e ie ed om Wo ldClim .2.1 (Fick and Hijmans, 2017) da ase using he TopoWeigh edDis unc ion implemen ed in he opoDis ance (Fick and Hijmans, 2017;Wang, 2020)R package. We assumed a linea unc ion o weigh aspec changes (hFunc ion pa ame e ) and an exponen ial unc ion o weigh he slope ( Func ion pa ame e ) (e.g. Nogue ales e al., 2021). Clima ic dissimila i ies we e es ima ed be ween popula ions o e alua e an isola ion by en i onmen scena io (IBE CLI ). En i onmen al da a was ex ac ed om he 15 CHELSA bioclima ic a iables o each o he 28 sampling si es, as well as om 500 andom poin s co e ing ou s udy a ea o a oid po en ial biases esul ing om only conside ing condi ions a ocal si es. Due o collinea i y among he bioclima ic a iables, we an a PCA using he ade4 R package and summa ized he en i onmen al a ia ion in he h ee fi s axes accoun ing o ca. 82% (PC1 = 43.79%; PC2 = 23.96%; PC3 = 14.31%) o o al a ia ion. The en i onmen al dissimila i y ma ix was ob ained by calcula ing he Euclidean dis ances o PC sco es be ween pai s o sampling si es (O ego e al., 2021). To es he ela i e ole o hos use in popula ion gene ic s uc u e in C. cac o um, a hos plan dis ance ma ix (IBH) was buil conside ing he hos species whe e mo hs we e collec ed. We cons uc ed a bina y ma ix coded wi h each hos species pe sampling si e, and an a PCA using he ade4 R package o build a dis ance ma ix wi h o al PCs a ia ion. IBH was ob ained by calcula ing he Euclidean dis ances o PC sco es be ween pai s o sampling si es. Rela ionships be ween explana o y dis ance ma ices based on landscape he e ogenei y (IBRs, IBD WTD ,IBE CLI ), and he con empo a y pa e ns o hos use (IBH) wi h he gene ic di e en ia ion be ween C. cac o um popula ions (F ST )we e e alua ed using uni a ia e and mul iple ma ix eg essions wi h andomiza ion using he unc ion MMRR (Wang, 2013)as implemen ed in R. An ini ial ull model was cons uc ed conside ing all significan explana o y e ms iden ified p e iously in uni a ia e analysis, and a final bes -fi model was selec ed using a backwa d-s epwise p ocedu e by p og essi ely emo ing non significan a iables un il all e ained e ms wi hin he model we e significan (Wang, 2013;O ego e al., 2014). The esul was he minimal mos adequa e model o explaining he a iabili y in he esponse a iable, whe e only he significan explana o y e ms we e e ained. Po eda-Ma ı ´nez e al. 10.3389/ e o.2023.1260857 F on ie s in Ecology and E olu ion on ie sin.o g05 2.8 Popula ion gene ic di e si y and clima e/habi a s abili y Gi en ha sui able en i onmen s end o suppo la ge e ec i e popula ion sizes (Ca na al e al., 2009), we es ed whe he he egions exhibi ing high gene ic di e si y we e hose whe e sui able clima ic condi ions emained s able h ough ime, om he LGM o he p esen . To es his hypo hesis, en i onmen al s abili y maps we e cons uc ed by a e aging cu en and LGM sui abili y maps ollowing Soley- Gua dia e al. (2019). Fo each sampling si e, alues o clima e/ habi a s abili y we e ex ac ed o each o he h ee scena ios conside ing (i) only clima ic a iables (S abili y ENV ), (ii) clima e- based habi a sui abili y maps o hos species (S abili y HOST ), and (iii) combining clima ic a iables and clima e-based habi a sui abili y maps o hos species (S abili y ENV-HOST ). Ras e calcula ions we e conduc ed using he as e R package. Nuclea and mi ochond ial di e si y we e cha ac e ized o each o he 28 sampling si es. Gene ic di e si y es ima es we e no malized o h ee indi iduals pe sampling loca ion o a oid po en ial bias esul ing om une en popula ion sample size ( ange: 3–11 indi iduals). Fo nDNA, expec ed he e ozygosi ies (H E ) and nucleo ide di e si y (p) we e calcula ed using he di eRsi y (Keenan e al., 2013) R package, and DNAsp 6 (Rozas e al., 2017), espec i ely. We also es ima ed pand haplo ype di e si y (H d ) o m DNA da a using DNAsp. Co ela ions be ween s abili y alues and es ima es o nuclea (H E ,p) and m DNA (H D ,p) popula ion gene ic di e si y we e es ed using linea eg essions. Longi ude and la i ude we e also included as explana o y ac o s o accoun o po en ial geog aphical clines o gene ic di e si y (Guo, 2012). 3 Resul s 3.1 Genomic da a We success ully geno yped 138 indi iduals o C. cac o um using ddRAD sequencing which a e ep esen a i e o popula ions eeding on bo h exo ic O. ficus-indica (71 indi iduals), and na i e hos plan s O. megapo amica (34), O. iopla ense (11), O. ela a (10), O. bonae ensis (4), O. penicillige a (3), O. quimilo (3), and O. anacan ha (2) along h ee biogeog aphic egions. A e fil e ing s eps, he a e age numbe o pai ed-end eads e ained pe indi idual was 1,859,116, o which an a e age o 1,253,277 eads mapped o he C. cac o um e e ence genome. A e disca ding loci in linkage disequilib ium and unde selec ion, we eco e ed a o al o 3,506 biallelic SNPs wi h an a e age co e age o 30× (Table S2). A o al o 143 indi iduals om he same 28 loca ions we e success ully sequenced o a agmen o 790 bp o he COI gene. The analysis o m DNA da a e ealed 66 haplo ypes, 171 a iable si es, 97 o which we e pa simony-in o ma i e sample-wide. 3.2 Assessing popula ion genomic s uc u e Clus e ing analysis wi h sNMF indica ed K = 6 as he mos likely numbe o gene ic clus e s. While all gene ic g oups we e spa ially s uc u ed (No h, Cen al, Wes e n, Sou heas , Sou hwes and Eas ), no gene ic s uc u e based on hos plan s was obse ed. Popula ion gene ic s uc u e analyses also e ealed a ce ain deg ee o admix u e among geog aphic g oups; some No he n and Cen al indi iduals exhibi ed a ela i ely high deg ee o admix u e wi h bo h he Wes e n and Eas e n clus e s (Figu es 2A–C). In e es ingly, No he n, Cen al and Wes e n popula ions appea ed o o m a unique clus e (Cen al lineage om he ea e ) when assuming K = 3, whe eas bo h Sou he n (Sou h lineage om he ea e ) and Eas e n (Eas lineage om he ea e ) popula ions emained as sepa a ed panmic ic g oups. Wi h K = 4, he Sou he n popula ion appea ed subdi ided in o wo geog aphic uni s; Sou heas e n and Sou hwes e n. Finally, wi h K = 5, No he n popula ions appea ed sepa a ed om he Cen al ones (Figu e S2). In he PCA he Eas e n lineage appea s sepa a ed om he Cen al and Sou h lineages along PC1, whe eas he Sou h lineage can be dis inguished om he Cen al lineage along PC2 (Figu e 2D). The median-joining ne wo k ob ained wi h m DNA sequence da a e ealed s ong spa ial gene ic s uc u e consis en wi h he clus e s iden ified wi h nDNA da a (Figu e 2E). The Cen al haplog oup iden ified wi h a s a -like pa e n had one o he mos equen haplo ypes and was sepa a ed by se en mu a ional s eps om he Eas e n haplog oup, whe eas a leas ou mu a ional s eps sepa a ed he Cen al and No he n haplog oups. Wes e n, sou heas e n and Sou hwes e n haplog oups appea ed close o he Cen al haplog oup sepa a ed by ew mu a ional s eps (Figu e 2E). 3.3 Species ee econs uc ion Phylogene ic ela ionships in e ed wi h SNAPP we e consis en wi h he hie a chical spa ial gene ic s uc u e obse ed wi h sNMF (Figu es 2B,C;Figu e S2). The mos ances al spli co esponded o he sepa a ion o Sou h and Eas lineages om he Cen al lineage. Sou hwes e n and Sou heas e n clus e s, as well as Wes , Cen al and No h clus e s di e ged subsequen ly om hei espec i e lineages, as obse ed in clus e ing analyses. 3.4 Demog aphic in e ence using coalescen -based simula ions The scena io conside ing ull asymme ic in e lineage gene flow (model 8) was he mos suppo ed among he nine demog aphic models es ed wi h Fas simcoal2 (Table S3). Unde his model, di e gence among he h ee main lineages o C. cac o um was es ima ed o ha e occu ed du ing he La e Pleis ocene (Figu e 3; Table 1). Specifically, his model in ol es an ini ial spli (T DIV1 )o he Cen al lineage ~75 kya (conside ing wo gene a ions pe yea in na i e hos species, as epo ed by Va one e al., 2014), and a mo e ecen spli (T DIV2 ) sepa a ing he Sou h and Eas lineages occu ing ~19 kya. Demog aphic simula ions es ima ed an e ec i e popula ion size close o 250,000 o he Cen al lineage (NE CENTER ), which was six imes highe han he alue ob ained o he Eas lineage (NE EAST = 41,551). Es ima es o gene flow a es a ied among lineages wi h he highes a e o con empo a y gene Po eda-Ma ı ´nez e al. 10.3389/ e o.2023.1260857 F on ie s in Ecology and E olu ion on ie sin.o g06 FIGURE 3 Visual ep esen a ion o he mos suppo ed demog aphic model (Model 8) es ima ed using Fas simcoal2. Pa ame e es ima es included iming o popula ion di e gence (T DIV1 and T DIV2 ), his o ical and con empo a y e ec i e popula ion sizes (NE ANC1 and NE ANC2 ), and asymme ic a es o gene flow (M SC ,M CS ;M ES ,M SE ;M CE , and M EC ). A B D E C FIGURE 2 Popula ion gene ic s uc u e and phylogene ic ela ionships o he C. cac o um popula ions. (A) Pie cha s ep esen he a e age ances y coe ficien o indi iduals belonging o each o 28 sampling si es dis ibu ed along h ee biogeog aphic egions. (B) Species ee econs uc ed by SNAPP; numbe s in nodes deno e pos e io p obabili y. (C) Resul s o indi idual assignmen in gene ic clus e s using sNMF. Ve ical ba s ep esen he ances y coe ficien o each indi idual o he co esponding gene ic clus e . (D) PCA using he wo majo axes o genomic a ia ion. (E) Median- joining ne wo k ob ained wi h m DNA da a. Image shows an adul male o C.cac o um.. Po eda-Ma ı ´nez e al. 10.3389/ e o.2023.1260857 F on ie s in Ecology and E olu ion on ie sin.o g07 flow es ima ed o occu om he Cen al o he Eas lineage (M CE = 1.50 × 10 −5 mig a ion a e by gene a ion) and he lowes om he Cen al lineage o he Sou h lineages (M CS = 2.71 × 10 −8 mig a ion a e by gene a ion) (Table 1). 3.5 Ecological niche modeling h ough ime Ecological niche models (ENMs) o he ocal species (C. cac o um) and hos species exhibi ed an o e all good pe o mance as sugges ed by he ela i ely high AUC sco es (Tables S4,S5). Sui abili y maps o hos species we e conco dan wi h hei espec i e cu en dis ibu ions (Figu e S3). Among hos species, he a eas o highe sui abili y o Opun ia megapo amica we e p edic ed in cen al and sou he n A gen ina; o O. ela a in cen al- eas e n, and o O. iopla ense in eas e n A gen ina (in he a ea o influence o he Pa ana-La Pla a i e basin). Opun ia anacan ha and O. quimilo ha e mo e es ic ed sui abili y a eas in no hwes e n and cen al A gen ina, espec i ely. P ojec ions du ing he LGM indica ed changes in hos dis ibu ion (Figu e S3) ha sugges an inc eased habi a sui abili y o almos all hos plan s, excep o O. anacan ha whose sui abili y a eas ha e been es ic ed o he cen al-eas e n egion du ing he LGM. Ins ead, habi a sui abili y o O. megapo amica and O. ela a appea s ex ended in o he no he n egion, whe eas he p ojec ions indica e an ex ended dis ibu ion in cen al A gen ina o O. quimilo. Finally, sui able habi a s o O. iopla ense ex ended o no he n and cen al A gen ina. Models o C. cac o um (Se 1: Clima ic CURR ,Se 2:Hos CURR ,Se 3: Clima ic-Hos CURR ) indica ed simila p edic ed dis ibu ions (Figu e S4). Howe e , he mos suppo ed model o C. cac o um acco ding o AIC alues was one ha conside ed as inpu a iables he p edic ed dis ibu ion o he hos plan s (Hos CURR ). Addi ional in o ma ion on he esul s o ENMs o bo h C. cac o um and hos species a e de ailed in Tables S4,S5. The p esen -day dis ibu ion in e ed by he mos suppo ed scena io (Hos CURR ) was in line wi h he con empo a y ange o he mo h. Unde his model, highe sui abili y a eas we e p edic ed in bo h Chaco and no h o Pampa biogeog aphic egions whe eas a eas o lowe sui abili y we e p edic ed in he sou he n po ions o he Pampas and he Mon e biogeog aphic egions (Figu e 1B). O he biogeog aphic p o inces (Yungas, Puna and P epuna) had ex emely low sui abili y alues o C. cac o um. P ojec ions o ENM o LGM p edic ed high en i onmen al sui abili y along he Chaco biogeog aphic p o ince o he mos sou he n po ion o Pampas, wi h a eas o low sui abili y loca ed in eas e n A gen ina, a confluencezoneo LaPla a i e basin (Figu e 1C). 3.6 Landscape genomic analyses Es ima es o gene ic di e en ia ion (F ST ) among sampling si es anged om 0.023 o 0.448 o nDNA da a (Table S6). Uni a ia e ma ix eg essions indica ed ha nuclea gene ic di e en ia ion was significan ly co ela ed wi h all dis ance ma ices excep wi h he cu en pa e n o hos plan use (IBH) consis en wi h clus e ing analysis, and bo h explana o y e ms o landscape he e ogenei y: Clima ic CURR and Clima ic-Hos CURR (Table S7). Ye , only Clima ic-Hos LGM was significan ly e ained in he bes -fi model in he backwa d-s epwise p ocedu e (Table 2). 3.7 Popula ion gene ic di e si y and clima e/habi a s abili y We ound highe nuclea and mi ochond ial di e si y in No he n, Cen al and Eas e n popula ions, han in Sou he n TABLE 1 Pa ame e s in e ed om coalescen simula ions wi h Fas simcoal2 unde he bes -suppo ed demog aphic model (Model 8; Figu e 3). Model 8 95% Confidence in e al Pa ame e Poin es ima e Lowe bound Uppe bound NE ANC1 3,504,318 2,576,620 3,858,949 NE ANC2 11,746 3,773 21,102 NE CENTER 245,957 208,570 380,730 NE EAST 41,551 35,596 54,151 T DIV2 75,513 61,524 120,375 T DIV1 19,207 17,202 23,310 M ES 1.18 × 10 −6 8.86 × 10 −7 1.75 × 10 −6 M SE 4.68 × 10 −8 5.41 × 10 −11 1.65 × 10 −6 M CS 2.71 × 10 −8 3.55 × 10 −11 6.83 × 10 −8 M SC 6.37 × 10 −6 4.37 × 10 −6 7.32 × 10 −6 M CE 1.50 × 10 −5 1.23 × 10 −5 1.67 × 10 −5 M EC 2.98 × 10 −5 3.78 × 10 −11 1.42 × 10 −6 Fo each pa ame e , poin es ima e and lowe and uppe 95% confidence in e als a e shown. T DIV1 and T DIV2 , popula ion di e gence; NE ANC1 and NE ANC2 , his o ical and con empo a y e ec i e popula ion sizes; M SC ,M CS ,M ES ,M SE ,M CE , and M EC , asymme ic a es o gene flow among lineages. Po eda-Ma ı ´nez e al. 10.3389/ e o.2023.1260857 F on ie s in Ecology and E olu ion on ie sin.o g08 popula ions (Figu e 4;Table S8). Linea eg ession analyses e ealed ha nuclea gene ic di e si y was significan ly co ela ed wi h la i ude bu no wi h longi ude (Table 3). Likewise, es ima es o nuclea gene ic di e si y o bo h H E and pwe e posi i ely co ela ed wi h he h ee en i onmen al s abili y es ima es (S abili y ENV , S abili y HOST and S abili y ENV-HOST ) o he LGM (~21 Kya) o he p esen . Con e sely, mi ochond ial gene ic di e si y es ima ed wi h H D and pwe e no co ela ed wi h geog aphic a iables no en i onmen al s abili y (Table 3;Figu e 4). 4 Discussion Nuclea SNPs and m DNA da a e ealed significan popula ion s uc u e in C. cac o um. Landscape genomic analyses p o ided suppo o he hypo hesis ha habi a connec i i y based on a combina ion o clima ic condi ions and habi a sui abili y, influenced by shi s in hos species dis ibu ions du ing LGM, was he main ac o shaping popula ion gene ic s uc u e. Al hough pa e ns o hos species use appea ed o be non influen ial on C. cac o um di e sifica ion, shi s in he dis ibu ion o Opun ia hos s, media ed by clima ic changes du ing he Qua e na y, had a di ec influence on he dis ibu ion and genomic a ia ion o he ocal species. Such shi s in hos dis ibu ions may ha e gene a ed agmen ed anges ha led o educ ions o gene flow among mo h popula ions, p omo ing gene ic di e en ia ion. Ou esul s also suppo ed he habi a s abili y hypo hesis, whe eby popula ions inhabi ing egions o g ea e habi a sui abili y du ing he las 21,000 yea s ha bo ed mo e gene ic di e si y han hose ha pe sis ed in a eas o lowe habi a s abili y. The hypo hesis ha he geog aphically widesp ead cul i a ion o O. ficus-indica acili a ed con empo a y gene flow among o he wise geog aphically dis an popula ions can be ejec ed, on he basis o he esul s o popula ion genomic analyses e ealing limi ed gene ic admix u e and a hie a chical pa e n o popula ion di e en ia ion mainly conco dan wi h geog aphy. 4.1 E ec s o Qua e na y landscape composi ion on popula ions di e en ia ion Gene ic clus e ing analyses e ealed h ee majo lineages defined by geog aphy: Cen al, Eas and Sou h (Figu es 2 and 3). Coalescen - based demog aphic modeling sugges ed ha hese lineages began o di e ge du ing he La e Qua e na y, app oxima ely 75 kya. The ea lies spli in ol ed popula ions o he Chacoan biogeog aphic domain (he ea e “Chaco”), including popula ions om No h, Wes and Cen al A gen ina, om he Pampean biogeog aphic egion (he ea e “Pampa”) ha includes he Sou h and Eas lineages (Figu es 1,2). P io o his spli , he egional con ac ion o sub opical and opical biomes made he Chaco and he Pampa biomes mo e alike (O iz-Jau eguiza and Clade a, 2006). Subsequen ly, du ing he Qua e na y, majo clima ic shi s ook place in he a ea, leading o he isola ion o he Chacoan xe ic woodlands om he ypical g asslands o he Pampean egion. Such clima ic changes we e p omo ed by he opog aphic econfigu a ion o Andean and sub-Andean Piedmon and he upli o he eas e n o og aphic sys ems (known as he Sie as Pampeanas) associa ed wi h he Pe ipampasic o ogenic a c (O iz- Jau eguiza and Clade a, 2006;Spe anza e al., 2007;Cala ayud- Masca ell e al., 2022). These e en s p oduced a ain-shadow e ec ha esul ed in he ex eme xe ic condi ions ha p e ailed in his a ea. Clima e changes du ing he Qua e na y, which included cold, d y glacial cycles al e na ing wi h wa m, mois in e glacial pe iods, a ec ed hese egions causing he expansion and con ac ion o open woodlands and g asslands (O iz-Jau eguiza and Clade a, 2006). This dynamic likely p omo ed changes in flo is ic composi ion o such egions, a ec ing he dis ibu ions o hos species and consequen habi a sui abili y o C. cac o um,p omo ing he di e gence be ween he Cen al lineage (Chaco) and Eas , Sou h lineages (Pampean) (Figu es 2,3). The spli be ween Eas and Sou h lineages was es ima ed o ha e occu ed app oxima ely 19 kya. Regional clima e models sugges ha his iming is coinciden wi h an inc ease in p ecipi a ion along he eas e n oo hills o he Andes (Cook and Vizy, 2006;O iz- Jau eguiza and Clade a, 2006). These opog aphical ea u es oge he wi h Qua e na y clima ic dynamics in he egion main ained a agmen ed dis ibu ion o C. cac o um hos species om he Pampa (Mou elle and Ezcu a, 1997). This likely limi ed dispe sal o C. cac o um, es ic ing he Eas lineage o he mo e humid en i onmen s wi h deep e ile soils ypical o he Pampa, and he Sou h lineage o he d ie g asslands o he no he n bo de o he Pa agonian s eppe, bo h egions wi h di e se pa e ns o hos species composi ion (Va one e al., 2014). TABLE 2 Resul s o mul iple ma ix eg ession wi h andomiza ion (MMRR) es ing o he ela ionships be ween popula ion gene ic di e en ia ion (F ST ) and explana o y a iables ep esen ing weigh ed opog aphic dis ance (isola ion by dis ance, IBD WTD ), clima ic dissimila i y (isola ion by en i onmen , IBE CLI ), hos plan dis ance (*IBH), and al e na i e isola ion by esis ance (IBR) scena ios. Va iable R 2 p Explana o y e m LGM ENV-HOST 0.416 16.365 0.001 Rejec ed e ms Clima ic dissimila i y (IBE CLI )−0.476 0.831 Topog aphic dis ance (IBD WTD )−2.058 0.301 IBR –Hos CURR 3.251 0.278 IBR –Clima ic LGM 4.808 0.117 IBR –Hos LGM −5.834 0.087 Fla landscape –IBR NULL −3.654 0.278 IBR scena ios consis o cu en and pas habi a sui abili y based on (i) only clima ic a iables (*Clima ic CURR , and Clima ic LGM ), (ii) clima e-based habi a sui abili y maps o hos species (Hos CURR and Hos LGM ), and (iii) combining clima ic a iables clima e-based habi a sui abili y maps o hos species (*Clima ic-Hos CURR and Clima ic-Hos LGM ). IBR NULL ep esen s an IBR scena io whe e all pixel alues a e equal o 1 (fla landscape). R 2 , eg ession coe ficien ; , -s a is ic; p, significance le el. *Clima ic CURR , Clima ic-Hos CURR , and cu en pa e n o hos use (IBH) we e excluded om he ini ial ull model due o hese e ms esul ed non-significan in hei espec i e uni a ia e analysis. Po eda-Ma ı ´nez e al. 10.3389/ e o.2023.1260857 F on ie s in Ecology and E olu ion on ie sin.o g09