Ecology

Based on the extensive camera trap study, the very first information about the occupancy, trapping rate, activity pattern, group size, social structure and vital rates of the critically endangered Western Derby eland in its last refugium, the NKNP in Senegal, is presented here. The first estimation of abundance since 2006 is also provided7.Spatiotemporal behaviour pattern of WDE in the parkThe results of the CT survey in the NKNP highlight the substantially lower occupancy and trapping rate of WDE in comparison to other large ungulates in the park. According to the current results, the WDE occupied less than 5% of the park area during the dry season, being exclusively within the zone of Mont Assirik, and more specifically the Mansa Fara marsh, which can be thus designated as the core area of the WDE distribution. The trapping rate of the roan antelope, which is considered the most abundant antelope species in the NKNP, was 4.04, i.e. more than 11 times higher than that of WDE in the present study. Even the Western hartebeest, which is considered a rare species in the NKNP, had a trapping rate of 0.61, which is ca. twice as high that that of the WDE (see Rabeil et al.8 for further details, and additional ungulate species).In the zone of Mont Assirik, the trapping rate of WDE increased to 2.42, but the trapping rates of other antelope species remained higher still (4.6 for roan antelope and 3.39 for Western hartebeest8). The WDE distribution is therefore strongly localised within an area which seems to also be attractive for the other species, including the incidental records of elephant. The Mont Assirik zone, and more specifically the Mansa Fara marsh area, is therefore the crucial zone within the park for the WDE, and it appears to support a larger number of other antelope species as well. This zone should therefore be considered as a key conservation area, potentially very sensitive to targeted poaching, and thus crucial for efficacy of targeted law enforcement actions.When looking at the diurnal activity pattern, the WDE were active before midnight, approximately 3 h after sunset, in the morning, approximately 2 h after sunrise, and then again in the afternoon, with the peak activity during the hottest part of the day. This activity pattern is different from the typical bimodal activity pattern, which has peaks at dawn and dusk, as reported for most African grazing and browsing herbivores, seen as a behavioural thermoregulation strategy to avoid heat stress41,42,43. Instead, the WDE, being a large body-sized browsing antelope19,44, must stay active throughout the day to seek discretely distributed food, and fulfil foraging requirements by feeding while moving. The WDE appears to be well-adapted to tolerate such high temperatures, similar to kudu45, roan46, and giraffe47. Such behaviour pattern enable the law enforcement patrols, as well to tourists, to detect herds of WDEs and monitor them, thereby increasing their protection against poaching.Individual identification and recapturesThe individual identification of animals was more successful during the daytime, as the light conditions mostly did not allow for the proper visualization of the stripes on the flanks during the night captures (as similarly reported in Jůnek et al.18). When the ID is targeted to be successful during the night (as for the leopards and tigers), the camera traps are often set to the video mode to ensure a higher possibility of identification48. However, the activity of the WDE is not predominantly nocturnal, and the captures were distributed over both the daylight and night hours, and therefore the results are considered representative for the whole period.The AD animals were more likely to be identified in the present study because of their larger body size, resulting in better visibility of their stripes. The higher identification rate of larger individuals also likely contributed to the higher probability of recaptures, which were only recorded for individuals of 2Y and older.Overall, the identification success rate was comparable, maybe even slightly higher, than the previous camera trap study performed on the Eastern Giant eland in Chinko, CAR, specifically in the dataset from the dry season13, which corresponds to the observation period in the present study as well.In the NKNP, recaptures of individuals were recorded, whereas there were none reported in Chinko13. The recapture rate of the WDE in NKNP, with mostly short distances between the capture-recapture sites, even after the long-time gaps between the captures, confirm again that the WDE likely inhabit a relatively limited area of the park.Group size and social structureThe mean group size recorded in the NKNP during the present study was slightly larger than that within Chinko; however, the maximum group size was smaller in NKNP (32 vs. 41 individuals). Mixed herds were the largest in terms of the number of individuals, in both studies. The average group size has been reported as 20–30 individuals49, but Derby elands may form large herds of over 100 individuals in the late dry season14. Similarly, a large herd was reported within NKNP in 2006, having 69 individuals7, and a herd of around 60 WDE was also recently reported by patrols in 2020 (GIE Niokolo, personal communications). It is important to highlight that the results from the present study reflect the number of individuals per event based on visible individuals within the scope of the camera, and that the real group sizes may actually be larger.No adult males were present in the mixed herd in two cases within the present study; however, there were always 2YM and a few unidentified individuals, suggesting that the herd should not be considered as a pure “nursery herd”, as known for sexually dimorphic antelope species50.Calves are born in the NKNP during the period comparable to that of Bandia, Fathala and Chinko, i.e. during the early dry season16. The higher proportion of calves in the dry season corresponds with the nursing period of six months for WDE44. Given a pregnancy length of nine months, the WDE mating season in NKNP peaks in January/February, which also corresponds with the formation of large herds with multiple males, as similarly seen in Chinko and Cameroon13,14.Vital ratesThe sex ratio of the WDE in the NKNP was female-biased. The skewed adult sex ratio reflects the lower survival rate of males in comparison with females, typical for polygynous species51. This result also corresponds with the findings from other Derby eland populations, namely from Chinko, where the bias towards females in the adult sex ratio was even more pronounced (0.67:113). A similar ratio was found in the hunting reserves within Cameroon35, but also in the semi-captive population, without hunting and without predators34. As the ratio in NKNP was less skewed than that within Chinko and Cameroon, a lower or zero selectivity for males by hunters/poachers is expected.The population of WDE in the NKNP showed a lower proportion of adults versus other age categories compared to the demographic structure of the WDE in the semi-captive breeding facilities of the Bandia and Fathala reserve33,52, and to those of the Eastern subspecies of Derby eland in the Central African Republic13 (see Table 3). The data from the present study also showed a surprisingly high breeding rate (likely close to 100%), as well as a high survival rate of yearlings. This combination of demographic characteristics should be highly favourable, and likely to lead to a significant population growth rate; however, this does not seem to be the case of the WDE population in the NKNP (please refer to further discussion about population size).In this context, the population of WDE in the NKNP was explored deeper, to examine possible scenarios of changes within the population structure. The changes in vital rates between two years of monitoring (2017 and 2018) were examined, by taking advantage of the possible recognition of the age category until two years of age, and the knowledge of the life tables of the enclosed, non-predated WDE population in the Bandia reserve34. Life tables were created for each year, and for males (M) and females (F) separately, according to the standard structure2, and based on two scenarios: a) only the observed number of JUV and 1Y (nx), and modelled 2Y (model ‘JUV + 1Y’); b) the observed number of JUV and 2Y (nx) (model ‘JUV + 2Y’). Then, estimations of animals in age categories based on two parameters were calculated: (i) based on the mortality rate (qx) known from the Bandia reserve (Senegal), and (ii) based on the recorded number of animals (NAD), to calculate the estimation of mortality rate (for details, see Additional file 1: Table S2).The resulting values demonstrated that with survival rates comparable to a population without predation and poaching, the number of adults would be twice or three times higher than currently detected in the present study. Yet, considering the recorded number of adult individuals, the annual adult survival rate was considerably low, i.e. 59–69% in males and 67–82% for females. To conclude, the demographic structure of WDE in NKNP showed a high breeding rate, moderate juvenile survival, high survival rate of yearlings, and a low survival rate of adults.Juvenile survival is one of the most fluctuating vital rate parameters, sensitive to population density, stochastic environmental variation, and predation53,54,55. Given the high proportion of juveniles within the population, and the breeding rate higher than that in Cameroon (74%14) and within the captive population (77%34), the juvenile survival rate does not seem to negatively affect the population growth in the NKNP. High breeding rates could be a more robust determinant of population change than AD mortality53, and it is therefore possible that the WDE population size is stable in the NKNP, or even increasing, despite the low adult survival rates. On the other hand, the relatively low numbers of AD individuals in the population indicates low survival rates, which may lead to the decline and final crash of the population54. It is acknowledged that data from two consecutive years was used in the present study, which were not comparable due to different CT settings, and that long-term monitoring, which accounts for variability in vital rates, would be a conservation essential to identify the trend and population change.Based on the present findings of WDE spatiotemporal behaviour and estimates of vital rates, several explanations about multiple processes interacting in the environmental, anthropogenic and conservation context of the park, which inherently affect the small population of WDE, can be inferred. One explanation may suggest that a low proportion of AD WDE and higher JUV survival rates may reflect the influence of growing populations of apex predators in the NKNP, specifically the population of lions56, which may preferentially target the adult individuals57. The age-sex structure also encourages the interpretation that the adult animals are exposed to human-related factors, which prevents them from expanding from the core area of their distribution, exacerbating male-male competition in the limited space34. The poaching activity was also highlighted as an existing threat to WDE populations35. However, law enforcement has been substantially intensified in the core and south-eastern part of the NKNP since 201758, and lion-conservation actions are specifically supported. Thus, the predator populations may have started to grow, which is confirmed by the relative high trapping rate of lions in this core area8. Hence, increased predation may interfere with other environmental factors and consequently affect the WDE population dynamics at the level of AD individuals55,59.A complementary scenario may highlight other factors, specifically, those which maintain the WDE population within a certain spatial extent of the park, i.e. Mont Assirik and Mansa Fara marsh zone. This area can be delimited either ecologically by specific unidentified resources, or by anthropogenic factors, namely a highly frequented trade road crossing the park, wild bushfires, and intensive livestock encroachment in a large band from the borders of the park, inwards (up to 10 km). There is also a vast area in the central part of the park that offers an important space with a supposed carrying capacity for large herbivore populations. This area is, however, outside of the zone of intensified law enforcement, and suffers from inadequate surveillance in the long-term, due to the absence of tracks and therefore being difficult for rangers to access. This area certainly represents an attractive zone for targeted illegal hunting actions. These limiting factors constrain large mammals to concentrate within the zone of Mount Assirik and Mansa Fara marsh, which, in turn, makes animal populations vulnerable to any potential environmental or man-induced incidents, like bush fire.Population sizeThe estimated population size of 195 individuals corresponds with the range of most recent estimates of the WDE population size in the NKNP, i.e. 100–200 (approximately 170) individuals6,7,60. Given the fact that the model contains only the data for AD animals (as no other age category had recapture records), it may be considered that this estimate refers to the number of adult individuals in the population. With regards to Table 3, showing that adults are likely to form 43 to 44% of the whole population, it may be inferred that the actual number of WDE in the NKNP could be higher, even up to 300 individuals, if the data are corrected for the 22% of unidentified individuals. The WDE density estimate of 0.138 individuals/km2 was comparable to densities of Eastern Derby eland in CAR (densities ranging between 0.04 and 0.16 individuals/km2), in Chinko13, and ranging between 0.002 and 0.1 individuals/km2 in the northern CAR61, as well as in Cameroon, with densities ranging between 0.002 and 0.08 individuals/km262. On the other hand, in comparison to other antelope species, the estimated WDE density falls within the range of densities of large herbivores reported from many other sites in African protected areas63, where lower values correspond to the larger areas and are also associated with large browsers, i.e. to the type of diet. Maximum densities of a healthy undisturbed DE population were estimated at about 0.5 individuals/km249, and can reach up to 1.19 individuals/km2 in intensively surveyed hunting zones in Northern CAR61. Thus, the density of WDE in the NKNP could be potentially higher. More

