Frugivore-frugtstørrelsesforhold mellem palmer og pattedyr afslører tidligere og fremtidige defaunationspåvirkninger

palmeartsfordeling og frugtstørrelsesdata

artsfordelinger for alle (~2500) palmearter blev opnået fra Verdenschecklisten over Palmer30 (hentet Juni 2015), en kurateret og kvalitetskontrolleret tjekliste, der giver tilstedeværelse-og frugtstørrelsesdata, der er fraværsdata i omfanget af “botaniske lande” som defineret af den internationale arbejdsgruppe om Taksonomiske databaser (TDG, https://www.tdwg.org/) 31. Botaniske lande svarer generelt til lande, men opdeler større lande i stater eller provinser (f.eks. USA) eller bevarer sammenhængende geografiske enheder, der er politisk opdelt (f. eks. Borneo og Ny Guinea) som enkeltenheder31. Vi brugte verdens tjekliste over palmer snarere end artsforekomster (dvs.kun tilstedeværelsesoptegnelser fra Global Biodiversity Information Facility, GBIF), da sidstnævnte er taksonomisk ufuldstændig (næsten 30% af alle palmearter er ikke repræsenteret blandt poster, der er adgang til via GBIF, hentet januar 2020). Derudover er begivenhedsregistre ofte geografisk og taksonomisk forudindtaget, hvilket resulterer i store huller og usikkerheder i Global begivenhedsinformation om planter (inklusive palmer)70.

af de 198 botaniske lande, hvor palmer forekommer, inkluderede vi kun dem i de tre store biogeografiske riger (n = 129), modificeret fra ref. 71: (1) Afrotropics, der omfatter hele Afrika syd for Sahara (n = 51), (2) Neotropics (n = 42), som omfatter Central-og Sydamerika og Caribien, og (3) Indo-Australien (n = 36), som omfatter Indo-Malaya, Australien og Papua Ny Guinea. Botaniske lande i andre biogeografiske regioner (dvs. Palearctic, Nearctic og Oceanien) blev udelukket, da de enten har depauperate palmeflora eller er så fysisk isolerede (f.eks. Stillehavsøerne), at spredningsfiltre kan have uforholdsmæssigt formet både palme-og frugivore-samlinger. Også udelukket var botaniske lande, for hvilke eksisterende pattedyrsfrugivorer ikke længere er til stede, bestemt af fraværet af sparsomme pattedyr efter at have krydset pattedyrs geografiske rækkevidde med den botaniske landeformfil (se nedenfor). Kun et botanisk land blev udelukket af denne grund—Reunion Island, hvor to indfødte arter af frugtflagermus oprindeligt blev udryddet fra øen i det tidlige attende århundrede eller tidligere72.

maksimale og median frugtstørrelsesværdier blev beregnet på tværs af alle palmearter i hvert botanisk land. For at kvantificere frugtstørrelse for hver art brugte vi PalmTraits database v1.0 (ref. 29), som har samlet palmefunktionelle trækdata for prøver fra vilde populationer fra den primære litteratur, monografier, artsbeskrivelser og herbariumprøver. Specifikt brugte vi gennemsnitlig frugtlængde (snarere end bredde eller diameter), da dette var den mest konsekvent registrerede frugtstørrelsesegenskab i litteraturen29. Maksimal frugtstørrelse blev defineret som den 95.percentil af frugtstørrelser af palmer i hvert land for at reducere indflydelsen fra outliers. Da oplysninger om palmefrugtlængde ikke var tilgængelige for 519 palmearter (~20% af den globale palme-artsrigdom), udfyldte vi huller i frugtlængden af disse arter ved hjælp af gennemsnittet af congeneriske arter. Dette var ikke muligt for en art, Butyagrus nabonnandii, en naturligt forekommende steril hybrid, som vi således udelukkede fra vores Analyse. Vi udelukker desuden tre arter fra overvejelse-kokosnød (Cocos nucifera) og coco De mer (Lodoicea maldivica) og nipa—palmen (Nypa fruticans) – da de har frugter, der sjældent spredes af dyr24. Maksimale og median frugtstørrelser på tværs af botaniske lande var ens, når huller i egenskabsdata blev udelukket (Pearson korrelationskoefficient >0,99; supplerende Fig. 1).

Frugivore artsfordeling og kropsstørrelsesdata

for at bestemme effekten af end-Pleistocene frugivore-udryddelser på nutidens makroøkologiske forhold definerede vi to forskellige faunale scenarier, der repræsenterer frugivore-samlinger, som de er i dag, og som de var i det sene Pleistocæn—de “nuværende” og “nuværende naturlige” scenarier efter Faurby et al.9. Det” nuværende ” scenarie betragter kun eksisterende frugivores, mens det kontrafaktiske “nuværende naturlige” scenarie betragter alle frugivores siden den sidste Interglacial (130 Kya), inklusive alle eksisterende arter såvel som for nylig uddøde takser. Mens klimaændringer kan have spillet en rolle i udryddelsen af nogle af disse frugivores39,73, er størstedelen af udryddelserne af disse arter tilskrevet menneskelig ankomst og aktiviteter39,40,74,75.

“nuværende” og “nuværende naturlige” rumlige fordelinger af både eksisterende og sene Pleistocene–til-Holocene uddøde pattedyrsfrugivorer var baseret på taksonomi og geografiske rækkevidde kort som samlet i PHYLACINE v1.2 datasæt76. Områdekort var tilgængelige som Behrmann cylindriske rastere med lige areal, som vi derefter overlejrede med botaniske landpolygoner for at generere pattedyrsartslister, der kan sammenlignes med datasættet for tilstedeværelse af palm–fravær. Aktuelle områdekort over eksisterende pattedyr i PHYLACIN-datasættet var baseret på IUCN-områdekort (http://www.iucnredlist.org/) med minimal modifikation (f. eks., irettesættelse). Imidlertid blev de nuværende naturlige geografiske områder af eksisterende skatter systematisk evalueret, og ændringer blev foretaget, hvis der var stærke tegn på menneskeskabt påvirkning (se Supplerende i refs. 9,76). Dette omfattede også menneskelige medierede områdeudvidelser, som IUCN-områdekort for nogle takser ikke korrigerer for. For uddøde skatter, nuværende naturlige intervaller blev estimeret baseret på fossile forekomster og fordelingen af eksisterende skatter,som den uddøde skat ofte forekom inden for fossile samlinger9, 76. Dette blev udført ved at identificere gitterceller, der indeholdt mindst 50% af den eksisterende taksa (med en fossil rekord), der opstod sammen med en uddød takson i fossile samlinger/steder (datakilder og yderligere detaljer om metode i refs. 9,76). Områdekort med nul-besatte gitterceller (f.eks. arter, hvis intervaller ikke er kortlagt af IUCN) blev udelukket fra vores Analyse. Pattedyrstjeklister for hvert botanisk land blev også manuelt inspiceret for at rette fejl, der undertiden blev introduceret på grund af den grove opløsning af områdekort (f. eks. på øer, der ligger tæt på fastlandskystlinjer).

vi brugte datasættet Pattedyrdiet77 til at identificere eksisterende arter, der enten primært eller sekundært var sparsomme. Synonymer af binomiale navne i Pattedyrdiætdatasættet blev korrigeret ved hjælp af IUCN-rødlisten som en taksonomisk ref. 78. Artsnavne, der ikke kunne forenes med taksonomien i PHYLACIN-datasættet, blev udeladt, ligesom taksa, hvis intervaller ikke overlapper hinanden med nogen eksisterende palmeafgift. Vi bestemte, i hvilken grad uddøde pattedyr var sparsomme baseret på indirekte beviser, da sparsommelighed ikke ligefrem er forbundet med morfologiske træk såsom tandproteser79. I erkendelse af denne kompleksitet overvejede vi tre frugivore klassifikationer med varierende grad af konservativitet. Under vores” liberale ” klassificering klassificerede vi enhver uddød takson som en frugivore, hvis taksonen overvejende var planteædende. Denne definition inkluderer obligatoriske græssere, der kan tjene som frødispergeringsmidler, men sandsynligvis kun har en mindre frugtkomponent i deres kost sammenlignet med ikke-græssere. Under vores” konservative ” klassificering inkluderede vi desuden kun en uddød takson som en frugivore, hvis den tilhørte familier, hvor 50% af de eksisterende arter overvejende er sparsomme. Denne definition, i modsætning til den liberale klassifikation, udelukker de fleste græsere, men også mange potentielle frugivorer, hovedsageligt fordi eksisterende slægtninge kan være dårlige analoger, da uddøde slægtninge kan variere meget i kropsstørrelse og dermed deres fodringsøkologi, såsom i store diprotodont (kænguruer) og pilosan (jorden dovendyr) planteædere (supplerende Tabel 7). For vores mellemliggende” standard ” klassificering, vi inkluderet taksa fra familier, der har en grænseandel af frugivores (til 40%). Vi inkluderede desuden uddøde takser, der kommer fra familier eller ordrer, der ikke har eksisterende arter (f.eks. litopternider og toksodoner) eller kommer fra familier, hvor eksisterende takser kan være dårlige analoger til kost af uddøde takser, hvis de tilgængelige paleontologiske beviser (f. eks. Under denne definition, jorden sloths80, litopternids81,82, toksodons80, 83,nogle heste—medlemmer af slægten Hippidion, men ikke Ækve80, nogle cingulids80, 84, nogle diprotodonter (makropodider, men ikke vombatidae og Diprotodontidae)85, og murid gnavere blev betragtet som formodede frugivores (supplerende Fig. 4, Supplerende Tabel 7).

vi beregnede median og maksimal frugivore kropsstørrelse (her defineret som den 95.percentil) for hvert botanisk land under det nuværende scenario. For det nuværende naturlige scenario blev dette udført for de tre uddøde frugivore klassifikationer separat. Kropsstørrelse (dvs.kropsmasse) information blev opnået fra PHYLACINDATASÆTTET. I alt betragter vi de nuværende og nuværende naturlige intervaller på henholdsvis 1759 og 1772 eksisterende skatter ud af 1930–eksisterende skatter i pattedyrdietdatasættet efter taksonomisk forsoning. I alt 177 Pleistocæn formodede frugivores og 23 uddøde eller uddøde i naturen blev desuden inkluderet i vores nuværende naturlige scenario.

miljødata

da frugtstørrelse også kan afspejle storskala variation i klima og/eller tidligere klimaændringer37 blev værdierne for det gennemsnitlige nuværende klima og klimaændringer siden det sidste glaciale maksimum (LGM, ~21 Kya) beregnet for hvert botanisk land i ArcGIS (version 10.1, ESRI, Redlands, CA, USA). Dette blev gjort ved at overlejre botaniske landpolygoner med de tilsvarende klimatiske rastere og tage gennemsnittet af overlappende rasterceller. Vi brugte CHELSA database86 til både nuværende og LGM klima (v1.2, 30 bue anden opløsning), som har fordele i forhold til Verdenklim data87 da CHELSA-algoritmen anvender korrektioner for fine orografiske effekter på Nedbør.

for det nuværende klima fokuserede vi på seks bioklimatiske variabler: gennemsnitlig årlig temperatur, temperatur sæsonbestemthed, gennemsnitlig temperatur i det koldeste kvartal, årlig nedbør, Nedbør sæsonbestemthed og nedbør i det tørreste kvartal. Fordi nogle af disse bioklimatiske variabler er meget kollinære med hinanden, og da vi primært var interesserede i den relative betydning af det nuværende klima og ikke indflydelsen af specifikke klimatiske variabler i sig selv, opsummerede vi variationen på tværs af de seks bioklimatiske variabler ved hjælp af en hovedkomponentanalyse (PCA) og brugte de første tre hovedkomponenter som klimatiske forudsigelsesvariabler. Dette blev udført på globalt plan og for hver biogeografisk region (“Afrotropics”, “Neotropics” og “Indo-Australien”) separat. De første tre PCA-akser (“Climate PC”) forklarede over 90% af variabiliteten inden for de bioklimatiske variabler blandt botaniske lande på både global og regional skala (supplerende tabel 8).

for tidligere klima beregnede vi størrelsen af ændringen i årlig nedbør (“LGM Prec. Anom.”) og gennemsnitlig årlig temperatur (“LGM Temp. Anom.”) siden det sidste glaciale maksimum for hvert botanisk land. Dette blev beregnet ved at trække den nuværende årlige nedbør og den gennemsnitlige årlige temperatur fra ensemblets gennemsnit af den gennemsnitlige årlige nedbør og den årlige temperatur på tværs af seks forskellige GCM-modeller fra Paleoclimate Modeling Intercomparison Project88 (https://pmip3.lsce.ipsl.fr/), der er blevet behandlet af CHELSA-algoritmen: CCSM4, CNRM-CM5, CESS-FGOALS-g2, IPSL-CM5A-LR, MIROC-ESM og MR-CGCM3. Højere værdier repræsenterer således en højere størrelse af klimaændringer siden det sidste glaciale maksimum.

statistiske analyser

for at evaluere forholdet mellem frugivore krop og frugtstørrelse monterede vi begge OLS lineære modeller af både median og maksimal størrelse ved hjælp af nutidens (“Climate PC1”, “Climate PC2” og “Climate PC3”) og tidligere klima (“LGM Prec. Anom ” og ” LGM Temp. Anom.”) som forudsigelsesvariabler såvel som henholdsvis maksimal og median frugivore kropsstørrelse. Skøn over kropsstørrelse for nuværende og nuværende naturlige scenarier, og hver biogeografisk region blev modelleret separat. Krop-og frugtstørrelsesmål blev log-transformeret for at forbedre normaliteten i rester, og alle forudsigelsesvariabler blev centreret og skaleret ved at standardisere dem til enhedsvarians.

rumlig autokorrelation blev regnet med anvendelse af rumlige autoregressive (SAR) modeller. Vi brugte Sar-fejlmodeller (SARerror), som har vist sig at være mere robuste end andre SAR-modeltyper89. SAR-modeller forsøger at redegøre for rumlig autokorrelation mellem enheder inden for en given kvarterstruktur. Imidlertid, på grund af den uregelmæssige afstand og størrelse af botaniske lande i vores undersøgelse, adjacency – eller afstandsbaserede kvarterer var ikke passende. I stedet definerede vi kvarteret i et givet botanisk land ved hjælp af en indflydelsessfære tilgang. Indflydelsessfæren for hver fokal enhed er defineret som en cirkel af radius svarende til fokalenhedens Afstand til centroid af sin nærmeste nabo. Hvor indflydelsessfæren af to enheder overlapper hinanden, betragtes de to enheder som naboer. Denne kvarterdefinition ser ud til at være mere passende sammenlignet med afstandsbaserede eller nærmeste nabo-definitioner, da den tager højde for størrelsen og geometrien i hvert botanisk land. Rumlige vægte var rækkestandardiserede. Kvarterer, rumlige vægte og rumlige SAR-modeller blev implementeret ved hjælp af ‘spatialreg’ r–pakken v. 1.1-3 (ref. 90), og vi testede desuden for enhver resterende rumlig autokorrelation ved hjælp af Morans i-test.

vi brugte en multi-model-gennemsnitlig tilgang til at estimere effektstørrelser for hver forudsigelsesvariabel på tværs af et sæt kandidatmodeller91. Den primære fordel ved modelgennemsnit er, at den tegner sig for modelusikkerhed på tværs af et sæt kandidatmodeller,i modsætning til mere traditionelle trinvise modelvalgsmetoder, der forsøger at identificere den bedste model92, 93. Vi definerede vores kandidatmodeller som sæt af almindelige mindst kvadrater (OLS) lineære regressionsmodeller med log-transformeret palmefrugtstørrelse som responsvariabel, hvor hver model indeholder en anden kombination af log-transformeret frugivore kropsstørrelse, nutidens klima og tidligere klimaændringer som kovariater. Effektstørrelserne for hver forudsigelsesvariabel blev derefter beregnet ved gennemsnitskoefficientværdier på tværs af alle kandidatmodeller vægtet med mængden af statistisk støtte, som hver af disse modeller har (dvs.Akaike-vægten)91. For at reducere virkningen af modelvalgsforstyrrelse på model-gennemsnitlige koefficienter blev koefficientværdien for en forudsigelsesvariabel i en model, hvor den er fraværende,indstillet til nul91, 92. Dette sikrer, at effektstørrelserne af forudsigere, der kun er i modeller med lave Akaike-vægte, er nedvægtet92. Akaike-vægte blev beregnet ved hjælp af AIC-værdier, der er korrigeret for lille prøvestørrelse (AICC). Korrelation mellem forudsigelsesvariabler var generelt lav (supplerende Fig. 3), og variansinflationsfaktorer for alle forudsigelsesvariabler i fulde OLS-modeller (dvs.modeller indeholdende alle forudsigelsesvariabler) var for det meste under 4. For at reducere indflydelsen af kollinearitet på modelgennemsnitlige estimater beregnede vi desuden modelgennemsnitlige koefficienter, hvor koefficienter er blevet standardiseret af deres delvise standardafvigelser før gennemsnitet32,93 (supplerende tabeller 2 og 5). Model gennemsnit af både OLS og SAR modeller blev udført ved hjælp af ‘dredge’ og ‘model.avg ‘funktioner fra’ MuMIn ‘ r pakke v. 1.42.1 (ref. 94). Standardisering af koefficienter for OLS-modeller før modelgennemsnit blev udført ved at indstille ‘beta’ – argumentet i ‘ mudder ’til ‘ partial.sd’.

ud over Model-gennemsnitlige standardiserede effektstørrelser sammenligner vi den forklarende effekt af kropsstørrelse under “nuværende” og “nuværende naturlige” scenarier. Fordi summen af Akaike-vægter91 har vist sig at være en inkonsekvent metrisk af variabel betydning93,95, beregnede vi i stedet andelen af varians, der blev nedbrudt på tværs af forudsigelsesvariabler96, som implementeret i r-pakken ‘relaimpo’ (v2.2-3)97. For enkelhedens skyld blev dette udført på fulde OLS-og SAR-modeller og ikke på tværs af alle modeller testet i model-gennemsnitsproceduren. For at bestemme andelen af varians forklaret af forudsigelsesvariabler i SAR-modellerne efter rumlig autokorrelation er taget i betragtning, en rumlig model med kun et skæringspunkt blev først monteret, og andelen af varians i resterne af den kun intercept-rumlige model forklaret og en fuld OLS-model med alle forudsigelsesvariabler blev derefter monteret på resterne. I betragtning af at modeller til nuværende og nuværende naturlige scenarier kun adskiller sig i kropsstørrelsesforudsigelsesvariablen, vi leverer både R2 og pseudo–R2 til fulde OLS-og SAR-modeller under begge scenarier.

simulering af frugivore-udryddelse

for at projicere ændringer i frugtstørrelse under et fremtidigt defaunationsscenarie simulerede vi frugivore-samlinger givet sandsynligheder for, at arter inden for forskellige IUCN-Rødlistekategorier vil uddø globalt i de næste 100 år98,99. Vi brugte to forskellige sæt udryddelsessandsynligheder til at undersøge omfanget af udryddelsesrisiko, som pattedyrsfrugivorer potentielt står over for (supplerende Tabel 9). Repræsenterer en høj defaunation scenario, vi brugte udryddelse sandsynligheder afledt af Davis et al.99, der oversatte kriterier for rødliste100 for truede kategorier (dvs.kritisk truet (CR), truet (EN) og sårbar (Vu)) til udryddelsesrater for disse kategorier og ekstrapolerede disse estimater for at udlede satser for ikke-truede kategorier (dvs. næsten truet (NT) og mindst bekymring (LC)) under den antagelse, at satser vil stige eksponentielt med alvorligheden af kategorien. Udryddelsessandsynligheder for hver kategori over en 100-årig tidsramme blev derefter beregnet ved hjælp af disse satser under en konstant udryddelsesproces ved hjælp af følgende ligning: I = 1-eksp (- rit), hvor ri er udryddelsesgraden for en takson i IUCN kategori i og t er lig med tiden (supplerende Tabel 9). Disse værdier er potentielt høje, da de anvender udledte estimater af udryddelsesrisiko ensartet, uanset de kriterier, der er anvendt på en takson, og overvejer ikke eksplicit variation i livshistorien og dens virkning på udryddelsesrisiko. Som et potentielt lav defaunationsscenarie udviklede vi en kontinuerlig Markov-kæde (CTMC) model i form af en overgangshastighedsmatrice (k), hvis elementer beskriver øjeblikkelige overgangshastigheder mellem tilstødende Rødlistekategorier:

$${\bf= \ left$$
(1)

hvor ri, J er den øjeblikkelige overgangshastighed fra IUCN kategori i til kategori j.

vi monterede vores ctmc-model ved hjælp af maksimal sandsynlighed for to datasæt, der kompilerede ændringer (eller mangel på dem) i rødlistekategorier på tværs af alle IUCN-evaluerede pattedyrarter over en 12-årig periode101, samt hovdyr og kødædere over en 33-årig periode102. Alle arter, der er blevet klassificeret som uddøde i naturen (f.eks.) og CE-arter, der er markeret som muligvis uddøde eller uddøde i naturen (dvs. CR (PE) eller CR (KIRKESTOL) kategorier) blev grupperet i en enkelt kategori, der repræsenterer formodede udryddelser. For enkelhedens skyld er denne kategori blevet betegnet som eks i EKV. (1). Vi antog desuden, at ændringer i rødliste kun kan forekomme mellem tilstødende rækker ved et givet uendeligt lille tidsinterval, Kurt. For eksempel skal en art af LC-bevaringsstatus for at overgå til en mere alvorlig (EN) status gennemgå (dog ikke udelukkende) status af mellemliggende sværhedsgrad (dvs.NT og VU). For Di Marco et al.102 datasæt, fordi der ikke blev registreret nogen formodet uddød skat for at være vendt tilbage til en mindre alvorlig rødlistekategori, blev cr sat til 0 for at forhindre overfitting. Dette betød, at eks-staten var en absorberende tilstand (dvs. denne kategori, men ikke fra denne kategori). Ved at antage, at satserne ikke varierer gennem tiden, kan overgangshastighedsmatricen eksponentieres for at give sandsynligheder for overgang mellem kategorier for enhver vilkårlig tid t (dvs. overgangssandsynlighedsmatricen, P):

$${\bf{P}} (t) = \ eksp ({\bf{s}} t)\\ = \ left,$$
(2)

hvor Pi, j (t) er sandsynligheden for, at en taksa starter med status i ender med status j efter en given tid t.

vi antog derefter ændringer (eller mangel på dem) i rødlistekategorier for hver art for at være en uafhængig realisering af en CTMC-proces. Dermed, sandsynligheden for at observere et givet sæt overgangshastighedsændringer givet en tidsramme t, kunne beregnes som den sammensatte overgangssandsynlighed på tværs af arter

$${\mathcal{L}} ({\bf{K}}, t / {\bf{D}})=Pr ({\bf{D}}| {\bf{K}},t)={\prod ^{i}} {\prod ^{j}} {P}_{i, j} (t)\gange {n}_{i, j},$$
(3)

hvor ni, j er antallet af takser, der startede i rødliste kategori i og sluttede i rødliste kategori j efter tid t

Ratestimater for begge datasæt blev opnået ved at maksimere EKV. (3) Brug af funktionen ‘optim’ i R. Robusthed af hastighedsestimater blev testet ved hjælp af parametrisk bootstrapping (supplerende Fig. 5). Sandsynligheden for udryddelse for hver kategori (dvs.sandsynligheden for overgang til AB-kategorien over en 100-årig periode) blev derefter beregnet ved hjælp af gennemsnittet af estimater for maksimal sandsynlighed for de to datasæt ved hjælp af EKV. (2) (Supplerende Tabel 9). Disse værdier er potentielt lave estimater af udryddelsesrisiko, da de er baseret på tidligere ændringer i Rødlistestatus, mens menneskeskabte virkninger kan intensiveres i fremtiden14.

udryddelse blev stokastisk simuleret for alle pattedyrsfrugivorer 1000 gange, og de maksimale kropsstørrelser (95.percentil) i botaniske lande blev registreret. Arter kategoriseret som datamangel blev konservativt antaget at have den samme udryddelsessandsynlighed som dem med LC-status. For simuleringer, hvor alle pattedyrsfrugivorer i et botanisk land uddøde, blev maksimal kropsstørrelse sat til nul. Middelværdier for maksimale kropsstørrelser på tværs af simuleringer blev derefter brugt sammen med globale OLS-modeller med maksimal frugtstørrelse parametreret ved hjælp af alle klimatiske variabler og den aktuelle maksimale kropsstørrelse på frugivore-samlinger for at beregne forventede ændringer i frugtstørrelse under et standardscenarie. Denne analyse blev gentaget ved anvendelse af absolutte maksimale værdier for frugt og kropsstørrelse (supplerende Fig. 6).

Rapporteringsoversigt

yderligere information om forskningsdesign er tilgængelig i Nature Research Reporting Summary, der er knyttet til denne artikel.

Skriv et svar

Din e-mailadresse vil ikke blive publiceret.

Previous post 5 kraftfulde sange, der vil trøste folk i deres sorg
Next post 1