From the phytosociological relevés performed in each sampling area it is evident that hives were positioned in grasslands rich in Alpine herbaceous species (Table S1). In fact, among the 169 identified species, 85% were herbaceous species common in meadows (of Arrhenatherion elatioris and Triseto flavescentis-Polygonion bistortae phytosociological alliance) and acidophilus pastures (Siversio-Nardetum). 15% of the species were trees and shrubs (not abundant in the floristic relevés of the apiary areas considered), including some of beekeeping interest such as: Rhododendron ferrugineum, Castanea sativa and Rubus idaeus. From the MDS biplot (Fig. 2) elevation is the main ecological variable that differentiates sampling areas. In particular, the relevés of stations B and F are characterized by a floristic composition which is different from the areas at higher elevation (characterized by a higher presence of microthermal alpine species). This is due to the separation between the sub-montane belt and the high mountain belt vegetation on the 1.300 m a.s.l. line in the study area25.Figure 2MDS of the phytosociological relevés. Capital letters indicate the six sampling areas, the 1.300 m a.s.l. contour line that separates sub-montane belt and high mountain belt vegetation is highlighted in red.Full size imageAlthough the beehives were positioned in mountain grasslands, melissopalynological analysis presented a different picture. The pollen of numerous species detected through the floristic relevés were found in the honey samples via melissopalynological analysis, although the latter did not totally overlap with the floristic characterization of the area, in particular from a “quantitative” point of view. In fact, the floristic relevés showed a relative richness of herbaceous species (Table S1) peculiar of mountain grasslands that would seem promising for the production of wildflower honeys. Conversely, in the melissopalynological analysis the species considered interesting but not predominant in the botanical description were relevant (Fig. 3 and Table S2).Figure 3MDS of the melissopalynological analysis of the six samples (dots) of mountain wildflower honeys produced in the stations considered. The crosses are the pollens found in the honey samples, the most important are indicated.Full size imageThe premises to produce wildflower honey is that the botanical species contributing must be different and sometimes very numerous, without any of them assuming a dominant character. However, this was not fully evident in our research: although it was possible to identify more than seventy species through melissopalynological analysis and even more through the floristic characterization of the areas, most of them were defined as minor or sporadic pollen (Table S2). Even though apiaries were in mountain grasslands, the most relevant role was played by some woody species/shrubs: Rubus (presumably Rubus idaeous L., identified in the floristic relevés) and rhododendron (Rhododendron ferrugineum L.) for the mountain/subalpine belt and Castanea and Ericaceae (heather) in the submountain belt. Following the rules to define ‘‘unifloral honey’’, three of the wildflower honeys could be defined unifloral or bifloral:

Rhododendron unifloral: honey A (Rhododendron 47.18%), and honey C (Rhododendron 62.93%);

Raspberry unifloral: honey B (Rubus 67.12%)

Raspberry and Rhododendron bifloral: honey D (Rhododendron 34.27% and Rubus 34.74%) as well as honey E (Rubus 44.25%, Rhododendron 34.14%).

Honey F, due to the contribution of pollen from Tilia genus (that was detected only in this sample as an important sporadic pollen, 3.5%) Castanea (96.4% in honey F, but it should be noted that chestnut pollen is an overrepresented pollen) and in the second count Ericaceae (32.45%, that was considered a secondary pollen together with Rubus, with a percentage of 38.59% in honey F) differed from the other honeys (Fig. 3).Rubus pollen was anyway present in good amounts in all the samples considered, and was a dominant pollen in honey B, a secondary pollen in honeys C, D, E and F and a minor pollen in honey A. Sorbus and Tilia pollens were detected only in honey F, while no rhododendron was detected in honey F. Honey D was characterized by a percentage higher than the “rare pollen” category of some important alpine essences, such as Liliaceae, Centaurea, Campanulaceae, Anthyllis f., Polygonum bistorta, Lotus alpinus and Potentilla/fragaria (Table S2).Although wildflower honeys are intrinsically characterized by a high variability compared with unifloral honey, this shows the importance of the formal characterization of honey to obtain a product which satisfies consumer expectations, and it was demonstrated that the botanical origin of honey cannot be based on the claims of local beekeepers by considering the predominant flowers surrounding the hive.Although honeybees are considered supergeneralists in their foraging choices, there are certain key species or plant groups that are particularly important in honeybee foraging2, and many were identified in the botanical characterization of the area, including Rubus idaeus L., Calluna vulgaris L., rhododendron and some present in the broad-leaved woods mentioned such as chestnut (Castanea sativa Mill.) or plants of Tilia genus. In the research work by Hawkins et al.2, Rubus fruticosus L. was among the frequently found species and tree pollen belonging to Castanea sativa L. as well as, for example, species of Malus, Salix and Quercus spp, was frequently seen. These kinds of preferences could relate to the ease of availability and abundance of the plant, the quality and abundance of the nectar and pollen and/or specific nutrients or trace elements provided by these species or neurological aspects (as will be discussed further). As referred by beekeepers, over the last decades the production of mountain wildflower honey, that often does not meet the characteristics expected and presents flavours that are reminiscent of other kinds of honey such as rhododendron or linden or chestnut, is becoming more and more critical and this was absolutely confirmed by this study.This could be linked to the fragmentation of an important habitat of the Alps—mountain grasslands (meaning pastures and meadows) for anthropic and climatic reasons8,9. Honeybees from the same colony forage across areas spanning up to several hundred square kilometres, and at linear distances as far as 9 km from the hive41. Onlooker bees are those in charge of finding nectar sources and of giving instructions to the employed bees, the other foraging bees, that communicate the necessity to look for new resources of food to the onlookers through continuous dance communication42. Among the onlookers, there is a difference between the bees that scout for different nectar sources or recruit to well known floral resources43 and there is an optimal ratio of scouts to recruits, for the most effective collective foraging41. However, this balance may change based on the structure of the landscape in which the bees forage for food44,45,46. Theoretical models47,48 and empirical tests49 suggest that when resources are concentrated into a small number of highly rewarding patches, colonies perform best with few scouts and many recruits, while when resource patches are small, evenly distributed, and easy to locate, successful colonies invest more in scouting than in recruitment. This is strictly linked to climate and social changes in the mountains: mountain grasslands are no longer evenly distributed and easily localizable, as they are scattered among expanding areas of shrublands and forests9 and, for the above-mentioned reasons, it is more efficient for the colony to invest in more recruiters than scouters, as recruiters will identify a small number of highly rewarding patches, such as raspberry or rhododendron shrublands or linden and chestnut woods, that are highly rewarding and very different in quality.This overlaps with individual and collective honeybee behaviour driven by proximate physiological mechanisms that involve the tryptophan metabolism via kynurenine pathway that is one of main neuroprotective mechanisms. In this research, many of the differences/similarities among the samples might be attributed to metabolic alterations within this pathway, represented by relative amounts of kynurenic acid. However, different quinoline structures have also been identified (Fig. 4). Neurotransmitters play a central role in several of the biological processes that honeybees require to perform activities such as foraging behaviour50. A considerable amount of literature highlights the involvement of the neuroprotective kynurenine pathway (KP) final product kynurenic acid (KinA) in the regulation of the stress-related hormone dopamine in the honeybee as well as in other animal species51,52. The major known source of dietary KynA are pollen and nectar produced by sweet chestnuts53,54 and it has been verified that this compound is found in high concentrations in chestnut flowers55. This is coherent with the results of this study: chestnut pollen was found in honey F, produced in the lower station where chestnuts also appear in the floristic relevés, and KynA was found to be a dominant compound in honey F. Interestingly, chestnut pollen was found as sporadic pollen in all the other samples, even those produced in the highest apiary stations (Table S2).Figure 4Kinurenic acid and 3-hydroxyquinaldic acid structure and content in the six honey samples, performed in triplicate. The box diagram representing the median with distribution interval between 25 and 75%.Full size imageFurther, KinA may possess positive properties in a number of pathologies of the gastrointestinal tract, especially colitis, colon obstruction or ulceration56,57. It has been proposed that KinA may also possess antioxidative properties56,57,58,59. This was confirmed by this study, since the wildflower honey with a high component of chestnut pollen was the one with the highest antioxidant properties at the FRSA test (66.61 ± 4.77%), even if lower than manuka honey (84.21 ± 1.04%), a dark honey that is a well-known nutraceutical product and has recently attracted attention for its biological properties, especially for its antioxidant and anti-microbial capacities60. Honey A showed the lowest power (22.40 ± 0.28%) while the other honeys ranked around 40% (Fig. 5). Interestingly, metabolomic analysis revealed the presence of 3-hydroxyquinaldic acid (Fig. 4), which is a kynurenic acid isomer and, although its function has not been elucidated in detail, a few literature data indicate its role as a precursor of naturally occurring peptide antibiotics from the quinomycin family61.Figure 5Results of the FRSA test. Capital letters represent the six honey samples considered. Manuka honey was used as a control.Full size imageIn order to evaluate the ability of honey to induce wound closure, a scratch wound assay was performed (Fig. 6)62. Scratch assay creates a gap in confluent keratinocyte monolayer to mimic a wound. It has already been demonstrated that honeys are able to induce wound closure63 to different extents depending on honey origins and properties.Figure 6The scratch wound test in keratinocytes, HaCaT cells, exposed to honeys. (a) The digitalized pictures of scratched cells after 24 h exposure to 0.5% (w/v) of honeys. (b) The closing percentage wound values after 24 h exposure. Statistics on bars indicate differences compared to the control (CTRL) condition determined by a One-Way ANOVA followed by Dunnett’s test (****p  More

Soil and litter samplingMineral soil (0–15 cm) was collected at the Elizabeth Woods site, a 120-year-old deciduous forest in West Virginia, US (39° 32′ 50.6″ N, − 80° 00′ 00.4″ W). Soils were collected from four 20 × 20 m plots dominated by either AM-associated trees (i.e. Liriodendron tulipifera and Acer saccharum), or ECM-associated trees (i.e. Quercus rubra, Quercus velutina and Carya ovata). These sites have been characterized previously as Culleoka-Westmoreland silt loam soils at the AM sites and Dormont and Guernsey silt loams at the ECM sites40. Soils were also characterized by C:N ratios 11.7 and 14.1 for the AM and ECM soils respectively, with a pH of 6.8 for both soils. Soils with the same mycorrhizal status were pooled and homogenized, air-dried at room temperature for ~ 24 h and sieved through 2.0 mm mesh before the initiation of the experiment. Uniformly 13C labeled litter ( > 97 atom % 13C) from Quercus robur (i.e., ECM substrate) and Liriodendron tulipifera (i.e. AM substrate) leaves (Isolife BV, Wageningen, NL) were incubated in soil mesocosms in a factorial design with five replicates for each treatment combination (2 soil types × 2 substrate types), along with five replicate controls (no 13C substrate addition) for each soil type. The 13C enriched substrates were dried and ground to a powder and added in a suspension of 0.5 ml sterile water to 20 g of soil at a concentration of 400 ug 13C g−1 soil. The control soils received 0.5 ml sterile water additions. These incubations were well mixed and kept at 60% water-holding capacity for the 21-day period at room-temperature18. Chemical characteristics of soils and plant substrates are provided in Table S1.DNA processing and qSIPFor quantitative stable isotope probing, DNA was extracted, quantified, ultracentrifuged, fractionated and sequenced as described in18,26. DNA was extracted using a MoBio PowerSoil HTP Kit following the manufacturer’s instructions. For stable isotope probing, 5 ug of DNA was loaded into a 5-ml ultracentrifuge tube with ~ 3.5 ml of a saturated cesium chloride (CsCl) solution and ~ 900 ml gradient buffer (200 mM Tris, 200 mM KCl, 2 mM EDTA). DNA was separated via ultracentrifugation at 127,000g for 72 h using a TLN-100 rotor in an Optima Max bench top ultracentrifuge (Beckman Coulter, Fullerton, CA, USA). Tubes were fractionated into ~ 25 fractions of 150 µl each, and the density of each fraction was measured with a Raichart AR200 digital refractometer. DNA was purified using an isopropanol precipitation method. The 16S rRNA gene was subsequently quantified and sequenced in samples containing DNA, within the density range 1.660–1.735 gml−1 (~ 10 fractions per sample). To quantify the 16S rRNA gene, quantitative PCR was performed in triplicate using a QuantStudio 5 applied biosystems (Thermo Fisher Scientific) and primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACVSGGGTATCTAAT-3′)41. The PCR program used was as follows: 95 °C for 2 min followed by 45 cycles of 95 °C for 30 s, 64.5 °C for 30 s and 72 °C for 1 min. Libraries were sequenced on an Illumina MiSeq instrument (Illumina, Inc., San Diego, CA, USA) using a 300-cycle v2 reagent kit. Fungal 18S rRNA gene copies in each fraction were also quantified using primers 1380F (5′-CCCTGCCHTTTGTACACAC-3′) and 1510R (5′-CCTTCYGCAGGTTCACCTAC-3′). The PCR program used was as follows: 98 °C for 3 min followed by 40 cycles of 98 °C for 45 s, 60 °C for 45 s and 72 °C for 30 s. DNA fractions were amplified for fungal ITS rRNA genes using primers ITS4F (5′-AGCCTCCGCTTATTGATATGCTTAART-3′) and 5.8SF (5′-AACTTTYRRCAAYGGATCWCT-3′)42 and 300-bp paired-end read chemistry on an IlluminaMiSeq (Illumina, Inc., San Diego, CA, USA). The PCR program used was as follows: 95 °C for 6 min followed by 35 cycles of 95 °C for 15 s, 55 °C for 30 s, and 72 °C for 1 min. DNA fractions were then sequenced using a 500 cycle v2 reagent kit.Files came pre-split and joined multiple paired ends that we combined to pick operational taxonomic units (OTU). Open reference OTUs were picked at 97% identity using SILVA 128 release database for Bacteria and RDP database for Fungi. Taxa were analyzed at the ‘OTU’ level from the QIIME L7 table. Calculation of 13C excess atom fraction (EAF) was performed for each taxon as described previously18,19. Briefly, using the CsCl density gradient data, a weighted average density (WAD) was computed for each taxon’s DNA extracted from control soils that did not receive an isotopically enriched substrate. This natural abundance WAD was then compared to the taxon’s WAD following incubation with the 13C enriched material. The change in WAD can be used to quantify the amount of isotope incorporated into the DNA17,18. Preliminary data analysis revealed an effect of ultracentrifuge tube on estimation of phylotype weighted average density, probably a consequence of slight differences in CsCl density gradients between tubes. This technical error was corrected as previously described18,19. In addition to the samples subjected to qSIP analysis we also extracted and analyzed fungal and bacterial OTU’s from control soils where the DNA was extracted prior to incubation.FTICR-MS and lipidomic analysesSoil from substrate-incubated and controls mesocosms were processed and analyzed with Fourier transform ion cyclotron resonance mass spectrometry (FTICR-MS), using a 12 T Bruker SolariX FTICR mass spectrometer at the Environmental Molecular Sciences Laboratory in Richland, WA, as described in Fudyma et al.43. Briefly, 100 mg of dried soil or litter substrate was extracted using an adjusted Folch extraction44. Extraction was performed on each sample by sequentially adding 2 ml MeOH, followed by a 5 s vortex; 4 ml CHCl3, followed by a 5 s vortex; sonication at 25 °C for 1 h (CPX3800 Ultrasonic Bath, Fisherbrand); addition of 1.25 ml of H2O, followed by a slight mix to achieve bi‐layer separation; and incubated at 4 °C overnight. The top, aqueous layer (metabolite—polar) was pipetted off into 1 ml glass vials and stored at − 80 °C until FTICR‐MS. The bottom, chloroform layer was dried down and stored in 50:50 methanol:chloroform until lipidomics analysis.A standard Bruker electrospray ionization (ESI) source was used to generate negatively charged molecular ions in the metabolite fraction. Samples were then introduced directly to the ESI source. The instrument settings were optimized by tuning on a Suwannee River fulvic acid (SRFA) standard, purchased from International Humic Substances Society (IHCC). Blanks (HPLC grade methanol) were analyzed at the beginning and end of the day to monitor potential carry over from one sample to another. The instrument was flushed between samples using a mixture of water and methanol. One hundred and forty‐four individual scans were averaged for each sample and internally calibrated using an organic matter homologous series separated by 14 Da (CH2 groups). The mass measurement accuracy was less than 1 ppm for singly charged ions across a broad m/z range (m/z 300– 800). Data analysis software (Bruker Daltonik version 4.2) was used to convert raw spectra to a list of m/z values, applying the FTMS peak picker module with a signal-to noise ratio (S/N) threshold set to 7 and absolute intensity threshold set to the default value of 100. Chemical formulae were then assigned using in-house software following the compound identification algorithm that was described in Tolić et al.45. Peaks below 200 and above 800 were dropped to select only for calibrated and assigned peaks. Chemical formulae were assigned based on the following criteria: S/N  > 7 and mass measurement error  800 were not detected in our samples. The m/z values represent the molecular mass (in Dalton) of the detected ions since all detected ions were singly charged ions. While our results do not represent a quantitative characterization of OM, the values presented are relative differences and should be representative of the samples. Finally, we would like to acknowledge that we were not able to see any clear evidence of 13C label in our FTICR-MS analysis of the soil samples. The lack of 13C label in our FTICR-MS analysis of the soil samples even though they received labeled substrate could be either due to the fact that most of the labeled substrates produced by microbial activities were of low molecular weight, which cannot be detected by FTICR-MS and/or the leftover labeled substrate was of low abundance compared to the organic compounds previously present in the soil matrix. As such, we used the FTCIR-MS data to identify shifts in the overall composition of the chemical compounds in each soil.Lipids in the chloroform fraction were analyzed by LC‐MS/MS in both positive and negative ESI modes using a linear trap quadropole (LTQ) Orbitrap Velos mass spectrometer (Thermo Fisher Scientific), as described in detail previously46. Lipid species were identified using the LIQUID tool46 followed by manual data inspection. Confidently identified lipid species were quantified using MZmine47 and the peak intensities were normalized by linear regression and central tendency (i.e., identifying a central or typical value for a probability distribution) using InfernoRDN.Statistical analysisAll data analyses were performed using R 3.2.048. To examine the effects of soil type, substrate type and their interaction in the bacterial, fungal and chemical composition of DOM and the lipid pool; Bray–Curtis distance matrices were compared with permutational multivariate analysis of variance (PerMANOVA) and visualized with Principle Coordinate Analysis (PCoA) using vegan package49. PerMANOVA analysis were run on the relative abundance and on the 13C EAF of individual microbial taxa, separately for both bacterial and fungal communities.The analyses for FTICR-MS were performed separately for control and incubated soils using all assigned molecular formulae remaining after quality filtering31. In all cases, we applied a Z-score standardization before calculating Bray–Curtis distance matrices49. We analyzed the results from FTICR-MS as resulting from the decomposition of the added substrates for two reasons. First, this is a fully factorial design where individual soil samples were split to either receive AM poplar or ECM oak litter substrate. Thus, each soil sample starts with the same characteristics and the changes at the end of the incubation period should reflect the processing of litter. Second, we excluded molecular formulae present in the litters and thus, the differences we report in each soil type are derived from this processing (or the lack of it).We calculated aggregated indices that characterize both the composition and the physicochemical properties of the microbial (both bacteria and fungi) and the SOM and lipid pool34,36. For bacterial and fungal communities, we quantified Shannon–Weaver diversity index for each sample H′ = (-{sum }_{i=1}^{S} pi ln(pi)) (where pi is the proportion of species I) using the relative abundance of individual microbial taxa50. To find the percent of substrate assimilation by individual taxa, we calculated the proportion of C assimilated by each group as previously described18,51 as a percent. For SOM and lipid molecular formulae, we separately calculated weighted means of formula-based characteristics (i.e. m/z, Aromaticity Index—AI; H/C, O/C, and Nominal Oxidation State of Carbon-NOSC) as the sum of the product of the single-formula information (i.e. m/zi, AIi, H/Ci and NOSCi) and the relative intensity (Ii) divided by the sum of all intensities (e.g., m/z sample1 = ({sum }_{i=1}^{S})(m/zi ·Ii)/Σ(Ii)). With these metrics we obtained sample-level information related to the molecular size (i.e. m/z), the molecular bioavailability (i.e. higher H/C ratio), the molecular reactiveness (i.e. lower AI) and the energetic rewards from molecular oxidative degradation (i.e. higher NOSC) of the SOM, which allows to infer the potential of decomposition products to form stable SOM12,31,35. Detailed information of the calculated indices can be found in the literature31,35,36.We further tested the effects of soil type, substrate type and their interaction on each index using the “lm” function from the “stats” package. In these analyses, P values were approximated by an F test using Type II ANOVA tests with Kenward-Roger Degrees of Freedom52. When interactions between soil and substrate type were found at P  More

1.Stanfield, B. J., Bliss, J. C. & Spies, T. A. Land ownership and landscape structure: A spatial analysis of sixty-six Oregon (USA) Coast Range watersheds. Landsc. Ecol. 17, 685–697 (2002).Article 

Google Scholar 
2.Spies, T. et al. Using an agent-based model to examine forest management outcomes in a fire-prone landscape in Oregon, USA. Ecol Soc 22, 25. https://doi.org/10.5751/ES-08841-220125 (2017).Article 

Google Scholar 
3.Zald, H. & Dunn, C. J. Severe fire weather and intensive forest management increase fire severity in a multi-ownership landscape. Ecol. Appl. 28, 1068–1080 (2018).Article 

Google Scholar 
4.Ager, A. A. et al. Network analysis of wildfire transmission and implications for risk governance. PLoS ONE 12, e0172867. https://doi.org/10.1371/journal.pone.0172867 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
5.Abatzoglou, J. & Williams, A. P. Impact of anthropogenic climate change on wildfire across western US forests. Proc. Natl. Acad. Sci. 113, 11770–11775 (2016).ADS 
CAS 
Article 

Google Scholar 
6.Sheehan, T., Let, D. B. & Ferschweiler, K. Projected major fire and vegetation changes in the Pacific Northwest of the conterminous United States under selected CMIP5 climate futures. Ecol. Model. 317, 16–29 (2015).Article 

Google Scholar 
7.Spies, T. A. et al. Examining fire-prone forest landscapes as coupled human and natural systems. Ecol. Soc. 19, 9. https://doi.org/10.5751/ES-06584-190309 (2014).Article 

Google Scholar 
8.Watkins, T. H. Untrammeled by man: The making of the Wilderness Act of 1964. Audubon 91, 74–90 (1989).
Google Scholar 
9.Huffman, D. W., Roccaforte, J. P., Springer, J. D. & Crouse, J. E. Restoration applications of resource objective wildfires in western US forests: A status of knowledge review. Fire Ecol. 16, 18. https://doi.org/10.1186/s42408-020-00077-x (2020).Article 

Google Scholar 
10.Charnley, S., Spies, T. A., Barros, A. M. G., White, E. M. & Olsen, K. A. Diversity in forest management to reduce wildfire losses: Implications for resilience. Ecol. Soc. 22, 1. https://doi.org/10.5751/ES-08753-220122 (2017).Article 

Google Scholar 
11.Lake, F. K. & Long, J. W. Fire and tribal cultural resources. Report No. PSW-GTR-274, (USDA USFS Pacific Southwest Research Station, Albany, CA, 2014).12.Binkley, C. S., Aronow, M. E., Washburn, C. L. & New, D. Global perspectives on intensively managed plantations: Implications for the Pacific Northwest. J. For. 103, 61–64 (2005).
Google Scholar 
13.Palaiologou, P. et al. Fine-scale assessment of cross-boundary wildfire events in the western United States. Nat. Hazards Earth Syst. Sci. 19, 1755–1777. https://doi.org/10.5194/nhess-19-1755-2019 (2019).ADS 
Article 

Google Scholar 
14.Ager, A. A., Palaiologou, P., Evers, C., Day, M. A. & Barros, A. M. Assessment of wildfire transmission from national forests to communities in the Western United States. 52 (USDA Forest Service, 2017).15.Steelman, T. U. S. wildfire governance as a social-ecological problem. Ecol. Soc. 21, 3. https://doi.org/10.5751/ES-08681-210403 (2016).Article 

Google Scholar 
16.Charnley, S., Kelly, E. C. & Fischer, A. P. Fostering collective action to reduce wildfire risk across property boundaries in the American West. Environ. Res. Lett. 15, 025007 (2020).ADS 
Article 

Google Scholar 
17.USDA Forest Service. Towards shared stewardship across landscapes: An outcome-based investment strategy. Report No. FS-118, (USDA Forest Service, Washington, DC, 2018).18.USDA Forest Service. National Cohesive Wildland Fire Management Strategy. http://www.forestsandrangelands.gov/strategy/index.shtml (2015).19.Marsik, M. et al. Regional-scale management maps for forested areas of the Southeastern United States and the US Pacific Northwest. Sci. Data 5, 1–13 (2018).Article 

Google Scholar 
20.Franklin, J. F. & Dyrness, C. T. in General Technical Report PNW-GTR-008 427 (U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station, Portland, OR, 1973).21.Simpson, M. Central Oregon Area Ecology and Forest Health Program (ed Pacific Northwest Region USDA Forest Service) (Bend, OR, 2013).22.MTBS. MTBS Data Access: Burned areas boundaries. https://www.mtbs.gov/index.php/direct-download. (2020).23.Picotte, J. J. et al. Changes to the monitoring trends in burn severity program mapping production procedures and data products. Fire Ecol. 16, 1–13 (2020).Article 

Google Scholar 
24.Meddens, A. J. H., Kolden, C. A., Lutz, J. A., Abatzoglou, J. & Hudak, A. T. Spatiotemporal patterns of unburned areas within fire perimeters in the northwestern United States from 1984 to 2014. Ecosphere 9, e02029 (2018).Article 

Google Scholar 
25.USGS. (USGS Gap Analysis Program (GAP), 2016).26.Gaines, L., Hemstrom, M., Kagan, J. & Salwasser, J. Integrated landscape assessment project final report. 62 (The Institute for Natural Resources, Oregon State University, Corvallis, Or, 2013).27.Bond, W. J. & Keeley, J. E. Fire as a global ‘herbivore’: The ecology and evolution of flammable ecosystems. Trends Ecol. Evol. 20, 387–394 (2005).Article 

Google Scholar 
28.Manly, B., McDonald, L. & Thomas, D. Resource Selection by Animals (Chapman & Hall, 1993).Book 

Google Scholar 
29.Bajocco, S., Pezzatti, G. B., Mazzoleni, S. & Ricotta, C. Wildfire seasonality and land use: When do wildfires prefer to burn?. Envrion. Monit. Assess. 164, 445–452 (2010).CAS 
Article 

Google Scholar 
30.Bajocco, S. & Ricotta, C. Evidence of selective burning in Sardinia (Italy): Which land cover classes do wildfires prefer?. Landsc. Ecol. 23, 241–248 (2008).Article 

Google Scholar 
31.Barros, A. M. G. & Pereira, J. M. C. Wildfire selectivity for land cover type: Does size matter?. PLoS ONE 9, e84760 (2014).ADS 
Article 

Google Scholar 
32.R Package ‘phuassess’ (2016).33.Fattorini, L., Pisani, C., Riga, F. & Zaccaroni, M. The R package “phuassess” for assessing habitat selection using permutation-based combination of sign tests. Mamm. Biol. 83, 64–70 (2017).Article 

Google Scholar 
34.Fattorini, L., Pisani, C., Riga, F. & Zaccaroni, M. A permutation-based combination of sign tests for assessing habitat selection. Environ. Ecol. Stat. 21, 161–187 (2013).MathSciNet 
Article 

Google Scholar 
35.R: A Language and Environment for Statistical Computing v.3.5.3 (R Foundation for Statistical Computing, Vienna, Austria, 2019).36.ArcGIS Desktop: Release 10 (Environmental Systems Research Institute, 2011).37.MATLAB Release 2019a v. 2019a (The Mathworks, Inc., 2019).38.Collins, B. & Stephens, S. Fire scarring patterns in Sierra Nevada wilderness areas burned by multiple wildland fire use fires. Fire Ecol. 3, 53–67 (2007).Article 

Google Scholar 
39.Reilly, M. J. et al. Cumulative effects of wildfires on forest dynamics in the eastern Cascade Mountains USA. Ecol. Appl. 28, 291–308 (2018).Article 

Google Scholar 
40.Johnston, J. D., Kilbride, J. B., Meigs, G. W., Dunn, C. J. & Kennedy, R. E. Does conserving roadless wildland increase wildfire activity in western US national forests?. Environ. Res. Lett. 16, 084040 (2021).ADS 
Article 

Google Scholar 
41.Schultz, C. A., Thompson, M. P. & McCaffrey, S. M. Forest service fire management and the elusiveness of change. Fire Ecol. 15, 1–15 (2019).Article 

Google Scholar 
42.Ager, A. A., Houtman, R., Day, M. A., Ringo, C. & Palaiologou, P. Tradeoffs between US national forest harvest targets and fuel management to reduce wildfire transmission to the wildland urban interface. For. Ecol. Manag. 434, 99–109 (2019).Article 

Google Scholar 
43.NWCG. Guidance for Implementation of Federal Wildland Fire Management Policy (2009).44.Franklin, J. F. et al. Extent and Distribution of Old Forest Conditions on Washington Department of Natural Resources-Managed Forest Lands in Eastern Washington (Washington Department of Natural Resources, 2007).45.Stephens, S. L. et al. Fire and climate change: Conserving seasonally dry forests is still possible. Front. Ecol. Environ. 18, 354–360 (2020).Article 

Google Scholar 
46.Long, J., Lake, F. K., Lynn, K. & Viles, C. Tribal ecocultural resources and engagement. Report No. General Technical Report PNW-GTR-966, 851-917 (USDA – USFS, 2018).47.Scott, J. H. & Burgan, R. E. Standard fire Behavior Fuel Models: A Comprehensive Set for Use with Rothermel’s Surface Fire Spread Model. Report No. RMRS-GTR-153, 72 (USDA Forest Service, Rocky Mountain Research Station, 2005).48.Fernandes, P. M., Pacheco, A. P., Almeida, R. & Claro, J. The role of fire-suppression force in limiting the spread of extremely large forest fires in Portugal. Eur. J. For. Res. 135, 253–262 (2016).Article 

Google Scholar 
49.WADNR, W. D. o. N. R. Forest Health Assessment and Treatment Framework (RCW 76.06.200) (Washington State Department of Natural Resources, 2020).50.Collins, B. M. & Stephens, S. L. Managing natural wildfires in Sierra Nevada wilderness areas. Front. Ecol. Environ. 5, 523–527 (2007).Article 

Google Scholar 
51.Holden, Z. A., Morgan, P., Rollins, M. G. & Kavanagh, K. Effects of multiple wildland fires on ponderosa pine stand structure in two southwestern wilderness areas, USA. Fire Ecol. 3, 18–33 (2007).Article 

Google Scholar 
52.Hunter, M. E., Iniguez, J. M. & Farris, C. A. (U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station, 2014). More

1.Mittelbach, G. G. et al. Evolution and the latitudinal diversity gradient: speciation, extinction and biogeography. Ecol. Lett. 10, 315–331 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
2.Mannion, P. D., Upchurch, P., Benson, R. B. J. & Goswami, A. The latitudinal biodiversity gradient through deep time. Trends Ecol. Evol. 29, 42–50 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
3.Kinlock, N. L. et al. Explaining global variation in the latitudinal diversity gradient: meta‐analysis confirms known patterns and uncovers new ones. Glob. Ecol. Biogeogr. 27, 125–141 (2018).Article 

Google Scholar 
4.Wiens, J. J., Graham, C. H., Moen, D. S., Smith, S. A. & Reeder, T. W. Evolutionary and ecological causes of the latitudinal diversity gradient in hylid frogs: treefrog trees unearth the roots of high tropical diversity. Am. Nat. 168, 579–596 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
5.Wiens, J. J., Sukumaran, J., Pyron, R. A. & Brown, R. M. Evolutionary and biogeographic origins of high tropical diversity in old world frogs (Ranidae). Evolution 63, 1217–1231 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
6.Jansson, R., Rodríguez-Castañeda, G. & Harding, L. E. What can multiple phylogenies say about the latitudinal diversity gradient? A new look at the tropical conservatism, out of the tropics, and diversification rate hypotheses: phylogenies and the latitudinal diversity gradient. Evolution 67, 1741–1755 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
7.Economo, E. P., Narula, N., Friedman, N. R., Weiser, M. D. & Guénard, B. Macroecology and macroevolution of the latitudinal diversity gradient in ants. Nat. Commun. 9, 1778 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
8.Stephens, P. R. & Wiens, J. J. Explaining species richness from continents to communities: the time‐for‐speciation effect in Emydid turtles. Am. Nat. 161, 112–128 (2003).PubMed 
Article 
PubMed Central 

Google Scholar 
9.Wiens, J. J. & Donoghue, M. J. Historical biogeography, ecology and species richness. Trends Ecol. Evol. 19, 639–644 (2004).PubMed 
Article 
PubMed Central 

Google Scholar 
10.Morley, R. J. Cretaceous and Tertiary climate change and the past distribution of megathermal rainforests. In Tropical Rainforest Responses to Climatic Change 1–31 (Springer Berlin Heidelberg, 2007). https://doi.org/10.1007/978-3-540-48842-2_1.11.Jablonski, D., Roy, K. & Valentine, J. W. Out of the tropics: evolutionary dynamics of the latitudinal diversity gradient. Science 314, 102–106 (2006).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Condamine, F. L., Sperling, F. A. H., Wahlberg, N., Rasplus, J.-Y. & Kergoat, G. J. What causes latitudinal gradients in species diversity? Evolutionary processes and ecological constraints on swallowtail biodiversity: phylogeny and latitudinal diversity gradient. Ecol. Lett. 15, 267–277 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
13.Rolland, J., Condamine, F. L., Beeravolu, C. R., Jiguet, F. & Morlon, H. Dispersal is a major driver of the latitudinal diversity gradient of Carnivora: dispersal and the latitudinal gradient of Carnivora. Glob. Ecol. Biogeogr. 24, 1059–1071 (2015).Article 

Google Scholar 
14.Fischer, A. G. Latitudinal variations in organic diversity. Evolution 14, 64 (1960).Article 

Google Scholar 
15.Rolland, J., Condamine, F. L., Jiguet, F. & Morlon, H. Faster speciation and reduced extinction in the tropics contribute to the mammalian latitudinal diversity gradient. PLoS Biol. 12, e1001775 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
16.Rosenzweig, M. L. Species diversity in space and time. (Cambridge University Press, 1995). https://doi.org/10.1017/CBO9780511623387.17.Allen, A. P., Gillooly, J. F., Savage, V. M. & Brown, J. H. Kinetic effects of temperature on rates of genetic divergence and speciation. Proc. Natl Acad. Sci. 103, 9130–9135 (2006).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Schemske, D. W. Ecological and evolutionary perspectives on the origins of tropical diversity. In Foundations of tropical forest biology (eds Chazdon, R. & Whitmore, T.) 163–173 (University of Chicago Press, Chicago, IL, 2002).19.Janzen, D. H. Why mountain passes are higher in the tropics. Am. Nat. 101, 233–249 (1967).Article 

Google Scholar 
20.Wahlberg, N. et al. Nymphalid butterflies diversify following near demise at the Cretaceous/Tertiary boundary. Proc. R. Soc. B. 276, 4295–4302 (2009).PubMed 
PubMed Central 
Article 

Google Scholar 
21.Aduse-Poku, K. et al. Systematics and historical biogeography of the old world butterfly subtribe Mycalesina (Lepidoptera: Nymphalidae: Satyrinae). BMC Evol. Biol. 15, 167 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
22.Kozak, K. M. et al. Multilocus species trees show the recent adaptive radiation of the mimetic Heliconius butterflies. Syst. Biol. 64, 505–524 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
23.Chazot, N. et al. Renewed diversification following Miocene landscape turnover in a Neotropical butterfly radiation. Glob. Ecol. Biogeogr. 28, 1118–1132 (2019).Article 

Google Scholar 
24.Chazot, N. et al. Priors and posteriors in Bayesian timing of divergence analyses: the age of butterflies revisited. Syst. Biol. 68, 797–813 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
25.Wahlberg, N., Wheat, C. W. & Peña, C. Timing and patterns in the taxonomic diversification of Lepidoptera (Butterflies and Moths). PLoS ONE 8, e80875 (2013).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
26.Heikkilä, M., Kaila, L., Mutanen, M., Peña, C. & Wahlberg, N. Cretaceous origin and repeated tertiary diversification of the redefined butterflies. Proc. R. Soc. B. 279, 1093–1099 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
27.Condamine, F. L., Nabholz, B., Clamens, A.-L., Dupuis, J. R. & Sperling, F. A. H. Mitochondrial phylogenomics, the origin of swallowtail butterflies, and the impact of the number of clocks in Bayesian molecular dating: mito-phylogenomics of swallowtail butterflies. Syst. Entomol. 43, 460–480 (2018).Article 

Google Scholar 
28.Espeland, M. et al. A comprehensive and dated phylogenomic analysis of butterflies. Curr. Biol. 28, 770–778 (2018). e5.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
29.Ree, R. H. & Smith, S. A. Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. Syst. Biol. 57, 4–14 (2008).PubMed 
PubMed Central 
Article 

Google Scholar 
30.Crisp, M. & Cook, L. Do early branching lineages signify ancestral traits? Trends Ecol. Evol. 20, 122–128 (2005).PubMed 
Article 
PubMed Central 

Google Scholar 
31.Meseguer, A. S. & Condamine, F. L. Ancient tropical extinctions at high latitudes contributed to the latitudinal diversity gradient*. Evolution 74, 1966–1987 (2020).Article 

Google Scholar 
32.Ziegler, A. et al. Tracing the tropics across land and sea: Permian to present. Lethaia 36, 227–254 (2003).Article 

Google Scholar 
33.Meng, J. & McKenna, M. C. Faunal turnovers of Palaeogene mammals from the Mongolian Plateau. Nature 394, 364–367 (1998).ADS 
CAS 
Article 

Google Scholar 
34.Archibald, S. B., Bossert, W. H., Greenwood, D. R. & Farrell, B. D. Seasonality, the latitudinal gradient of diversity, and Eocene insects. Paleobiology 36, 374–398 (2010).Article 

Google Scholar 
35.Baker, W. J. & Couvreur, T. L. P. Global biogeography and diversification of palms sheds light on the evolution of tropical lineages. I. Historical biogeography. J. Biogeogr. 40, 274–285 (2013).Article 

Google Scholar 
36.Liu, Z. et al. Global cooling during the Eocene-Oligocene climate transition. Science 323, 1187–1190 (2009).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
37.Eldrett, J. S., Greenwood, D. R., Harding, I. C. & Huber, M. Increased seasonality through the Eocene to Oligocene transition in northern high latitudes. Nature 459, 969–973 (2009).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
38.Saupe, E. E. et al. Climatic shifts drove major contractions in avian latitudinal distributions throughout the Cenozoic. Proc. Natl Acad. Sci. USA 116, 12895–12900 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Hawkins, B. A. & DeVries, P. J. Tropical niche conservatism and the species richness gradient of North American butterflies. J. Biogeogr. 36, 1698–1711 (2009).Article 

Google Scholar 
40.Mayr, G. Two-phase extinction of “Southern Hemispheric” birds in the Cenozoic of Europe and the origin of the Neotropic avifauna. Palaeobio. Palaeoenv. 91, 325–333 (2011).Article 

Google Scholar 
41.Veizer, J. & Prokoph, A. Temperatures and oxygen isotopic composition of Phanerozoic oceans. Earth-Sci. Rev. 146, 92–104 (2015).ADS 
CAS 
Article 

Google Scholar 
42.Zhang, Z. et al. Aridification of the Sahara desert caused by Tethys Sea shrinkage during the Late Miocene. Nature 513, 401–404 (2014).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
43.Feakins, S. J. et al. Northeast African vegetation change over 12 m.y. Geology 41, 295–298 (2013).ADS 
Article 

Google Scholar 
44.Jacobs, B. F. Palaeobotanical studies from tropical Africa: relevance to the evolution of forest, woodland and savannah biomes. Philos. Trans. R. Soc. Lond. B 359, 1573–1583 (2004).Article 

Google Scholar 
45.Jaramillo, C. Cenozoic plant diversity in the Neotropics. Science 311, 1893–1896 (2006).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Stebbins, G. L. Flowering plants: evolution above the species level. (Harvard University Press, 1974). https://doi.org/10.4159/harvard.9780674864856.47.Wahlberg, N. & Wheat, C. W. Genomic outposts serve the phylogenomic pioneers: designing novel nuclear markers for genomic DNA extractions of Lepidoptera. Syst. Biol. 57, 231–242 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
48.Philippe, H. et al. Resolving difficult phylogenetic questions: why more sequences are not enough. PLoS Biol. 9, e1000602 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Nee, S. Birth-Death models in macroevolution. Annu. Rev. Ecol. Evol. Syst. 37, 1–17 (2006).Article 

Google Scholar 
50.Ricklefs, R. E. Estimating diversification rates from phylogenetic information. Trends Ecol. Evol. 22, 601–610 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
51.Crisp, M. D. & Cook, L. G. Explosive radiation or cryptic mass extinction? Interpreting signatures in molecular phylogenies. Evolution 63, 2257–2265 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
52.Lambert, A. & Stadler, T. Birth–death models and coalescent point processes: the shape and probability of reconstructed phylogenies. Theor. Popul. Biol. 90, 113–128 (2013).PubMed 
MATH 
Article 
PubMed Central 

Google Scholar 
53.Rabosky, D. L. Extinction rates should not be estimated from molecular phylogenies: estimating extinction from molecular phylogenies. Evolution 64, 1816–1824 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
54.Quental, T. B. & Marshall, C. R. Diversity dynamics: molecular phylogenies need the fossil record. Trends Ecol. Evol. 25, 434–441 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
55.Burin, G., Alencar, L. R. V., Chang, J., Alfaro, M. E. & Quental, T. B. How well can we estimate diversity dynamics for clades in diversity decline? Syst. Biol. 68, 47–62 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
56.Louca, S. & Pennell, M. W. Extant timetrees are consistent with a myriad of diversification histories. Nature 580, 502–505 (2020).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
57.Sohn, J.-C., Labandeira, C. C. & Davis, D. R. The fossil record and taphonomy of butterflies and moths (Insecta, Lepidoptera): implications for evolutionary diversity and divergence-time estimates. BMC Evol. Biol. 15, 12 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
58.de Jong, R. Fossil butterflies, calibration points and the molecular clock (Lepidoptera: Papilionoidea). Zootaxa 4270, 1 (2017).Article 

Google Scholar 
59.Edger, P. P. et al. The butterfly plant arms-race escalated by gene and genome duplications. Proc. Natl Acad. Sci. USA 112, 8362–8366 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.Allio, R. et al. Genome-wide macroevolutionary signatures of key innovations in butterflies colonizing new host plants. Nat. Commun. 12, 354 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.Stamatakis, A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
62.Drummond, A. J., Suchard, M. A., Xie, D. & Rambaut, A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29, 1969–1973 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
63.Lanfear, R., Frandsen, P. B., Wright, A. M., Senfeld, T. & Calcott, B. PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol. Biol. Evol. 34, 772–773. https://doi.org/10.1093/molbev/msw260 (2016).64.Smith, S. A. Taking into account phylogenetic and divergence-time uncertainty in a parametric biogeographical analysis of the Northern Hemisphere plant clade Caprifolieae. J. Biogeogr. 36, 2324–2337 (2009).Article 

Google Scholar 
65.Beeravolu Reddy, C. & Condamine, F. An extended maximum likelihood inference of geographic range evolution by dispersal, local extinction and cladogenesis. bioRxiv. https://doi.org/10.1101/038695 (2016).66.Rabosky, D. L. Automatic detection of key innovations, rate shifts, and diversity-dependence on phylogenetic trees. PLoS ONE 9, e89543 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
67.Rabosky, D. L. et al. BAMMtools: an R package for the analysis of evolutionary dynamics on phylogenetic trees. Methods Ecol. Evol. 5, 701–707 (2014).Article 

Google Scholar 
68.Rabosky, D. L., Mitchell, J. S. & Chang, J. Is BAMM flawed? theoretical and practical concerns in the analysis of multi-rate diversification models. Syst. Biol. 66, 477–498 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
69.Dudas, G. et al. Virus genomes reveal factors that spread and sustained the Ebola epidemic. Nature 544, 309–315 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
70.Morlon, H., Parsons, T. L. & Plotkin, J. B. Reconciling molecular phylogenies with the fossil record. Proc. Natl Acad. Sci. 108, 16327–16332 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.Morlon, H. et al. RPANDA: an R package for macroevolutionary analyses on phylogenetic trees. Methods Ecol. Evol. 7, 589–597 (2016).Article 

Google Scholar  More

We have analyzed the inbreeding-purging process in four captive populations of different ungulate species with effective sizes ranging 4–40 and with available pedigrees as well as survival and productivity records. This allows us to explore the role of inbreeding and purging in determining the evolution of fitness traits in a range of scenarios relevant in the context of conservation.In A. lervia (Ne ≈ 4), purging is expected only for the most severely deleterious alleles (those giving dNe  > 1, which implies d  > 0.25 as, for example, in completely recessive alleles with deleterious homozygous disadvantage s  > 0.5). Thus, it could be that purging has not been detected for this species because such severely deleterious alleles had been purged during the demographic decline in the wild, before the foundation of the captive population. This would be consistent with the low and non-significant inbreeding load estimated in this species. It is also possible that these estimates are non-significant due to the relatively small number of individuals available.G. cuvieri and N. dama have significant initial inbreeding loads that, adding up the direct and maternal components, is about 1.25 in both cases, which is on the order of other estimates published for captive populations (Ralls et al. 1988). Since in both species Ne  > 10, purging should be efficient against less severely deleterious alleles than in A. lervia (d  > 0.1). Purging is detected for both species with very low P values. This result is in agreement with Moreno et al. (2015), who suggested that purging had occurred in G. cuvieri as they found an increased juvenile survival parallel to an increased inbreeding coefficient. The relative contribution of severe and mild deleterious effects to the inbreeding load of populations is under a scientific debate with direct implications in conservation biology (Ralls et al. 2020, Kyriazis et al. 2021, Pérez-Pereira et al. 2021). The large d estimates obtained in our analysis indicate that a substantial fraction of the initial inbreeding load is being purged under modest effective population sizes, implying that such substantial fraction is due to relatively severe deleterious mutations in these two populations. As far as we are aware, these are the first estimates of this purging parameter obtained in managed, non-experimental populations. Previous estimates of d were obtained in D. melanogaster bottlenecked populations, first for egg-to-pupae viability in lines with Ne = 6 or 12 under noncompetitive conditions (d = 0.09, Bersabé and García-Dorado 2013), and second in lines with higher Ne ≈ 40–50 under more competitive conditions, giving a larger estimate of d, of the order of that estimated in these two ungulate endangered species (d ≈ 0.3, López-Cortegano et al. 2016).Regarding G. dorcas, given its larger population size, purging is expected even against alleles with mild recessive component of the deleterious effect (d  > 0.025). However, although a significant (if modest) inbreeding load was estimated, no significant purging was detected. Nevertheless, the number of equivalent complete generations by the end of the pedigree (EqG = 7) was smaller than our proposed minimum number of generations required to detect purging (tm = 10). This suggests that, due to the large size of this population, more generations are needed to detect purging.The results above support the use of tm to get an approximate idea about when a pedigree is too shallow for purging to be detected. Should the number of generations available be larger than tm, IP predictions could additionally be computed to search the d values that can be expected to produce detectable purging. Supplementary Fig. S3 shows that the true number of generations required to detect purging becomes increasingly larger than tm for alleles with smaller d values, as they suffer weaker purging each time they are exposed in homozygosis. The tm approach helps to understand the failure of many studies to detect purging. Such is the case of the extensive meta-analyses on 119 zoo populations by Boakes et al. (2007), where the median Ne value was 22.6 while the median number of generations was t = 3 meaning that, for most species, at least 5 more generations were needed before purging could be detectable. On the contrary, and in agreement with this tm approach, purging was experimentally detected in lines of D. melanogaster with Ne = 43 (i.e., tm ≈ 10) where, after an initial period of inbreeding depression, fitness experienced a substantial recovery beginning between generations 10 and 20 (López-Cortegano et al. 2016).A reason why detecting purging in captive populations is challenging is that a fitness rebound can also be due to adaptation to captive conditions or to environmental effects, such as those derived from improved husbandry (Clifford et al. 2007). In fact, this might have been the case in Speke’s gazelle breeding program, where the observed rebound of fitness was first ascribed to purging (Templeton and Read 1984, 1998), while Kalinowski et al. (2000) suggested that husbandry improvements could also be responsible for these findings. Our estimates of d and δ, however, are based on the association between the fitness trait and purged inbreeding at the individual level (Wi, gi) which, in our data, is mainly expressed within cohorts while average survival showed little variation through time. In addition, the analyses included temporal factors (YOB or POM) that should have removed confounding effects from adaptation to captivity or improved husbandry. Therefore, adaptive processes or time-dependent environmental factors are not expected to have biased our IP estimates.For productivity, the estimates of inbreeding load were high (overall inbreeding load ~5, P value  More

Insects and bamboosFive internodes (length: mean ± SD = 44.8 ± 1.1 cm, n = 5; diameter in the middle part of internodes: 21.4 ± 0.8 mm, n = 5) of five living mature culms of P. simonii bamboo were sampled at Kawaminami, Miyazaki Prefecture, Japan [32°9′ N, 131°29′ E] on 6 June, 2019. Per internode, four semi-cylindrical strips (ca. 15 × 2 cm) were made and stored at − 25 °C until use.To obtain fungus-free larvae of D. bucculenta, we sampled five beetle eggs from P. simonii bamboo obtained at Toyota, Aichi Prefecture, Japan [35°9′ N, 137°13′ E] on 9 May, 2019 in the laboratory from ovipositing females collected at Kawaminami on 10 and 11 April, 2019. The eggs were immersed in 99.5% ethanol for 10 s followed by 70% ethanol for 10 s for surface sterilization and then individually placed on potato dextrose agar (PDA) (Difco, Detroit, MI, USA) plates. The plates were incubated at 25 °C in the dark until 30 days after larval hatching to confirm the absence of the formation of yeast or other microbial colonies. Consequently, all five larvae hatched successfully and aseptically.The bamboo used in this study was morphologically identified using the literature29. This is native to the study areas and no other host bamboo species are distributed there29. Therefore, no voucher specimen of this bamboo has been deposited in a publicly available herbarium. No specific permits were required for the described field studies. The location is not privately-owned or protected in any way. The field studies did not involve endangered or protected species. All applicable international, national, and/or institutional guidelines for the care and use of animals and plants were followed. This study is reported in accordance with ARRIVE guidelines.Component analyses of bamboo tissuesFor YP and LP, the yeast W. anomalus originating from D. bucculenta in Kawaminami (strain: DBL05Kawaminami) was cultured on a 9-cm PDA plate to obtain enough biomass for further experiments. Afterwards, yeast cells were suspended in ca. 10 mL of sterilized water, and were inoculated on the inner surface of the autoclaved internode strips using an autoclaved tissue paper immersed with the yeast suspension. For LP, additionally, the fungus-free 2nd instar larvae (weight: mean ± SD = 2.4 ± 0.4 mg, n = 5) were individually placed on the yeast-inoculated strips. Each of these yeast-inoculated and yeast-and-larva-inoculated strips was then put in a sterilized test tube (3.0 cm in diameter and 20 cm tall) with moistened cotton placed at the bottom. Each of the test tubes was covered with a sterilized polypropylene cap, sealed with Parafilm Sealing Film (Pechiney Plastic Packaging, Chicago, IL, USA) on which three small holes were made using a fire-sterilized insect pin to avoid oxygen shortage, and individually put in a plastic zipper bag. These yeasts and insects were incubated at 25 °C in the dark for 47 days for YP (n = 5), and 47 (n = 4) and 73 (n = 1) days until these larvae reached adulthood for LP (adult elytral length: mean ± SD = 9.2 ± 0.4 mm, n = 5). Microbial contamination was invisible to the naked eye.For FP, YP and LP, the inner surface (up to 0.3 mm in thickness, dry weight: 336 to 935 mg) of a strip was sampled using a small U-shaped gouge. In the case of FX, first, the pith of a strip was completely removed, and then xylem tissue (up to 0.5 mm in thickness, dry weight: 729 to 872 mg) was sampled using a small U-shaped gouge. These tissues were individually sampled from five strips derived from five different internodes for each tissue type.Samples were extracted by aqueous ethanol and hydrolyzed by sulfuric acid with reference to the literature30,31,32 as follows. Four types of samples were freeze-dried and pulverized using a rotor-speed mill (Fritsch, PULVERISETTE 14, 0.2 mm mesh). About 80 mg of powdered sample was extracted using 5-mL 80% ethanol aqueous solution (aq.) at 63 °C three times. The volume of the extracts was adjusted to 25 mL, filtered, and analyzed using ion exchange chromatography measurements (extractable sugar analysis). Their extracted residues were hydrolyzed using sulfuric acid as follows: 50-mg samples were immersed in 1.64-g 72% sulfuric acid aq. at 30 °C for 2 h, boiled in 39.4-g 3% sulfuric acid aq. for 4 h, and filtered to collect sulfuric acid residues as sulfuric acid lignin fractions. The volumes of the filtrates were fixed to 100 mL, passed through a sulfuric acid-removing filter (DIONEX OnGuard IIA), and submitted to ion exchange chromatography measurements (structural sugar analysis). For the uronic acid measurements, the sulfuric acid-removing filter was not used.Ion exchange chromatography measurements were conducted using a DIONEX ICS-3000 apparatus. The measurement conditions were as follows: column, CarboPac PA-1 (2.0 mm I.D. × 250 mm L, Dionex corp.); flow rate, 0.3 mL min−1; column temperature, 30 °C; injection volume, 25 µL; eluent, H2O (solvent A), 100 mM NaOHaq. (solvent B), aqueous solution containing 100 mM NaOH and 1.0 M CH3COONa (solvent C), and aqueous solution containing 100 mM NaOH and 150 mM CH3COONa (solvent D). The gradient conditions for monomers, dimers, and uronic acids were as follows: for monomers, with a gradient of B 0.5% C 0% 45 min, C 100% 10 min, B 100% 10 min, B 0.5% C 0% 20 min; for dimers, with a gradient of B 50% C 0% 50 min, C 100% 10 min, B 100% 10 min, B 50% C 0% 15 min; for uronic acids, with a gradient of D 100% 10 min. These extraction, hydrolysis, and measurement procedures were conducted using n = 5 samples. For the structural sugars, their yield was calculated as the dehydrated state. The values of other extractives % were calculated by the subtraction of total extractable sugars % from total extractives %.Elemental analysis (carbon, hydrogen, nitrogen) was conducted by 2400 CHNS Organic Elemental Analyzer (PerkinElmer Japan, Yokohama, Japan). About 1-mg dried samples were burned completely and the produced CO2, H2O, and N2 (after reduction of NOx species) gasses were quantified by a thermal conductivity detector.Means of components of bamboo tissues were compared among tissue types using the Steel–Dwass test after the Kruskal–Wallis test. Calculations were performed using R 3.5.133.Carbon assimilation testThe yeast W. anomalus (DBL05Kawaminami) was cultured aerobically in 20 mL of yeast nitrogen base (YNB) (Difco) containing 0.5% glucose at 25 °C in the dark for 2 days with shaking at 85 rpm. The culture media were centrifuged and cell pellets were suspended in sterile water, in which the OD600 was adjusted to 0.10. Fifty μL of the cell suspension was added into a tube (2 mL) with 1 mL of each of 14 different media containing YNB and one of the following carbon sources: d-glucose, d-galactose, d-mannose, d-xylose, l-arabinose, d-fructose, d-galacturonic acid, d-glucuronic acid, sucrose, cellobiose, starch from corn, xylan from corn, carboxymethyl cellulose, and no carbon source (n = 5 to 6). The concentration of each carbon source was 0.5 g L−1, except for xylan at 1.5 g L−1. The tubes were shaken at 85 rpm and incubated at 25 °C in the dark for 7 days. Afterwards, the presence of visible pellets of yeasts and OD600 were recorded to determine the growth of the strain. The degree of assimilation was scored according to the presence of the pellets and the difference in the turbidity increase (ΔOD600) between culture media containing no and a given carbon source as follows: no growth (without a pellet, ΔOD600  More

1.Hughes, T. et al. Coral reefs in the Anthropocene. Nature 546, 82–90 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Moberg, F. & Folke, C. Ecological goods and services of coral reef ecosystems. Ecol. Econ. 29, 215–233 (1999).Article 

Google Scholar 
3.Brander, L. M., Van Beukering, P. & Cesar, H. S. The recreational value of coral reefs: A meta-analysis. Ecol. Econ. 63, 209–218 (2007).Article 

Google Scholar 
4.Jackson, J. B. et al. Historical overfishing and the recent collapse of coastal ecosystems. Science 293, 629–637 (2001).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Bruno, J. F., Sweatman, H., Precht, W. F., Selig, E. R. & Schutte, V. G. Assessing evidence of phase shifts from coral to macroalgal dominance on coral reefs. Ecology 90, 1478–1484 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
6.Hughes, T. Catastrophes, phase shifts, and large-scale degradation of a Caribbean Coral Reef. Science 265, 1547–1551 (1994).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Hughes, T. P., Bellwood, D. R., Folke, C. S., McCook, L. J. & Pandolfi, J. M. No-take areas, herbivory and coral reef resilience. Trends Ecol. Evol. 22, 1–3 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
8.Gardner, T. A., Côté, I. M., Gill, J. A., Grant, A. & Watkinson, A. R. Long-term region-wide declines in Caribbean corals. Science 301, 958–960 (2003).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
9.Edgar, G. J. et al. Reef Life Survey: Establishing the ecological basis for conservation of shallow marine life. Biol. Conserv. 252, 108855 (2020).Article 

Google Scholar 
10.Chabanet, P., Bigot, L., Garnier, R., Tessier, E. & Moyne-Picard, M. Coral reef monitoring at Reunion island (Western Indian Ocean) using the GCRMN method. Proc. 9th Int. Coral Reef Symp. 2, 873–878 (2000).
Google Scholar 
11.Lang, J. C., Marks, K. W., Kramer, P. A., Kramer, P. R. & Ginsburg, R. N. AGRRA Protocols Version 5.4. (2010).12.Cortés, J. et al. The CARICOMP network of Caribbean Marine Laboratories (1985–2007): History, key findings, and lessons learned. Front. Mar. Sci. 5, 519 (2019).Article 

Google Scholar 
13.Dethier, M. N., Graham, E. S., Cohen, S. & Tear, L. M. Visual versus random-point percent cover estimations: ‘objective’ is not always better. Mar. Ecol. Prog. Ser. 96, 93–100 (1993).Article 

Google Scholar 
14.Beijbom, O., Edmunds, P. J., Kline, D. I., Mitchell, B. G. & Kriegman, D. Automated annotation of coral reef survey images. In 2012 IEEE Conference on Computer Vision and Pattern Recognition 1170–1177 (IEEE, 2012).15.Beijbom, O. et al. Towards automated annotation of benthic survey images: Variability of human experts and operational modes of automation. PloS One 10, e0130312 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
16.Williams, I. D. et al. Leveraging automated image analysis tools to transform our capacity to assess status and trends of coral reefs. Front. Mar. Sci. 6, 222 (2019).Article 

Google Scholar 
17.Roelfsema, C. et al. Benthic and coral reef community field data for Heron Reef, Southern Great Barrier Reef, Australia, 2002–2018. Sci. Data 8, 1–7 (2021).MathSciNet 
Article 

Google Scholar 
18.Gonzalez-Rivero, M. et al. Monitoring of coral reefs using artificial intelligence: A feasible and cost-effective approach. Remote Sens. 12, 489 (2020).Article 

Google Scholar 
19.Cairns, S. D. Stony corals (Cnidaria: Hydrozoa, Scleractinia) of Carrie Bow Cay, Belize. Smithson. Contrib. Mar. Sci. 21, 271–302 (1982).
Google Scholar 
20.Rutzler, K. & Macintyre, I. G. The Atlantic Barrier Reef Ecosystem at Carrie Bow Cay, Belize, 1: Structure and Communities (Smithsonian Institution Press, 1982). https://doi.org/10.5479/si.01960768.12.539.Book 

Google Scholar 
21.Rutzler, K. Caribbean coral reef ecosystems: Thirty-five years of smithsonian marine science in Belize. In Proceedings of the Smithsonian Marine Science Symposium (2009).22.McField, M. et al. Mesoamerican Reef Report Card. (2020).23.Cox, C. E. et al. Genetic testing reveals some mislabeling but general compliance with a ban on herbivorous fish harvesting in Belize. Conserv. Lett. 6, 132–140 (2013).Article 

Google Scholar 
24.Pebesma, E. J. & Bivand, R. S. Classes and methods for spatial data in R. R News 5, 9–13 (2005).
Google Scholar 
25.Pebesma, E. Simple features for R: Standardized support for spatial vector data. R J. 10, 439–446 (2018).Article 

Google Scholar 
26.R Core Team. A language and environment for statistical computing. (2020).27.Althaus, F. et al. A standardised vocabulary for identifying benthic biota and substrata from underwater imagery: The CATAMI classification scheme. PLoS One 10, e0141039 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
28.Goatley, C. H., Bonaldo, R. M., Fox, R. J. & Bellwood, D. R. Sediments and herbivory as sensitive indicators of coral reef degradation. Ecol. Soc. 21(1), 29 (2016).Article 

Google Scholar 
29.Connell, S., Foster, M. & Airoldi, L. What are algal turfs? Towards a better description of turfs. Mar. Ecol. Prog. Ser. 495, 299–307 (2014).Article 

Google Scholar 
30.Lozada-Misa, P., Schumacher, B. D. & Vargas-Angel, B. Analysis of benthic survey images via CoralNet: A summary of standard operating procedures and guidelines. Pacific
Islands Fish. Sci. Cent. Natl. Mar. Fish. Serv. https://doi.org/10.7289/V5%2FAR-PIFSC-H-17-02 (2017).Article 

Google Scholar 
31.Obura, D. & Grimsditch, G. Resilience Assessment of Coral Reefs: Assessment Protocol for Coral Reefs, Focusing on Coral Bleaching and Thermal Stress (Citeseer, 2009).
Google Scholar 
32.Broeke, J., Pérez, J. M. M., & Pascau, J. Image processing with ImageJ. (Packt Publishing Ltd, 2015).
Google Scholar 
33.Wood, S. N. Generalized Additive Models: An Introduction with R (CRC Press, 2017).MATH 
Book 

Google Scholar 
34.Fasiolo, M., Nedellec, R., Goude, Y. & Wood, S. N. Scalable visualization methods for modern generalized additive models. J. Comput. Graph. Stat. 29, 78–86 (2020).MathSciNet 
Article 

Google Scholar 
35.Oksanen, J. et al. Community ecology package. R Package Version 2, (2013).36.Arnold, S. N. & Steneck, R. S. Settling into an increasingly hostile world: The rapidly closing ‘“Recruitment Window”’ for corals. PLoS One. https://doi.org/10.1371/journal.pone.0028681 (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
37.Mumby, P. J. & Harborne, A. R. Marine reserves enhance the recovery of corals on Caribbean reefs. PLoS One 5, e8657 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
38.Adam, T. C., Burkepile, D. E., Ruttenberg, B. I. & Paddack, M. J. Herbivory and the resilience of Caribbean coral reefs: Knowledge gaps and implications for management. Mar. Ecol. Prog. Ser. 520, 1–20 (2015).Article 

Google Scholar 
39.Suchley, A., McField, M. D. & Alvarez-Filip, L. Rapidly increasing macroalgal cover not related to herbivorous fishes on Mesoamerican reefs. PeerJ 4, e2084 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
40.Arnold, S. N., Steneck, R. S. & Mumby, P. J. Running the gauntlet: Inhibitory effects of algal turfs on the processes of coral recruitment. Mar. Ecol. Prog. Ser. 414, 91–105 (2010).Article 

Google Scholar 
41.Box, S. J. & Mumby, P. J. Effect of macroalgal competition on growth and survival of juvenile Caribbean corals. Mar. Ecol. Prog. Ser. 342, 139–149 (2007).Article 

Google Scholar 
42.Williams, I. & Polunin, N. Large-scale associations between macroalgal cover and grazer biomass on mid-depth reefs in the Caribbean. Coral Reefs 19, 358–366 (2001).Article 

Google Scholar 
43.Newman, M. J., Paredes, G. A., Sala, E. & Jackson, J. B. Structure of Caribbean coral reef communities across a large gradient of fish biomass. Ecol. Lett. 9, 1216–1227 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
44.Mumby, P. J., Steneck, R. S., Roff, G. & Paul, V. J. Marine reserves, fisheries ban, and 20 years of positive change in a coral reef ecosystem. Conserv. Biol. https://doi.org/10.1111/cobi.13738 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
45.Aronson, R., Precht, W., Toscano, M. & Koltes, K. The 1998 bleaching event and its aftermath on a coral reef in Belize. Mar. Biol. 141, 435–447 (2002).Article 

Google Scholar 
46.Green, D. H., Edmunds, P. J. & Carpenter, R. C. Increasing relative abundance of Porites astreoides on Caribbean reefs mediated by an overall decline in coral cover. Mar. Ecol. Prog. Ser. 359, 1–10 (2008).Article 

Google Scholar 
47.Roff, G., Joseph, J. & Mumby, P. J. Multi-decadal changes in structural complexity following mass coral mortality on a Caribbean reef. Biogeosciences 17, 5909–5918 (2020).Article 

Google Scholar 
48.Graham, N. & Nash, K. The importance of structural complexity in coral reef ecosystems. Coral Reefs 32, 315–326 (2013).Article 

Google Scholar 
49.Aronson, R. B. & Precht, W. F. White-band disease and the changing face of Caribbean coral reefs. Ecol. Etiol. New. Emerg. Mar. Dis. 159, 25–38 (2001).
Google Scholar 
50.Aronson, R. B., Macintyre, I. G., Precht, W. F., Murdoch, T. J. & Wapnick, C. M. The expanding scale of species turnover events on coral reefs in Belize. Ecol. Monogr. 72, 233–249 (2002).Article 

Google Scholar 
51.McField, M. et al. Status of the Mesoamerican Reef after the 2005 coral bleaching event. Status Caribb. Coral Reefs Bleach. Hurric. In 45–60 (2005).52.Arias-González, J. E. et al. A coral-algal phase shift in Mesoamerica not driven by changes in herbivorous fish abundance. PLoS One 12, e0174855 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
53.Miller, S. Automatically Annotating 175,000+ Images with the CoralNet API. CoralNet (Accessed 23 August 2021); https://coralnet.ucsd.edu/blog/automatically-annotating-175000-images-with-the-coralnet-api/ (2020).54.Muller, E. M., Sartor, C., Alcaraz, N. I. & van Woesik, R. Spatial epidemiology of the stony-coral-tissue-loss disease in Florida. Front. Mar. Sci. 7, 163 (2020).Article 

Google Scholar 
55.Alvarez-Filip, L., Estrada-Saldívar, N., Pérez-Cervantes, E., Molina-Hernández, A. & González-Barrios, F. J. A rapid spread of the stony coral tissue loss disease outbreak in the Mexican Caribbean. PeerJ 7, e8069 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
56.Weil, E. et al. Spread of the new coral disease “SCTLD” into the Caribbean: implications for Puerto Rico. Reef Encount. 34, 38–43 (2019).
Google Scholar 
57.Heres, M. M., Farmer, B. H., Elmer, F. & Hertler, H. Ecological consequences of Stony Coral Tissue Loss Disease in the Turks and Caicos Islands. Coral Reefs 40, 609–624 (2021).Article 

Google Scholar 
58.Walton, C. J., Hayes, N. K. & Gilliam, D. S. Impacts of a regional, multi-year, multi-species coral disease outbreak in Southeast Florida. Front. Mar. Sci. 5, 323 (2018).Article 

Google Scholar  More