More stories

  • in

    A microfluidic platform for highly parallel bite by bite profiling of mosquito-borne pathogen transmission

    Device design and fabricationVectorchips with integrated PDMS elastomer membranes were fabricated using a combination of laser patterning and soft lithography. We selected PDMS as the material of choice due to its low cost, biocompatibility, optical transparency, and chemical stability. PDMS-based devices can be fabricated using reproducible and scalable fabrication techniques26. Such devices have routinely been used for molecular analysis of small sample volumes using microfluidics27, and are compatible with molecular assays and cell cultures21,22. We used laser patterning to create open-ended arrayed compartments (3−4 mm deep) to capture saliva droplets from mosquito bites in PDMS micro-wells (Fig. 1a–d). In order to obtain ultra-thin PDMS films that the mosquito mouth parts can pierce, we utilized spin-coating of uncured PDMS over a flat silicon surface. We obtained membrane-integrated devices after plasma-induced PDMS-PDMS bonding of the laser-patterned array with the PDMS film (Fig. 1e). Detailed steps describing chip fabrication are included in the methods section below. The versatile fabrication process can provide user-defined variation in the size and density of the individual compartments. We were able to fabricate chips with hole diameters ranging from 150 μm to 1 cm (Supplementary Fig. 1) showing the scales of operation for the laser-ablation process while maintaining membrane stability. Additionally, the membrane can be easily integrated with fluidic networks for direct interfacing with mosquito bites, enabling assays involving on-chip fluid exchange (Fig. 1f). The microfluidic compartments on the chips can hold feeding media such as blood or sugar water (Fig. 1g), collect saliva during biting events, and act as isolated reaction chambers for molecular assays.Fig. 1: Fabrication of Vectorchip for collection of mosquito saliva and molecular assays.a Schematic showing the fabrication process of Vectorchip. Laser patterning was used to obtain through-holes in blocks of PDMS. A thin sacrificial layer of photoresist (Microposit™ S1813 G2, DOW Inc., USA) was spread onto a silicon wafer followed by spin-coating a thin film of PDMS on the wafer surface. The laser-patterned body of the chip was then bonded to the thin film of PDMS using plasma-activated PDMS-PDMS bonding. Membrane-bonded chips were obtained by placing the wafer in an acetone bath, where the sacrificial layer of resist was dissolved. b A PDMS block with laser-generated through-holes (diameter 250 μm). c A PDMS chip bonded to a 1.6 μm thin PDMS membrane on a 4-inch silcon wafer. d A Vectorchip with around 3000 wells (dia: 250 μm) is shown. e Chip with well diameter 1.75 mm and a visible PDMS membrane (thickness 1.6 μm). f Microfluidic channel-integrated chips with a PDMS membrane on top. g These wells can be loaded with feeding solution, reaction reagents, and can store mosquito saliva after biting assays. Wells loaded with feeding media (10% sucrose laced with blue food color) shown in the image. Scale bars represent 1 mm in (b), 5 mm in (e), and 2 mm in (f, g).Full size imageMosquito-chip interactionsNext, we examined the ability of mosquitoes to pierce through the PDMS membranes. Mosquitoes were attracted to the chips using heat as a guiding cue although several other baiting methods could be used. We placed a camera above the chip and observed stylet insertion through 1.6 μm thick PDMS membranes (Fig. 2a, Supplementary Movies 1 and 2). Abdominal engorgement in mosquitoes was observed after biting, indicating that mosquitoes can successfully feed through these membranes. Four species of mosquitoes were tested (Aedes aegypti, Aedes albopictus, Culex tarsalis, Culex quinquefasciatus) and all demonstrated successful probing and feeding through 1.6 μm thick PDMS membranes.Fig. 2: Mosquito biting on Vectorchip.a Stylet entry through a PDMS membrane for an Aedes aegypti female. b Mosquito feeding success as a function of membrane thickness for two mosquito species—Aedes aegypti and Culex tarsalis. We demonstrate that by changing the membrane thickness we can turn biting on or off. Biting off is defined as when no mosquito in the cage was able to feed from the chip in a duration of 30 min. The feeding assays were repeated at least three times (n = 3) at every given membrane thickness to confirm the ability of mosquitoes to feed through the membrane. c A Culex tarsalis female mosquito biting through a 20 μm thick PDMS membrane. d Bending of the proboscis observed for an Aedes aegypti female, indicating failure in biting through the 20 μm membrane. e Fluorescent salivary droplets expectorated by an Aedes aegypti mosquito while probing. These mosquitoes were fed on rhodamine-laced sugar water resulting in fluorescent saliva. Three salivary droplets (1, 2, and 3) are encircled. f Timelapse images show a magnified view of fluorescent salivary droplet deposition (droplets numbered 1, 2, and 3 as shown in panel (e) over a period of 4 s). Scale bars represent 500 μm in (a) and 2 mm in (c, d).Full size imageWe also examined whether the membrane thicknesses can be engineered to selectively prevent the biting of some species while letting other species feed through the barrier. We loaded the laser-patterned microfluidic compartments with blood and warmed the chips to attract mosquitoes (Supplementary Fig. 2 and Supplementary Movie 3). We considered biting to be off when none of the mosquitoes in the cage demonstrated abdominal engorgement within 30 min since the start of the experiment. Tests were performed with Aedes aegypti (~45 mosquitoes per cage) or Culex tarsalis (~25 mosquitoes per cage) while varying the thickness of the PDMS membrane from 900 nm to 100 μm. We observed significant differences in the biting ability of the two species, where complete inhibition of feeding by Aedes aegypti was observed with membrane thickness at or greater than 15 μm (Supplementary Movies 4 and 5). In contrast, complete inhibition of feeding by Culex tarsalis was only observed for membranes that were 100 μm thick (Fig. 2b). This remarkable difference in the biting capacity of these two species has not been reported earlier. While some studies have focused on understanding the biting mechanics of a single mosquito species (Aedes aegypti or Aedes albopictus)28,29, we are not aware of any study which quantifies the differences in biting strength between different species and the evolutionary or biomechanical differences involved in this variation. Evolutionary pressure on mosquitoes biting different hosts could generate great variability of biting strength in mosquito species, which may be needed to probe different skin thickness. This intriguing biomechanical phenomenon can be advantageous for deploying devices where membranes of desired thickness can act as species-selective filters, providing a means to better track highly relevant species and pathogens in the field.In order to extract molecular information from mosquito bites, we exploited the release of saliva during probing and biting events. The volume of saliva expectorated by mosquitoes in oil has been reported previously to be around 0.03 nL/min for Culex pipiens30, and 6.8 nL in approximately an hour for Aedes aegypti31 using a forced salivation method. Aedes aegypti have been estimated to salivate around 4.7 nL during blood feeding events32. We demonstrate that feeding rhodamine-labeled sugar water to mosquitoes renders their saliva fluorescent such that the release of saliva droplets can be directly visualized using fluorescence imaging (Fig. 2e and Supplementary Movie 6). Time sequences (Fig. 2f) show magnified images of three locations on a membrane surface where short probing events result in the deposition of fluorescent salivary droplets. We estimate that the volume of salivary droplets expectorated during this process is approximately 0.66 nL (Supplementary Fig. 3). This calculated mean volume is considerably lower compared to the volume associated with blood-feeding events reported earlier32; however, this discrepancy likely reflects differences between probing events such as observed in Fig. 2e, f and full feeding events as reported by Devine et al.32. The deposited salivary droplets harbor several biomarkers in the form of mosquito salivary proteins33, active pathogens, mosquito cells, and/or nucleic acids that can be used for their identification14,20,34.Distribution of bites on chipAn important consideration for molecular surveillance of mosquitoes is the granularity at which the population is sampled. Methods can range from pooling strategies with low coverage where samples are pooled together at the population scale, to the individual collection of saliva from mosquitoes while recording species identity. Vectorchip potentially provides an adjustable method between these two strategies by increasing or decreasing the density of bite wells presented to a population of mosquitoes.We first analyzed the distribution of bites on Vectorchip by tracking mosquito activity during biting assays. In order to perform sucrose biting assays for molecular diagnostics, we filled the wells with 10% sucrose, and covered the open side of the chip with a glass slide. A resistor that acts as a heat source was placed atop the glass slide and the chip-resistor assembly was placed on top of the mesh ceiling of a mosquito cage. A part of the chip (50 percent of the total chip area) was protected by paper tape such that mosquitoes could not bite into that section of the device. The area of the chip protected by tape served as negative control for the downstream molecular assays. A camera (Raspberry Pi 3 camera module V2.1) was placed on the bottom of the cage to track and record individual mosquito activity on chip (Supplementary Fig. 4). Image recognition and mosquito tracking algorithms have been described by us previously25, where we used computer vision to analyze image sequences and robustly track the presence, movement, and biting behavior of individual mosquitoes at high resolution. In this study, we have used an imaging setup where the camera is attached to the bottom of the cage (Supplementary Fig. 4a). The obtained image datasets reveal mosquito trajectories on the chip as well as areas with prolonged interactions (Fig. 3a–d, Supplementary Fig. 4, and Supplementary Movie 7). Unique trajectories are identified as movements performed by an individual mosquito while in the field of view (e.g., landing, exploration, biting, and take-off by a single mosquito would constitute a trajectory). The movement of mosquitoes on a chip is shown in Fig. 3e–g and Supplementary Fig. 4. Depending on the duration of mosquito-chip interactions, multiple trajectories do not overlap (Fig. 3e), partially overlap (Fig. 3f), or completely overlap (Fig. 3g). Trajectory plots can indicate the likelihood that a set of wells has been probed by more than one mosquito. These results indicate that by modulating the time of interaction between the mosquito population and the chip, the user can control the degree of probing experienced by the chip. As an example, a chip with no overlap between trajectories should exclusively contain wells probed by single mosquitoes. Tracking mosquito movement on chip was used to locate trajectories and thus regions with unique mosquito presence.Fig. 3: Tracking mosquito activity on chip.a Image recognition software was used to identify and track mosquitoes on the chip. Wells are highlighted in yellow. b Trajectories of mosquito movement can be plotted and indicate the overlap between different movement tracks. c Location plot of mosquitoes on the chip are indicated by red circles. Regions with overlapping red circles are darker in color. d Trajectory and location plot can be used to identify regions with individual mosquito locomotion activity as compared to regions with more crowding. e−g Varying experimental parameters (symbols indicate the duration of experiments and number of mosquitoes) can be used to obtain either individual mosquito plots or population data. Three cases are discussed with e no overlap between multiple probing trajectories, f partial overlap, and g complete overlap. h Representative image for a chip where the presence of fluorescent droplets on wells has been indicated by a red dot. i Zoomed-in fluorescent image of the same chip showing wells with single, or multiple droplets. j Number of bites per well were calculated by counting fluorescent droplets deposited on the wells after biting events. k We defined the number of mosquitoes = Nm, the number of wells = Nw, probing frequency = η, and time = t. The probability of obtaining ‘k’ bites on a well for the parameters (Nm, Nw, η, t) utilized in (h−j) has been shown. Scale bars represent 5 mm in (a−c, e−g), 2 mm in (d), and 1 mm in (i).Full size imageWhile tracking provided macro-scale information about mosquito position and activity, we also utilized the release of fluorescent salivary droplets deposited by biting mosquitoes to resolve mosquito probing events at higher resolution. We quantified the number of fluorescent spots seen in every well of the PDMS chip by using fluorescent scans of the chip before and after biting (Fig. 3h). No spots were observed in the negative control region where access was blocked to mosquitoes, while several fluorescent spots observed on the open areas of the chip indicate the wells which have been probed. The images show the presence of wells with multiple bite marks, as well as single bite marks.We further examined the design parameters used to build Vectorchip that can dictate the resolution of this sampling strategy. The primary determinants of how often a single mosquito will bite a specific well are dependent on multiple factors, including the number of mosquitoes in the cage, the time mosquitoes have access to the chip and biting behavior (e.g., frequency), and the size and number of wells on the chip. Other additional factors like the feeding media may affect the number of biting events performed by a female mosquito. In this study, we focused on sucrose as the feeding media since it is a diet source that does not inhibit polymerase chain reaction (PCR) assays. We considered a simple mathematical model to arrive at a first approximation for the distribution of bites on a chip. In this model, we made a simple assumption that mosquitoes probe the chip at a fixed frequency (η). We also assumed that mosquito bites are random and independent events, since no spatial gradients have been created on the chip. Detailed analysis of this model can be found in Supplementary Section 1.We defined the number of mosquitoes = Nm, the number of wells = Nw, and probing frequency = η, and estimated the probability (P) that a well receives k bites in given time t as (Supplementary Note 1):$$P[k]=({N}_{m}eta t)!/(k!times ({N}_{m}eta t-k)!)times 1/{N}_{w}^{k}times {(1-1/{N}_{w})}^{{N}_{m}eta t-k}$$
    (1)
    This was approximated by the form (Supplementary Note 1):$$P[k]approx {{{{{{mathrm{e}}}}}}}^{-{N}_{m}eta t/{N}_{w}}times {({N}_{m}eta t/{N}_{w})}^{k}/k!$$
    (2)
    and represents a Poisson distribution of the form,$$P[k]={{{{{{mathrm{e}}}}}}}^{-lambda }times {(lambda )}^{k}/k!$$
    (3)
    with the expected rate of events (i.e., bites per well), λ = Nmηt/Nw. Probability distributions for the expected number of bites on a well were plotted while varying the number of mosquitoes, wells, and probing frequencies (Supplementary Fig. 5). The probing frequency of Aedes aegypti mosquitoes on a parafilm membrane has been estimated previously to be in the range of 0–1 bites per min35. Supplementary Fig. 5 shows a diverse phase space where a single bite per well is the most dominant outcome, i.e., λ was less than or equal to 1, when up to 50 mosquitoes interact with a 150-well Vectorchip for 15 min—in line with our experiments. In the context of a larger number of mosquitoes interacting with the chip, increasing the number of wells can ensure that the majority of wells will not get bitten more than once. We quantified bites on wells as indicated from fluorescent droplets (Supplementary Table 1), and compared it to the analytical model (Fig. 3j, k). The fitting parameter was η = 0.1 bites per minute, which is close to the range of probing frequencies reported recently35. With water as the feeding media and parafilm as the probing membrane, Jove et al. report that the majority of tested mosquitoes showed a probing frequency below 0.5 bites per minute35. While we observe and model slightly lower probing frequencies, this difference is possible due to the change in the membrane material from parafilm to PDMS, and other differences in the setup. The methods discussed above provide multiple ways to understand mosquito−chip interactions as well as help us define parameters to obtain data at single-bite resolution.PCR diagnostics on-chipHaving established that mosquitoes can bite through the membrane and deposit salivary droplets in microwells while biting on a Vectorchip, we analyzed these bites for their DNA and RNA content. Running PCR-based nucleic acid amplification directly on individual wells establishes a low-cost assay to detect bite content in a high-throughput manner. We utilized devices with a PDMS membrane thickness of 1.6 μm and well diameter of 1.75 mm to perform molecular analysis of salivary samples; providing 145 independent reactions per chip (within a chip area of 5 cm × 2.5 cm).Detection of mosquito mitochondrial DNA (mtDNA) and RNA from viruses in saliva was performed using on-chip PCR with end-point fluorescence readout. A detailed protocol for performing PCR in Vectorchip is provided in the Methods section and Supporting Information file (Supplementary Fig. 6). We tested the efficiency of PCR amplification in Vectorchip by manually loading a known concentration of DNA and RNA into the reaction wells (Fig. 4a–c). We performed the PCR for 42 cycles and detect DNA amplification from approximately 5 DNA copies in 4 μL of reaction mix (~1 copy/μL) (Fig. 4b). Simultaneously detection of Zika virus RNA using reverse transcription-PCR (RT-PCR), was successfully demonstrated from as low as 15 RNA copies in 4 μL of reaction mix (Fig. 4c). PCR reactions were also performed in 96-well plates to test if amplification accuracy and sensitivity were similar as compared to Vectorchips. Both qPCR amplification curves and end point fluorescence information was obtained from the well plates. The 96-well plate reactions were performed in DI water and also in the presence of dried sucrose (10%) to ensure that the amplification response provided a good comparison to the Vectorchip-based reactions using manually spiked concentrations, and to test if the presence of sucrose would inhibit the reactions (see Supplementary Fig. 7a, b).Fig. 4: On-chip PCR for detection of mosquito and pathogen.The symbols indicate manually spiked assays, uninfected mosquito assays, and Zika infected mosquito assays from top to bottom. a–c A spiked assay with successful amplification of Aedes aegypti DNA and Zika RNA on chip. a Indicates wells filled with PCR mix, with ROX dye providing the background fluorescence. Fluorescent wells indicate amplification of b mosquito DNA, and c Zika RNA on chip. Assays were performed with uninfected mosquitoes where d tracking patterns were collected and demonstrate extensive translocation activity on the chip. e PCR shows detection of mosquito DNA on the chip directly from biting. f Percentage of wells available for biting on Vectorchip that display a positive PCR outcome. The PCR assay was repeated 6 times on different chips (n = 6). The false positive rate (amplification in wells that were covered by tape) was close to zero (0.23%). Assays performed with Zika infected mosquitoes indicate the presence of both g Aedes aegypti DNA and h Zika virus RNA after bites on chip. i Percentage of wells available for biting that show a positive signature for Zika RNA, mosquito DNA or both. The PCR assay was repeated 3 times on different chips (n = 3). Scale bars represent 5 mm in (a−e) and 3 mm in (g, h). Source data are provided as a Source Data file.Full size imageVectorchip provides an opportunity for tracking both pathogens in salivary bites, as well as identifying host species by detecting DNA from the host biting the well. Mosquito DNA has previously been reported in salivary droplets deposited by feeding Culex mosquitoes20. Utilizing uninfected Aedes aegypti mosquitoes, we performed PCR assays on-chip to detect the presence of mosquito DNA in deposited saliva droplets. Figure 4d, e shows the tracking data and PCR detection of mtDNA for a chip, which was placed on a cage with 75 uninfected Aedes aegypti females for 45 min. Areas highlighted in green were protected by paper tape such that mosquitoes could not probe them and serve as negative controls. The tracking data obtained from the chip indicate significant activity of the mosquitoes on-chip as represented by several overlapping trajectories. PCR data confirms that mosquito DNA indeed can be detected in probed wells. Multiple experiments indicate that a subset (3−12%, n = 6) of the wells where mosquitoes are active test positive for the presence of mosquito DNA (Fig. 4f). The rate of false positives obtained from wells not accessible to mosquitoes was very close to zero (0.23%, n = 6). In order to further validate the rate of detection of mosquito DNA on Vectorchips, we performed a biting assay using a 96 well plate loaded with agarose and covered with a parafilm membrane. Agarose gel was used as feeding media during this test, to minimize the risk of cross-contamination between wells during the peeling off of the parafilm membrane (details in “Methods” section). We recorded the tracking and DNA amplification response from uninfected mosquito bites on the well plate (Supplementary Fig. 7c, d). Approximately 15% of the wells where mosquitoes were active tested positive for mosquito DNA, which is close to the fractions observed on Vectorchips. Since mosquito DNA can only be detected in a subset of wells showing mosquito activity, this suggests that the presence of detectable levels of mosquito DNA in saliva is a noisy physiological process, likely dependent on multiple variables such as duration of feeding and time of last bite. It is furthermore important to note that our tracking algorithm only detects the presence of mosquitoes, but does not provide information regarding if a well was bitten or not.Next, we utilized the device to perform assays on infected mosquitoes, for the detection of both mosquito species and deposited pathogens. Figure 4h shows results obtained when Aedes aegypti infected with Zika virus (ZIKV)36 interacted with the chip for 20 min. We subsequently performed RT-PCR on the chip and observed that Aedes aegypti mtDNA and ZIKV RNA could be detected in 24/118 and 37/118 wells, respectively. We summarize the detection rates of mosquito DNA and ZIKV RNA for three assays performed with ZIKV-infected mosquitoes in Fig. 4i. Interestingly, not all wells positive for ZIKV RNA showed positive for mosquito DNA and vice versa, highlighting the stochastic genetic composition of bite-derived nanoliter saliva droplets. A higher fraction of wells (mean ~ 22%, n = 3) showed positive for ZIKV RNA as compared to mosquito mtDNA (mean ~ 17%, n = 3), indicating a higher prevalence of viral genomic material in mosquito saliva as compared to mosquito mtDNA. Detection of mosquito DNA and virus RNA in these assays demonstrates the capacity of this tool to identify mosquito and pathogen species directly from mosquito bites on-chip.Detection of infectious viral particles directly from bitesWhile PCR-based end-point assays are an important strategy for multiplexed detection of vectors and pathogens in various settings, they cannot determine the presence or concentration of infectious virus particles in bites. Detection of infectious pathogen load in mosquito bites is important to understand the vector competence (including potential and efficiency of pathogen transmission through bites) of mosquito-pathogen systems and its dependency on factors such as local climate37, mosquito species38, physiology39, and presence of endosymbionts (e.g., Wolbachia)40. Typical methods to measure vector competence rely on manual “forced salivation” of individual mosquitoes12, where the viral loads may differ from biting events and mosquitoes are sacrificed preventing time-course measurements. Sugar feeding on filter papers has also been utilized recently, and provides promising results towards longitudinal sampling of viral transmission by individual mosquitoes20. As compared to filter paper-based approaches, the Vectorchip uses an integrated PDMS membrane which can enable purely bite-based collection and avoid possible contamination from e.g., excretion events. Furthermore, our method improves the throughput of the assay as mosquitoes do not need to be individually housed and sampled. Finally, filter paper assays are reliant on sampling nucleic acids, and cannot perform live cell culture assays to detect infectious pathogens.We tested the potential of Vectorchip to perform focus forming assays (FFA) and detect infectious viral particles directly as a result of mosquito bites on chip (Fig. 5). Vectorchip wells were filled with cell culture media (Dulbecco’s modified eagle medium (DMEM) supplemented with 10% fetal bovine serum (FBS). The efficiency of FFA on-chip was tested using manual spiking of viral particles into Vectorchip (Supplementary Fig. 8). In order to perform biting assays, Aedes aegypti mosquitoes were infected with Mayaro virus, an emerging zoonotic pathogen that causes a dengue-like disease41. The mosquitoes were able to bite through the membrane into the wells and transfer live Mayaro viral particles directly into the cell culture medium (Fig. 5a). These viral particles then infect a monolayer of vero cells cultured on the PDMS membrane. These infections can be identified via fluorescent antibodies specific to viral envelope proteins (Fig. 5c). Each fluorescent patch (focus) is attributed to an infection resulting from a single viral particle. Foci counted on two sample chips can be seen in Fig. 5d. No foci were visible from negative control samples. The distribution of viral foci on the chips indicates the heterogeneity of infectious particle dose in mosquito bites. Figure 5 shows that the Vectorchip supports the growth of a cell monolayer on the chip and can enable the direct collection of mosquito saliva through bites in cell culture media. Antibody assays can be used to directly quantify infectious viral particles in the wells. These results demonstrate the suitability of using the Vectorchip to probe the transmission of infectious viral particles.Fig. 5: Quantifying active viral particles using focus-forming assays.Focus forming assays on Vectorchip can directly quantify the number of active viral particles in mosquito bites, without the need for isolation of individual mosquitoes and manual salivary extraction. a Image shows Aedes aegypti mosquitoes infected with Mayaro virus biting on Vectorchips filled with DMEM cell culture media. b A monolayer of vero cells was cultured on Vectorchip membrane. The formation of vero cell monolayers on membranes was verified on 5 independent chips (n = 5). c FFA performed on a Vectorchip. Symbols indicate wells in the top row which were bitten by mosquitoes with no infection (negative control), and rows of wells bitten by mosquitoes infected with Mayaro virus. The green channel shows fluorescence from an antibody against a viral envelope protein. Every green island represents infected cells likely to have resulted from an active viral particle. No active viral particles were seen in the control wells. FFA on Vectorchips was verified on 3 chips (n = 3) for both uninfected and Mayaro infected mosquitoes. d FFA formed on two chips and viral foci were counted using the antibody fluorescence. Scale bars represent 2 mm in (a), 30 μm in (b), and 100 μm in (c).Full size imageFinally, we examined if the feeding media provided in Vectorchip can enable blood meal-mimic biting behavior through the addition of phagostimulants such as ATP, while maintaining accurate molecular analysis (Supplementary Note 2). We observed that the addition of 100 μM ATP, (which is a phagostimulant and used by mosquito sensory neurons to identify blood35), to the feeding media did not impact nucleic acid amplification reactions (Supplementary Fig. 9a) (n = 3). Interestingly, mosquitoes feeding on DMEM supplemented with 10% FBS and ATP (1 mM) led to a similar or higher number of abdominal engorgements in 45 min as compared to blood meals (n = 3) (Supplementary Fig. 9b). This is reasonable as DMEM (supplemented with FBS and ATP) contains a variety of the components used by mosquito sensory neurons to identify blood at relevant concentrations (Supplementary Note 2)35). In summary, our data indicate that Vectorchip can provide physiologically relevant mosquito feeding behavior, enabling accurate assessment of vector and pathogen identity through nucleic acid amplification reactions for field surveillance, as well as on-chip assessment of transmission of infectious viral particles in bites through focus forming assays. More

  • in

    Diet diversity and environment determine the intestinal microbiome and bacterial pathogen load of fire salamanders

    1.Ley, R. et al. Evolution of mammals and their gut microbes. Science 320, 1647–1651 (2008).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    2.Robinson, C. J., Bohannan, B. J. M. & Young, V. B. From structure to function: The ecology of host-associated microbial communities. Microbiol. Mol. Biol. Rev. 74, 453–476 (2010).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    3.Gill, S. R. et al. Metagenomic analysis of the human distal gut microbiome. Science 312, 1355–1359 (2006).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    4.Pryor, G. & Bjorndal, K. Symbiotic fermentation, digesta passage, and gastrointestinal morphology in bullfrog tadpoles (Rana catesbeiana). Physiol. Biochem. Zool. 78, 201–215 (2005).PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    5.Claus, S. P., Guillou, H. & Ellero-Simatos, S. The gut microbiota: A major player in the toxicity of environmental pollutants?. NPJ Biofilms Microbiomes 2, 16003 (2016).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    6.Mazmanian, S. K., Liu, C. H., Tzianabos, A. O. & Kasper, D. L. An immunomodulatory molecule of symbiotic bacteria directs maturation of the host immune system. Cell 122, 107–118 (2005).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    7.Alberdi, A., Aizpurua, O., Bohmann, K., Zepeda-Mendoza, M. L. & Gilbert, M. T. P. Do vertebrate gut metagenomes confer rapid ecological adaptation?. Trends Ecol. Evol. 31, 689–699 (2016).PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    8.Bourguignon, T. et al. Rampant host switching shaped the termite gut microbiome. Curr. Biol. 28, 649-654.e2 (2018).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    9.Amato, K. et al. Evolutionary trends in host physiology outweigh dietary niche in structuring primate gut microbiomes. ISME J. 13, 1 (2018).
    Google Scholar 
    10.Sullam, K. E. et al. Environmental and ecological factors that shape the gut bacterial communities of fish: A meta-analysis. Mol. Ecol. 21, 3363–3378 (2012).PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    11.Bolnick, D. I. et al. Individuals’ diet diversity influences gut microbial diversity in two freshwater fish (threespine stickleback and Eurasian perch). Ecol. Lett. 17, 979–987 (2014).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    12.Grond, K., Sandercock, B. K., Jumpponen, A. & Zeglin, L. H. The avian gut microbiota: Community, physiology and function in wild birds. J. Avian Biol. 49, e01788 (2018).Article 

    Google Scholar 
    13.Michel, A. et al. The gut of the finch: Uniqueness of the gut microbiome of the Galápagos vampire finch. Microbiome 6, 1–14 (2018).Article 

    Google Scholar 
    14.Delsuc, F. et al. Convergence of gut microbiomes in myrmecophagous mammals. Mol. Ecol. 23, 1301–1317 (2014).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    15.Carmody, R. N. et al. Diet dominates host genotype in shaping the murine gut microbiota. Cell Host Microbe 17, 72–84 (2015).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    16.Kohl, K., Amaya, J., Passement, C., Dearing, M. D. & Mccue, M. Unique and shared responses of the gut microbiota to prolonged fasting: A comparative study across five classes of vertebrate hosts. FEMS Microbiol. Ecol. 90, 883–894 (2014).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    17.Vences, M. et al. Gut bacterial communities across tadpole ecomorphs in two diverse tropical anuran faunas. Sci. Nat. 103, 25 (2016).Article 
    CAS 

    Google Scholar 
    18.Li, G. et al. Host-microbiota interaction helps to explain the bottom-up effects of climate change on a small rodent species. ISME J. 14, 1795–1808 (2020).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    19.Rawls, J., Mahowald, M., Ley, R. & Gordon, J. Reciprocal gut microbiota transplants from zebrafish and mice to germ-free recipients reveal host habitat selection. Cell 127, 423–433 (2006).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    20.Bletz, M. C. et al. Amphibian gut microbiota shifts differentially in community structure but converges on habitat-specific predicted functions. Nat. Commun. 7, 13699 (2016).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    21.Woodhams, D. C. et al. Host-associated microbiomes are predicted by immune system complexity and climate. Genome Biol. 21, 23 (2020).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    22.Adlerberth, I. & Wold, A. E. Establishment of the gut microbiota in Western infants. Acta Paediatr. Int. J. Paediatr. 98, 229–238 (2009).CAS 
    Article 

    Google Scholar 
    23.Wu, G. D. et al. Linking long-term dietary patterns with gut microbial enterotypes. Science 334, 105–108 (2011).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    24.Stuart, S. N. et al. Status and trends of amphibian declines and extinctions worldwide. Science 306, 1783 (2004).ADS 
    CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    25.Lips, K. R. et al. Emerging infectious disease and the loss of biodiversity in a Neotropical amphibian community. Proc. Natl. Acad. Sci. USA. 103, 3165 (2006).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    26.Bishop, P. et al. The amphibian extinction crisis -what will it take to put the action into the amphibian conservation action plan?. Surv. Perspect. Integr. Environ. Soc. 5, 97–111 (2012).
    Google Scholar 
    27.Kats, L. & Ferrer, R. Alien predators and amphibian declines: Review of two decades of science and the transition to conservation. Divers. Distrib. 9, 99–110 (2003).Article 

    Google Scholar 
    28.Chanson, J., Hoffman, M., Cox, N. & Stuart, S. The State of the World’s Amphibians. In Threatened Amphibians of the World 33–44 (Lynx Edicions, Barcelona, Spain, 2015)29.Rollins-Smith, L. A. & Woodhams, D. C. Amphibian immunity: Staying in tune with the environment. In Ecoimmunology ( eds Demas, G. & Nelson, R.) 92–143 (Oxford University press, Oxford, UK, 2011).30.Martel, A. et al. Recent introduction of a chytrid fungus endangers Western Palearctic salamanders. Science 346, 630 (2014).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    31.Birnie-Gauvin, K., Peiman, K. S., Raubenheimer, D. & Cooke, S. J. Nutritional physiology and ecology of wildlife in a changing world. Conserv. Physiol. 5, cox030 (2017).PubMed 
    PubMed Central 

    Google Scholar 
    32.Scheele, B. C. et al. Amphibian fungal panzootic causes catastrophic and ongoing loss of biodiversity. Science 363, 1459 (2019).ADS 
    CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    33.Whiles, M. R. et al. The effects of amphibian population declines on the structure and function of Neotropical stream ecosystems. Front. Ecol. Environ. 4, 27–34 (2006).Article 

    Google Scholar 
    34.Hocking, D. & Babbitt, K. Amphibian contributions to ecosystem services. Herpetol. Conserv. Biol. 9, 1–17 (2014).
    Google Scholar 
    35.Burton, T. M. & Likens, G. E. Energy flow and nutrient cycling in salamander populations in the Hubbard Brook experimental forest, New Hampshire. Ecology 56, 1068–1080 (1975).CAS 
    Article 

    Google Scholar 
    36.Reagan, D. P. & Waide, R. B. The Food Web of a Tropical Rain Forest (University of Chicago Press, 1996).
    Google Scholar 
    37.Stebbins, R. C. & Cohen, N. W. A Natural History of Amphibians (Princeton University Press, 1997).
    Google Scholar 
    38.Flecker, A. S., Feifarek, B. P. & Taylor, B. W. Ecosystem engineering by a tropical tadpole: Density-dependent effects on habitat structure and larval growth rates. Copeia 1999, 495–500 (1999).Article 

    Google Scholar 
    39.Beard, K., Vogt, K. & Kulmatiski, A. Top-down effects of a terrestrial frog on nutrient dynamics. Oecologia 133, 583–593 (2002).ADS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    40.Davic, R. & Welsh, H. On the ecological role of salamanders. Annu. Rev. Ecol. Syst. 12, 405–434 (2004).Article 

    Google Scholar 
    41.Reinhardt, T., Steinfartz, S., Paetzold, A. & Weitere, M. Linking the evolution of habitat choice to ecosystem functioning: Direct and indirect effects of pond-reproducing fire salamanders on aquatic-terrestrial subsidies. Oecologia 173, 281–291 (2013).ADS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    42.Buckley, D. & Alcobendas, M. Salamandra salamandra (Linnaeus, 1758). (2002).43.Fryxell, J. & Lundberg, P. Diet choice and predator—prey dynamics. Evol. Ecol. 8, 407–421 (1994).Article 

    Google Scholar 
    44.Deagle, B. E. et al. Studying seabird diet through genetic analysis of faeces: A case study on macaroni penguins (Eudyptes chrysolophus). PLoS ONE 2, e831 (2007).ADS 
    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    45.Botzler, R. G., Wetzler, T. F. & Cowan, A. B. Yersinia enterocolitica and yersinia-like organisms isolated from frogs and snails. Bull. Wildl. Dis. Assoc. 4, 110–115 (1968).CAS 
    Article 

    Google Scholar 
    46.Cooper, J. E., Needham, J. R. & Griffin, J. A bacterial disease of the Darwin’s frog (Rhinoderma darwini). Lab. Anim. 12, 91–93 (1978).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    47.Hird, D. et al. Enterobacteriacae and Aeromonas hydrophila in Minnesota frogs and tadpoles (Rana papiens). Appl. Environ. Microbiol. 46, 1423–1425 (1984).ADS 
    Article 

    Google Scholar 
    48.Olson, M., Gard, S., Brown, M., Hampton, R. & Morck, D. Flavobacterium indologenes infection in leopard frogs. J. Am. Vet. Med. Assoc. 201, 1766–1770 (1992).CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    49.Pearson, M. D. Motile Aeromonas septicaemia of farmed Rana spp. (1998).50.Green, S. et al. Identification and management of an outbreak of Flavobacterium meningosepticum infection in a colony of South African clawed frogs (Xenopus laevis). J. Am. Vet. Med. Assoc. 214(1833–8), 1792–1793 (1999).
    Google Scholar 
    51.Bernardet, J.-F. et al. Polyphasic study of Chryseobacterium strains isolated from diseased aquatic animals. Syst. Appl. Microbiol. 28, 640–660 (2005).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    52.Pasteris, S., Guidoli, M., Otero, M., Bühler, M. & Nader-Macías, M. In vitro inhibition of Citrobacter freundii, a red-leg syndrome associated pathogen in raniculture, by indigenous Lactococcus lactis CRL 1584. Vet. Microbiol. 151, 336–344 (2011).PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    53.Kirk, K. et al. Chryseobacterium angstadtii sp. nov., isolated from a newt tank. Int. J. Syst. Evol. Microbiol. 63, 4777–4783 (2013).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    54.Suzina, N. E. et al. Cytophysiological characteristics of the vegetative and dormant cells of Stenotrophomonas sp. strain FM3, a bacterium isolated from the skin of a Xenopus laevis frog. Microbiology 87, 339–349 (2018).CAS 
    Article 

    Google Scholar 
    55.Hallinger, M., Taubert, A. & Hermosilla, C. Endoparasites infecting exotic captive amphibian pet and zoo animals (Anura, Caudata) in Germany. Parasitol. Res. 119, 3659–3673 (2020).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    56.Deagle, B. E., Chiaradia, A., McInnes, J. & Jarman, S. N. Pyrosequencing faecal DNA to determine diet of little penguins: Is what goes in what comes out?. Conserv. Genet. 11, 2039–2048 (2010).Article 

    Google Scholar 
    57.Deagle, B. E., Thomas, A. C., Shaffer, A. K., Trites, A. W. & Jarman, S. N. Quantifying sequence proportions in a DNA-based diet study using Ion Torrent amplicon sequencing: Which counts count?. Mol. Ecol. Resour. 13, 620–633 (2013).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    58.Nakahara, F. et al. The applicability of DNA barcoding for dietary analysis of sika deer. DNA Barcodes 3, 200–206 (2015).Article 

    Google Scholar 
    59.Deagle, B., Kirkwood, R. & Jarman, S. Analysis of Australian fur seal diet by pyrosequencing prey DNA in faeces. Mol. Ecol. 18, 2022–2038 (2009).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    60.Pompanon, F. et al. Who is eating what: Diet assessment using next generation sequencing. Mol. Ecol. 21, 1931–1950 (2012).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    61.Thomas, A. C., Jarman, S. N., Haman, K. H., Trites, A. W. & Deagle, B. E. Improving accuracy of DNA diet estimates using food tissue control materials and an evaluation of proxies for digestion bias. Mol. Ecol. 23, 3706–3718 (2014).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    62.Deagle, B. & Tollit, D. Quantitative analysis of prey DNA in pinniped faeces: Potential to estimate diet composition?. Conserv. Genet. 8, 743–747 (2007).CAS 
    Article 

    Google Scholar 
    63.Ando, H. et al. Methodological trends and perspectives of animal dietary studies by noninvasive fecal DNA metabarcoding. Environ. DNA 2, 391–406 (2020).Article 

    Google Scholar 
    64.Deagle, B. et al. Molecular scatology as a tool to study diet: Analysis of prey DNA in scats from captive Steller sea lions. Mol. Ecol. 14, 1831–1842 (2005).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    65.Parsons, K., Piertney, S., Middlemas, S., Hammond, P. & Armstrong, J. DNA-based identification of salmonid prey species in seal faeces. J. Zool. 266, 275–281 (2005).Article 

    Google Scholar 
    66.Meekan, M., Jarman, S., McLean, C. & Schultz, M. DNA evidence of whale sharks (Rhincodon typus) feeding on red crab (Gecarcoidea natalis) larvae at Christmas Island, Australia. Mar. Freshw. Res. 60, 607–609 (2009).CAS 
    Article 

    Google Scholar 
    67.Guillerault, N., Bouletreau, S., Iribar, A., Valentini, A. & Santoul, F. Application of DNA metabarcoding on faeces to identify European catfish Silurus glanis diet. J. Fish Biol. 90, 2214–2219 (2017).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    68.Brown, D. S., Jarman, S. N. & Symondson, W. O. C. Pyrosequencing of prey DNA in reptile faeces: Analysis of earthworm consumption by slow worms. Mol. Ecol. Resour. 12, 259–266 (2012).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    69.Ferenti, S., Cicort-Lucaciu, A. S., Dobre, F., Paina, C. & Covaci, R. The food of four Salamandra salamandra populations from Defileul Jiului National Park (Gorj County). Olten. Stud. Si Comun. Stiintele Nat. 2008, 153–160 (2008).
    Google Scholar 
    70.Ferenti, S., David, A. & Nagy, D. Feeding-behaviour responses to anthropogenic factors on Salamandra salamandra (Amphibia, Caudata). Biharean Biol. 4, 139–143 (2010).
    Google Scholar 
    71.Lezău, O. et al. The feeding of two Salamandra salamandra (Linnaeus, 1758) populations from Jiului Gorge National Park (Romania), South West. J. Hortic. Biol. Environ. 1, 143–152 (2010).
    Google Scholar 
    72.Balogová, M., Maxinová, E., Orendáš, P. & Uhrin, M. Trophic spectrum of adult Salamandra salamandra in the Carpathians with the first note on food intake by the species during winter. Herpetol. Notes 8, 371–377 (2015).
    Google Scholar 
    73.Sebastiano, S., Antonio, R., Fabrizio, O., Dario, O. & Roberta, M. Different season, different strategies: Feeding ecology of two syntopic forest-dwelling salamanders. Acta Oecologica 43, 42–50 (2012).ADS 
    Article 

    Google Scholar 
    74.Lunghi, E. et al. Field-recorded data on the diet of six species of European Hydromantes cave salamanders. Sci. Data 5, 1–7 (2018).Article 

    Google Scholar 
    75.Lunghi, E. et al. What shapes the trophic niche of European plethodontid salamanders?. PLoS ONE 13, e0205672 (2018).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    76.Measey, G. Diet of feral Xenopus laevis (Daudin) in South Wales, UK. J. Zool. 246, 287–298 (1998).Article 

    Google Scholar 
    77.Le, D. T. T., Rowley, J. J., Tran, D. T. A. & Hoang, H. D. The diet of a forest-dependent frog species, Odorrana morafkai (Anura: Ranidae), in relation to habitat disturbance. Amphib. Reptil. 41, 29–41 (2020).Article 

    Google Scholar 
    78.Pamintuan, P. E. & Starr, C. K. Diet of the giant toad, Bufo marinus (Amphibia: Salientia), in a coastal habitat of the Philippines. Trop. AgricTrinidad 93, 323–327 (2016).
    Google Scholar 
    79.Plummer, M. & Farrar, D. Sexual dietary differences in a population of Trionyx muticus. J. Herpetol. 15, 175–179 (1981).Article 

    Google Scholar 
    80.Shetty, S. & Shine, R. Activity patterns of yellow-lipped sea Kraits (Laticauda colubrina) on a Fijian island. Copeia 2002, 77–85 (2002).Article 

    Google Scholar 
    81.Vincent, S., Herrel, A. & Irschick, D. Sexual dimorphism in head shape and diet in the Cottonmouth Snake (Agkistrodon piscivorus). J. Zool. 264, 53–59 (2004).Article 

    Google Scholar 
    82.Manenti, R., Conti, A. & Pennati, R. Fire salamander (Salamandra salamandra) males’ activity during breeding season: Effects of microhabitat features and body size. Acta Herpetol. 12, 29–36 (2017).
    Google Scholar 
    83.Keen, W. H. Feeding and activity patterns in the salamander Desmognathus ochrophaeus (Amphibia, Urodela, Plethodontidae). J. Herpetol. 13, 461–467 (1979).Article 

    Google Scholar 
    84.Forester, D. C. Parental care in the salamander Desmognathus ochrophaeus: Female activity pattern and trophic behavior. J. Herpetol. 15, 29–34 (1981).Article 

    Google Scholar 
    85.Harris, W. E. Spermatophore deposition behaviour in an explosive breeder, the Small mouthed salamander, Ambystom texanum. Herpetologica 64, 149–155 (2008).Article 

    Google Scholar 
    86.Anderson, T. & Mathis, A. Diets of two sympatric neotropical salamanders, bolitoglossa mexicana and B. rufescens, with notes on reproduction for B. rufescens. J. Herpetol. 33, 601 (1999).Article 

    Google Scholar 
    87.Shu, Y. et al. Comparison of intestinal microbes in female and male Chinese concave-eared frogs (Odorrana tormota) and effect of nematode infection on gut bacterial communities. MicrobiologyOpen 8, e00749 (2019).PubMed 
    Article 
    CAS 
    PubMed Central 

    Google Scholar 
    88.Zhou, J. et al. A comparison of nonlethal sampling methods for amphibian gut microbiome analyses. Mol. Ecol. Resour. 20, 844–855 (2020).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    89.Huang, C. & Liao, W. Seasonal variation in gut microbiota related to diet in Fejervarya limnocharis. Animals 11, 1393 (2021).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    90.Chang, C.-W., Huang, B.-H., Lin, S.-M., Huang, C.-L. & Liao, P.-C. Changes of diet and dominant intestinal microbes in farmland frogs. BMC Microbiol. 16, 33 (2016).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    91.Kohl, K. D., Cary, T. L., Karasov, W. H. & Dearing, M. D. Restructuring of the amphibian gut microbiota through metamorphosis. Environ. Microbiol. Rep. 5, 899–903 (2013).PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    92.Colombo, B. M., Scalvenzi, T., Benlamara, S. & Pollet, N. Microbiota and mucosal immunity in amphibians. Front. Immunol. 6, 111–111 (2015).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    93.Novoslavskij, A. et al. Major foodborne pathogens in fish and fish products: a review. Ann. Microbiol. 66, 1–15 (2015).Article 

    Google Scholar 
    94.Standish, I. et al. Yersinia ruckeri isolated from common mudpuppy necturus maculosus. J. Aquat. Anim. Health 31, 71–74 (2019).PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    95.Hird, D. W. et al. Enterobacteriaceae and Aeromonas hydrophila in Minnesota frogs and tadpoles (Rana pipiens). Appl. Environ. Microbiol. 46, 1423–1425 (1983).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    96.Heiman, M. L. & Greenway, F. L. A healthy gastrointestinal microbiome is dependent on dietary diversity. Mol. Metab. 5, 317–320 (2016).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    97.Amato, K. & Righini, N. The howler monkey as a model for exploring host-gut microbiota interactions in primates.https://doi.org/10.1007/978-1-4939-1957-4_9 (2015).98.Kartzinel, T. R., Hsing, J. C., Musili, P. M., Brown, B. R. P. & Pringle, R. M. Covariation of diet and gut microbiome in African megafauna. Proc. Natl. Acad. Sci. 116, 23588 (2019).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    99.Tiede, J., Scherber, C., Mutschler, J., McMahon, K. D. & Gratton, C. Gut microbiomes of mobile predators vary with landscape context and species identity. Ecol. Evol. 7, 8545–8557 (2017).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    100.Peig, J. & Green, A. J. New perspectives for estimating body condition from mass/length data: The scaled mass index as an alternative method. Oikos 118, 1883–1891 (2009).Article 

    Google Scholar 
    101.Vences, M. et al. Freshwater vertebrate metabarcoding on Illumina platforms using double-indexed primers of the mitochondrial 16S rRNA gene. Conserv. Genet. Resour. 8, 323–327 (2016).Article 

    Google Scholar 
    102.Amaral-Zettler, L. A., McCliment, E. A., Ducklow, H. W. & Huse, S. M. A method for studying protistan diversity using massively parallel sequencing of V9 hypervariable regions of small-subunit ribosomal RNA genes. PLoS ONE 4, e6372 (2009).ADS 
    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    103.Caporaso, J. G. et al. QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336 (2010).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    104.Amir, A. et al. Deblur rapidly resolves single-nucleotide community sequence patterns. mSystems 2, e00191-16 (2017).PubMed 
    PubMed Central 

    Google Scholar 
    105.Klindworth, A. et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 41, e1 (2013).CAS 
    PubMed 
    Article 
    PubMed Central 

    Google Scholar 
    106.Aguirre, A. A. et al. The One Health Approach to toxoplasmosis: Epidemiology, control, and prevention strategies. EcoHealth 16, 378–390 (2019).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    107.Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    108.Quast, C. et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596 (2012).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    109.Vavrek, M. J. Fossil: Palaeoecological and palaeogeographical analysis tools. Palaeontol. Electron. 14, 16 (2011).
    Google Scholar 
    110.McMurdie, P. J. & Holmes, S. phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8, e61217 (2013).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    111.Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    112.Jaccard, P. The distribution of the flora of the Alpine zone. New Phytol. 11, 37–50 (1912).Article 

    Google Scholar 
    113.Oksanen, J. et al. vegan: Community Ecology Package. R package version 2.5-5. 2019 (2020).114.Dray, S. & Dufour, A.-B. The ade4 package: Implementing the duality diagram for ecologists. J. Stat. Softw. 22, 1–20 (2007).Article 

    Google Scholar  More

  • in

    Precision management of pollination services to blueberry crops

    Experimental designWithin the new varieties obtained in the market the variety “Emerald” is a good model for pollination service management. In Argentina the plant is vigorous and very productive, and can produce fruit in autumn without reducing the spring production through the long flowering period41. As of 2020, around 260 hectares are planted with this blueberry variety in Entre Ríos province.Field work was conducted during the 2017 flowering period (from July 17th to August 23rd, 2017) on 9 blueberry plots (Vaccinium corymbosum var. Emerald) located in the region of Concordia (Entre Rios, Argentina). Three plots were pollinated under “precision (precision) honey bee management”, while 6 other plots were pollinated under conventional (conventional) honey bee management (two conventional plots per precision plot). We decided to double the number of control plots per farm (conventional bee management) in order to reflect the variability that these agroecosystems present in the study region.Each plot comprised a 1-ha area, planted with the Emerald variety. The plants were all the same age and subjected to the same agricultural management in respect to fertilization, pruning, water, etc. Six of the nine plots (2 precision plots and the corresponding 4 conventional ones) were open air, while 3 plots (one precision and two conventional) were covered with mesh. This mesh is used to prevent hail and wind damage. All the plots were pollinated by honey bees (Apis mellifera L.) at a density of 10 hives per cultivated ha (Fig. 5).Figure 5Map of the plots of blueberries studied. The fields with conventional and precision management of bees are detailed in each farm, and the location of the hives used in each particular case. One Precision (yellow) and two Conventional (the numbers stand out C1 and C2, in red) plots were sampled per farm.Full size imagePrecision honey bee managementOne month before the blooming season, we carefully selected 30 beehives based on: (a) number of combs covered with adult bees, (b) number of combs with brood, and (c) the presence of an active egg-laying queen, to be prepared for the pollination service. The selected stock of beehives to be used in the precision plots came from the same stock of bees to be used in the control plots, and thus, we could expect that any differences regarding to bees’ performance in the field would be related to the different management (precision) and not due to differences regarding the stocks of bees, neither the ecotypes nor the beekeeping management previous to ours. During the month prior to the pollination service, the selected beehives were prepared for the blooming period by being fed with a dietary supplement (nutritional technology exclusively licensed by Beeflow S.A. and Beeflow Corporation) to boost their immune systems and strength. When blooming reached 5–10% in the precision plots, the beehives were placed so as to surround each plot and managed individually. Beehives were checked during the blooming period and, if it was necessary, they were also fed. We also coordinated pesticide spray applications with growers to avoid bee losses during the pollination service. When 95% of the female flowers were senescent the hives were retired from the field and returned to their original apiary.Conventional bee managementThis group of beehives represented the conventional management performed by beekeepers in Concordia Region, which consists of beehives placed in big groups, between 5 to 50, and located either next to or up to hundreds of meters from the field, with high variation in terms of hive quality (see above), and with infrequent servicing. This management style is similar to those performed by many beekeepers around the globe.In summary, the main differences between the management of bee hives for pollination services in blueberry are: Start of the incentive diet, the hives for precision management began to be fed before the blooming period with a dietary supplement (sugar syrup 1: 1) to encourage the collection of pollen; Protein feeding, the hives for precision handling were fed during the entire flowering period of the blueberry with a diet especially rich in protein and essential amino acids (i.e.,: alpha arginine); Spatial arrangement, strategic and homogeneous location of the hives individually in contrast to distribution in medium groups (5 hives) to large (30 hives) located in strategic places according to the production advisor; Management against agrochemical applications on the plot, Preventive management for pesticides, fungicides and fertilizers applied by means of aqueous suspension, in this case through an agricultural spraying turbine for a tractor, was carried out only in the hives for precision pollination.Data collectionThe visitation rate to flowers was estimated simultaneously in precision and conventional plots, and at different times of the day. To do this, pollinator censuses were performed by recording the number of flowers visited by pollinators for a period of 5 min (i.e., no. visits . flower−1 . 5 min−1). Each census involved a different randomly selected group of flowers (between 10 and 30 floral units per census), within a randomized location within the plantation. These visits censuses were carried out in all the lots studied on the same day between 10 am and 6 pm, periodically changing the visiting hours to each lot to cover the entire hourly range throughout the 15 sampling dates distributed in between June 17 and August 23, 2017. Census of visits to flower were carried out on days with a favorable climate for foraging bees, sunny to partially sunny days at times where the temperature it is suitable for bee activity ( > 10° C). In total, we performed 70 pollinator censuses per plot distributed in 15 days, totaling 52.5 h of observation (i.e. 5.8 hs per plot × 9 plots).In each plot, 5 plants were tagged (i.e., 5 plants per plot * 9 plots, 45 plants in total), and within each plant 3 branches were randomly selected and tagged with paper tape (a total of 225 experimental branches). The number of open flowers between the tag and the end of the branch were counted in order to determine the number of flowers that produced a fruit (i.e., fruit set). This experimental design allowed the estimation of variation in flowers within and among plants and plantations, subjected to the different pollination practices.When the time of the harvest arrived, the number of fruits (corresponding to each tagged branch) developed from each plant was counted in relation to the total number of flowers originally present on the branch to estimate fruit-set. Then, the quality of all tagged fruits harvested was evaluated by quantifying individual weight and firmness. A total of ~ 6700 flowers and ~ 5200 fruits were analyzed.To estimate the firmness of the fruits we used the “compression force,” which is a widely used indicator, using a texture analyzer TA.XT Plus (Stable Micro Systems Ltd., United Kingdom) equipped with a 5 kg load cell and a 75 mm cylinder aluminum probe. It is defined as the force necessary to deform the fruit by 10% and is used as a quantifiable measure equivalent to the force produced when a fruit is pressed between the thumb and the index finger. Finally, to estimate individual fruit weight we used an electronic scale.Statistical analysisWe analyzed the effect of honey bee management practice (conventional vs. precision) on: (i) visitation frequency to blueberry flowers, (ii) fruit set, (ii) fruit weight, and (iii) fruit firmness, using general and generalized mixed-effects models. Data analysis was carried out using the lmer function from the “lme4” package42,43 of R software (version 2.15.1), applying separate analyses and distributions for each response variable (e.g. visits frequency, fruit set, and firmness, see below).For (i) visitation frequency (i.e., bee visits) as the response variable the model assumed a Poisson error distribution with a log-link function. We included the number of flowers observed in each census as an offset (i.e., a fixed predictor known in advance to influence insect visitation)44. For (ii) fruit set (i.e., proportion of flowers setting fruit) as response variable the model assumed a Binomial error distribution with a logit-link function. And finally, for (iii) fruit weight and (iii) fruit firmness as response variables the models assumed a Gaussian error distribution. For all models, honey bee management practice (conventional vs precision) was included as a fixed effect, each “plant” was nested within “field”, and each “field” was nested within “farm” as a random effect, following the hierarchical design of our experimental design (i.e. 3 farms, 3 fields within each farm, and 5 plants within each field). More

  • in

    Resistance, resilience, and recovery of salt marshes in the Florida Panhandle following Hurricane Michael

    Based on comparative analysis of aerial imagery, the marshes in the Florida Panhandle were overwhelmingly undamaged by Hurricane Michael. Of the 173,259 km2 of marsh analyzed across Bay, Gulf, and Franklin counties, only 1.9% (3371 km2) was classified as damaged after Hurricane Michael. Less than 2% of marshes were damaged in Bay and Franklin counties, and less than 5% of marshes were damaged in Gulf County. This result suggests that the majority of marshes in the study counties can withstand the effects of a category 5 hurricane and continue to provide coastal protection ecosystem services29. This also identifies marshes as more resistant than most other coastal defenses, including bulkheads. While this study did not examine bulkhead failure, similar studies found that 76% of bulkheads surveyed in North Carolina were damaged after Hurricane Irene29, and a large portion of bulkheads in the Florida Keys experienced significant damage after Hurricane Irma31,32.Damage to marshes in the Florida Panhandle from Hurricane Michael was spatially distributed largely due to the storm track. Shortly before making landfall in Mexico Beach, the storm turned northeastward, influenced by the southern edge of mid-latitude westerlies1. The highest hurricane wind speeds are typically associated with the front-right quadrant of the storm and weaken rapidly after landfall, which was also true for Hurricane Michael20. As expected, the majority of damaged marshes (4.3%) were in Gulf County, which is directly southeast of landfall. Given the sizable decrease in wind speed and storm surge as the storm lost intensity over land, about half of the damage (51.6% across all counties) occurred within 50 m of the coast, and the vast majority of damage (95.1% across all counties) occurred within 500 m of the coast.The types of damage were consistent with hurricane impacts in the Gulf of Mexico, with sediment deposition and wrack deposits among the most common morphological impacts from hurricane strikes21,22. The type of damage was also spatially distributed. Marshes in Bay County experienced the highest wind speeds and also experienced t he most damage from fallen trees, a largely wind-driven form of damage. Marshes in Franklin County experienced the largest inundations, as well as the most conversion to open water. Deposition of sediment or vegetation can be influenced by both wind and inundation, whether through aeolian transport3,4,5,7,8,23 or overwash12. Therefore, it is not surprising that deposition was both the most common as well as the most evenly spatially-distributed form of marsh damage. Vegetation loss, though less common than deposition, is also largely driven by high-velocity wind-driven currents; a particularly intense flow can pluck or denude a marsh of its vegetation24,25.Marsh damage is also likely closely tied to physical properties of the marsh, in addition to the spatial component of storm conditions. The vast majority (88.5%) of marsh damage occurred at elevations of less than 1 m (NAVD88), consistent with both average marsh elevation and areas of greatest susceptibility to inundation. Localized susceptibility to damage may also indicate differences in marsh soil composition or shear strength24,26 or previous physical disturbance27. Given this study’s reliance on visual observations to determine damage, it is important to consider that not all storm-induced marsh damage can be ascertained visually; subsurface processes such as shallow subsidence or expansion can greatly influence the marsh’s elevation and resilience14,28.It is also important to consider that, of the ~ 2% of marshes in the Florida Panhandle that were damaged, ~ 80% experienced deposition of sediment or vegetation on the marsh surface, which is a potentially less-permanent form of damage, especially compared to fallen trees or conversion to open water. Hurricanes regularly place thick sediment deposits on marshes in the Gulf of Mexico; short-term sedimentation rates in coastal Louisiana marshes after Hurricane Andrew increased for three months post-storm by 1–3 orders of magnitude compared to pre-storm rates, suggesting that resuspension and deposition continued well beyond the timeframe of the passing storm33. Deposition of sediment from storm overwash may actually benefit the system, increasing total marsh elevation and counteracting sea-level rise8. While storm-induced sedimentation may increase resilience, deposits that are too thick ( > 5–10 cm34) may cause plant mortality and ultimately reduce the stability of the marsh and its resilience to storm impacts.Much of the deposition on the marsh surface was rafted vegetation and wrack mats, which impacts marsh vegetation differently than sediment. This wrack deposition is consistent with past hurricane impacts in the Gulf of Mexico; wrack deposits from Hurricane Andrew completely buried vegetation in coastal Louisiana marshes, and areas of especially thick wrack were slow to recolonize, even a year after landfall35. A 1995 study, however, found that only 30% of wrack mats on a New England salt marsh damaged underlying vegetation36, and while wrack mats were a major cause of damaged low-marsh vegetation in a Virginia salt marsh, vegetation typically recovers37. This suggests that even the 2% of damaged marshes in the study area may be more resilient than originally thought.Though the vast majority of marshes in the Florida Panhandle were not visually damaged by Hurricane Michael, of the damaged marshes with aerial imagery from six months after the storm, only 16% recovered. Marshes exposed to less extreme environmental conditions during the storm often were more likely to recover; this was the case for the relationship between mean maximum wind speed and deposition recovery, as well as maximum inundation and vegetation loss and conversion to open water.Perhaps unsurprisingly, some forms of damage are easier to recover from than others. More than 45% of marshes experiencing vegetation loss had recovered within six months, whereas only 16% of marshes experiencing deposition of sediment or wrack, the most common form of damage, had visually recovered within 6 months. This is consistent with other studies; marshes in the Mississippi River Delta denuded by Hurricane Camile partially or completely recovered within one growing season, provided the root mat was not destroyed and the area was not permanently submerged12. This study is limited to 6 months after landfall, which is a relatively short time frame and does not include the spring and summer growing seasons. Given that other studies have shown that it can take over a year to recover marsh vegetation after disturbance (including storms and fire)38, and it is likely that more marsh recovery would be observed a full year after disturbance. Despite these considerations, this study provides important insight into the short-term recovery of vegetation (through recolonizing after plucking or growing through thin layers of deposition).Once vegetation is substantially disturbed or submerged, however, recovery becomes more difficult, and it is these areas that might be prioritized for more active recovery interventions. Kirwan et al. (2008) found that disturbed areas experience decreased vertical accretion, which in turn causes localized submergence of the marsh platform, as well as potential expansion of channel networks, further destabilizing the marsh. Less than 4% of marsh that had been converted to open water and none of the marsh experiencing channel widening or cutting had recovered within 6 months. This is consistent with damage identified from a survey of hurricanes impacting southern Louisiana over the last 50 years12. Since newly formed or widening ponds or channels indicate significant loss of the marsh platform, it is far more difficult for the marsh to both recover elevation and regrow vegetation. Research has shown that ponds created by hurricane impacts could maintain their shape for decades and eventually become a permanent feature39.This study found that though the vast majority of marshes in Bay, Gulf, and Franklin counties were not visually damaged during Hurricane Michael, of the 2% of marshes that did experience damage, 84% had not recovered within six months. This suggests that marshes may be largely resistant to storm impacts but not particularly resilient19, at least within the first six months of disturbance. Recovery rates for marshes on public land were significantly higher than private land. This speaks to the benefits of public land management and the need to create better incentives for marsh management on private lands. Additionally, the majority of marshes on public land are within larger wildlife refuges, state parks, or management areas and, as a result, may be less altered and more resilient to storm impacts.Future damage surveys incorporating on-the-ground assessments and finer-resolution aerial imagery are recommended to fully understand the impacts of intense storms on salt marshes. Since this study, by design, only quantifies damage and recovery visible from aerial imagery, it is important to consider that there may be additional damage to the marshes that may be small, subtle, or belowground and therefore not observable at the scale of the imagery. While this study is specific to the inundation and hydrodynamic conditions associated with Hurricane Michael, the implications about marsh performance and recovery provide an important base for similar future assessments of future storm impacts.Post-storm management actions may be key to increasing marsh resistance and resilience after storm impacts. As seen in this study, deposition is the one of the most common forms of marsh damage, but also one of the most easily reversible. While storm-induced sediment deposition is often largely beneficial to building marsh elevation31, wrack deposits may damage vegetation36,37, which in turn must recover by the next growing season to remain resilient to the next hurr icane season. Removing the thickest layers of wrack buildup to prevent vegetation death is a relatively straightforward management activity, which in turn will build the marsh’s resilience to future near-term storm impacts. Cost-effective restoration and recovery efforts should prioritize maintaining vegetation health, particularly in locations where the marsh platform is still intact post-storm. Post-storm wrack removal is a target for recovery funding that may be particularly beneficial.Results of this study inform policy, funding, and approaches for hazard mitigation, disaster recovery, adaption, and conservation. Very few marshes were damaged or failed during a direct hit from one of the strongest hurricanes to impact the continental United States, and marsh failure rate was generally much lower than that seen for other artificial defenses, even in lower-intensity storms30,31,32. A portion of damaged marshes also began to recover within six months with few or no interventions, which artificial defenses cannot. Managed marshes on public lands are more likely to be both resistant and resilient to storm damage. Previous work has shown that marshes significantly reduced hurricane damage to property15 and are a cost-effective method of coastal protection40. Combined with our results, which indicate that marshes are highly resistant to storm impacts, these studies identify the need for additional public and private incentives, including insurance incentives and hazard mitigation funding, to conserve and restore marshes for their benefits to both people and property. After a storm, allocation of recovery funds to this natural, national infrastructure is imperative to improving management and recovery trajectories, particularly given the potential for more intense hurricane landfalls as the climate changes19. More

  • in

    The flowering of Atlantic Forest Pleroma trees

    Study siteThe study covered the Brazilian Atlantic Forest domain (Fig. 2), which is located on the east coast of Brazil between latitudes 5° and 30° south, expanding over 500 km inland in the south. It consists of a total area of 1,085,151 km2 with limits defined by the Brazilian Ministry of the Environment48. The total area covered by Sentinel-2 tiles overlapping with the Atlantic Forest domain is ~2,006,959 km2. This latter area was used to compute the descriptive statistics of detections.DataSentinel 2 imagesThe pink or magenta blossoms of Pleroma trees were mapped using Sentinel-2 multi-spectral data with 10 m spatial resolution taken approximately every five days under the same viewing conditions. We used only Sentinel-2 images with Level-1C correction—which are orthoimage products, i.e. map projections of acquired images using a digital elevation model to correct ground geometric distortions—and delivered in images of 100 km × 100 km. Pixel radiometric measurements were provided in Top-Of-Atmosphere (TOA) reflectances (coded in 12 bits)49.In the analysis, 213 Sentinel-2 tiles covering the Brazilian Atlantic Forest domain were used, totaling 2,006,959 km2 which is equivalent to ~20 billion Sentinel-2 pixels with 10 m spatial resolution (Fig. 2a). Amongst the 213 selected tiles, 36 had 2 orbits to download to obtain the full tile image due to the overlapping orbit paths (called replicates in the following text).For each tile and replicate (213 + 36), the times series between 31 June 2016 and the 1 July 2020 was downloaded from the Google Cloud Storage Sentinel-2 repository (https://cloud.google.com/storage/docs/public-datasets/sentinel-2). To reduce the dataset size, we retained only images with less than 80% cloud cover; and, from the month outside the flowering months of the Pleroma trees (July to November), we kept only images with less than 25% cloud cover. The complete dataset was made up of 33,798 Sentinel-2 images.Four spectral bands available at 10 m spatial resolution were used: Red (665 nm), Green (560 nm), Blue (490 nm) and NIR (842 nm). A border of 120 pixels with NA values was added to the image to produce images of 10240 × 10240 pixels to ease automation of the image analysis workflow, which generally works with 2n × 2n size pixel images. In our case here, the deep learning analysis was made with 128 × 128 pixel images and an additional 8 × 8 border. Sentinel L1C reflectance values are in the range of 0–10000 and were converted to 8 bits (0–254) with the following rules : for Red, Green and Blue bands, we kept the minimum value between 2540 and the original pixel value, divided this value by 10 and converted the result to integer; and for the NIR band, we keep the minimum value between 2540 and the original pixel value divided by a constant equaling 3.937, divided this value by 10 and converted the result to integer. While it was not expected to have RGB pixel values for vegetation with reflectance above 2540, it occured frequently for the NIR values. Dividing the NIR band values by the constant 3.937 enabled scaling the full range of the original NIR values between 0 and 2540 without losing too much information. For each tile, all 4 bands were saved in one GeoTIFF of 8 bits to ease storage and processing. The size of the complete dataset was 5.59 teraoctets. The automatic download, scaling and conversion of the images to 8 bits took about 25 days (from 16 July 2020 to 3 August 2020 and from 10 September 2020 to 13 September 2020).Environmental dataFigure 12Environmental and climatic variables used in the study to analyse spatial distribution of Pleroma trees (a) elevation (m), (b) slope (°), (c) tree cover (%), (d) mean annual precipitation (mm yr−1), (e) annual mean of minimum temperatures (°C), and (f) maximal annual temperatures (°C).Full size image
    To test the association of Pleroma trees with elevation and slope, elevation data from the Shuttle Radar Topography Mission (SRTM) were used50 (Fig. 12a). Specifically, we used the 3 arc-seconds (~90 m) spatial resolution digital elevation database (version 4) provided by the CGIAR Consortium for Spatial Information51. This dataset, in comparison to the original NASA STRM dataset, has been processed to fill data voids. From this dataset, we used the variables elevation (m) and computed slope (°) considering the four neighbor pixels (Fig. 12b). To analyse the relationship between Pleroma trees presence and forest tree cover, we used the tree cover percentage for the year 2000 at 30 m of spatial resolution, which we obtained from the global forest cover dataset (Fig. 12c), which is based on Landsat time series52.The association of Pleroma trees with local climate was tested using the annual means of precipitation and air temperatures (Fig. 12d–f). The mean annual precipitation over the study period was computed from the CHIRPS v2p0 monthly precipitation dataset at 0.05° of spatial resolution produced by University of California, Santa Barbara (UCSB). CHIRPS data are global rainfall estimates from rain gauges and satellite observations53. The mean of maximum and minimum air surface temperatures over the study period were computed from the Aqua/AIRS L3 Daily Standard Physical Retrieval (AIRS-only) at 1° of spatial resolution V7.0 (AIRS3STD). AIRS, the Atmospheric Infrared Sounder on NASA’s Aqua satellite, gathers daily infrared energy emitted from Earth’s surface and atmosphere globally and provides 3D measurements of temperature and water vapor through the atmospheric column54. The annual mean of minimum and maximum air surface temperatures was calculated using the daily air surface temperature measured from the descending orbital pass, which occurs at 1:30 am local time (’SurfAirTemp_D’), and the ascending orbital pass, which occurs at 1:30 pm (’SurfAirTemp_A’).Additionally, maps produced by the Brazilian Institute of Geography and Statistics (IBGE) of the geomorphological units and rivers of Brazil were used to describe the spatial distribution of the blossom detections33.All environmental variables were resampled to a raster of 1280 × 1280 m spatial resolution using an average interpolation to match the resolution of the Pleroma tree detection dataset.ModelNeural network architectureThis detection model is a deep learning neural network (Fig. 13), more specifically an encoder with a VGG16-like structure35, that given an image (input image) return the probability of presence of Pleroma trees with flowers in the image. The model inputs are 4 bands RGB-NIR images made up of 136 × 136 pixels at 10 m of spatial resolution (Fig. 13). Sentinel-2 tiles of 10240 × 10240 pixels were cropped based on a regular grid of 128 × 128 pixels, and 4 neighbouring pixels were added on each side to create an overlap between the patches. The resulting images are 136 × 136 pixels in size. However, in the training, the presence or absence of blooming Pleroma was given only for the images of 128 × 128 pixels without consideration of the borders. This enable to avoidance of the border effect that is common in convolutional neural networks. Each image of 136 × 136 pixels goes through a data augmentation process that consists in random vertical and horizontal flips. No additional data augmentation necessary due to the natural data augmentation provided by atmospheric conditions and illumination. After data augmentation, the images were then fed to the detection encoder. The encoder was made up of 5 consecutive convolution and pooling blocks, one fully connected layer (dense 100) and a final output layer with a softmax activation that provided the probability of presence of blooming Pleroma trees in the image (Fig. 13). Additionally, one drop-out layer was used at the end of the fully connected layer to perform further implicit data augmentation and avoid overfitting during training. The model has a total of 25,448,686 parameters, of which 25,440,622 are trainable. The model was coded in R language55 with Rstudio interface to Keras and TensorFlow 2.256,57,58,59.Figure 13Architecture of the Pleroma blossom detection model.Full size imageNetwork trainingTo make the training sample, a manual sample was produced for the Sentinel-2 tile 23KMQ, in the area where we had previously made a high resolution map of blooming Pleroma6, and for five other tiles where flowering Pleroma were detected visually from high resolution Google Earth images (22JFQ, 22JGQ, 23KLP, 23KLQ and 23KNQ, respectively). What is identified in the Sentinel-2 images are forest stands dominated by Pleroma and not single individuals. Pleroma trees have a small stature (8–12 m height) and crown of less than 10 m and one tree alone cannot influence sufficiently the reflectance to be clearly detectable in Sentinel-2 images. However, they occur very frequently clumped together, a common behaviour of this pioneer Genus. These flowering Pleroma dominated forest stands were easy to identify visually in the Sentinel-2 images because they combined several very distinctive features. First, an intense magenta-to-deep-purple color, which is an unusual color for other land covers in this ecological domain. Second, these identified Pleroma pixels should be undoubtedly identified as forested pixels and have a green color outside the blooming season. Third, Pleroma dominated forests often formed continuous magenta-to-deep-purple patches of size ranging from some 10 m × 10 m pixels to more than thousands of pixels and the shape of the patches tend to present linear features, likely representing the border of the space that was colonized by the Pleroma trees. Fourth, individual crowns were not visible, and the texture of the patches was very smooth during the blooming season with sometimes some inclusions of tree crowns of green color. Finally, texture of the Pleroma dominated forest stands outside the blooming season shown a smooth green texture, more homogeneous and with less shade than other forests. The first sample was constituted of images of background and of blooming Pleroma dominated forest stands that were following the previously described criteria. From this sample, we train a first model and applied it to the complete time series of Sentinel-2. From the results of this model, we obtained a first map of Pleroma trees and were also able to identify the main detection errors of this model, mainly clouds and dirt roads proximity with some unidentified agriculture fields or sometimes Eucalyptus plantation. The results of this first model were checked visually to produce a second larger sample (which was used for the results presented in this study) made up of images containing blooming Pleroma dominant to monospecific forest stands, a set of background images without blooming Pleroma and images identified erroneously by the first model as containing blooming Pleroma. While a large majority of the detected pixels were undoubtedly forest stands dominated by Pleroma trees, some other isolated trees of the genus Handroanthus (Ipê in Brazil or Lapacho in Argentina) with pink flowers and large crowns covering several pixels of 10 m × 10 m were also detected and kept in the training sample. For these particular Handroanthus trees, crowns were visible during and sometimes also outside the blooming season, which was not the case for detected Pleroma dominated forest stands. Finally, as our model detected also large Handroanthus trees, we must acknowledge that other tree species with highly similar features could also potentially being detected.The final training samples comprised a total of 158,612 images of 136 × 136 pixels. Among them, 35,541 contained blooming Pleroma trees and 123,071 images contained only background. Among the background images, there was nine different images types: images without blooming Pleroma, i.e., background such as other land covers, urban structures, water surfaces and agriculture and other land uses (57,007), images with forests containing Pleroma but outside the flowering period (23,427), images with clearly identified detection errors mainly located in the east of the São Paulo state (12,965), clouds and detection errors in clouds (10,991), images clearly identified as detection errors near the State of Bahia (9030), other detection errors over Atlantic Forest (5843), images of forests without Pleroma trees during the season of blooming (2170), images with identified detection errors in the northern part of Atlantic Forest (1126) and images with no data (512). Of these images, 80% (126,890) were used for training and 20% (31,722) used for validation.During network training, we used a standard stochastic gradient descent optimization, a binary cross-entropy loss and the optimizer RMSprop60 with a learning rate of 1e-4. We used the accuracy (i.e. the frequency with which the prediction matches the observed value) as the metrics of the model. However, due to the imbalance between the number of blooming Pleroma and background images, the metric of the model was weighted by one for the background and, for the Pleroma, by the ratio between the number of background images and the number of images containing blooming Pleroma: that is, ~3.5. The network was trained for 5000 epochs, where each epoch was made of 61 batches with 2048 images per batch and the model with the best weighted accuracy was kept for prediction (epoch 4331 and weighted accuracy of 99.58%). The training of the models took approximately 9 hours using a Nvidia RTX2080 Graphics Processing Unit (GPU) with an 8 GB memory.PredictionTo avoid border effects, each 10240 × 10240 pixels Sentinel-2 image was cropped on a regular grid of 128 × 128 pixels (1280 × 1280 m), and 4 neighboring pixels were added on each side to create an overlap between the patches. The function gdal_retile61 was used for this operation. Prediction was then made for each subset image: for each image, the detection model returned 0 or 1 if a blooming Pleroma was found in the image. Then the results were spatialized again using the grid, but this time, each cell of the grid only received 1 value, the prediction, resulting in a raster of 80 columns and 80 rows and a spatial resolution of 1280 m, of the same extent as the Sentinel-2 image. The value of the pixels (1 or 0) indicated the presence or absence of blooming Pleroma trees in this squared area of 1280 m of side. Prediction using GPU of a single tile of 10240 × 10240 pixels took approximately 1 minute on a Nvidia GTX1080 with an 8 GB memory and 45 s on a Nvidia RTX2080 with an 8 GB memory. The prediction for the complete Sentinel-2 time series presented in this work took approximately 22 days using a Nvidia GTX1080 GPU—from the 30 October 2020 to the 20 November 2020.Spatio-temporal analysisTo analyse the seasonality of the detections, daily maps of flowering Pleroma detections were produced for the studied period on a grid overlapping the entire Atlantic Forest (projection UTM 23S and datum WGS84) with a spatial resolution of 1280 m to match the resolution of the predictions. For each day, each pixel of the grid was given a classification: observed with flowering Pleroma, observed without flowering Pleroma, observed with clouds (using the cloud cover mask for Sentinel-2 images of this day) or as non-observed (no image or NA data for the pixel on that day). These daily grids were use to produce the map of flowering Pleroma trees (number of detections of flowering Pleroma for each pixel along the time series), the map of the total number of observations per pixel and the map of the total number of observations without clouds.To analyse the seasonality of blooming, the detection results were aggregated by month—even with a 5-day frequency there were still too few observations to analyse each annual timing and duration of flowering, and changes of the flowering dates between years were not expected based on the existing botanical information of the species. For each pixel, the number of detections per month were divided by the total number of observations without clouds per month. This enabled to normalized the detection values between zero and one and made sense given that we were not interest in the number of detections but rather in the times of the year when the number of detections was the highest: the peak of the blooming.To find the characteristics of these time series—one or more blooming peaks and the days when the blooming begins, peaks and stops—the normalized time series of mean monthly observations of flowering Pleroma were filtered using the Fourier transform (FT) (Eq. 1). This decomposition was made the keeping only the annual, bis- and tris-annual frequencies that compose the blooming signal, and to provide a continuous representation of the discrete blooming observations. In other worlds, the Fourier transform of the normalized time series observations enabled to model and compute the values of blooming for each day of the year and better estimate the days blooming started, peaked, and ended. While initially, a decomposition with only annual and biannual frequencies was expected to fit well to the times series (as more than two peaks per year were not expected), it appeared that when the two peaks were close in time (such as in a 2–3 month interval), only annual and biannual frequencies were not sufficient to give a good model of the signal, and the triannual frequency was added to resolve this issue. Furthermore, it was assumed that other periods in the signal were only constituted by noise.The blooming signal was modelled by the following equation:$$begin{aligned} {widehat{bloom}}(t)& = bloom_0 + pow_0 ,left( p_{4} sin left( 2pi frac{1}{4} t + rho _4right) right. \ & quadleft. + p_{6} sin left( 2pi frac{1}{6} t + rho _6right) + p_{12} sin left( 2pi frac{1}{12} t + rho _{12}right) right) end{aligned}$$
    (1)
    with (p_4 + p_6 + p_{12}=1) and for (t=1,ldots ,12 times n), ({widehat{bloom}}) is the filtered blooming time series; (bloom_0) is as an estimate of the mean annual blooming; t is the time in months; (rho _4), (rho _6) and (rho _{12}) are the delay of signal components with periods of 4 months, 6 months and 12 months, respectively; (pow_0) is the power of the signal and (p_4), (p_6) and (p_{12}) are the relative proportions in the signal of the periods of 4 months, 6 months and 12 months, respectively.To ease optimisation and cohere with the biological significance of our model, some data cleaning and adjustments were made. First, pixels with less than 4 observations over the 4-year period were removed from the analysis. Second, isolated peaks with only 1 or 2 observations during the 4-year period and between months without Pleroma detection were set to 0. Third, all the values of the normalized blooming time series were multiplied by 10, which seemed to ease the convergence of the optimisation algorithm. Fourth, the first months before and after the blooming period were set to a negative value equal to − 0.15 × the maximum value of the pixel time series. This was made based on the assumption that blooming is quite fast (based on the observation data) and happens between the month when the blooming is first observed and the previous month (and when blooming is last observed and the next month), and it forced the model to go below the 0 value during this period. Fifth, a weight was added to each point corresponding to its value, as we were interested in estimating accurately the peak value. A weight value of 0 was set to the month with a 0 value, and a weight of 1 was set to the months with negative blooming values (pre- and post-peak months). Finally, to facilitate the optimization, the time series of values and weights was replicated 3 times (n = 3). While this did not change the periodicity of the signal, it enabled to estimates better the value of the first and last month of the time series, as well as to ease optimisation. The parameters (bloom_0), (pow_0), (p_4), (p_6), (p_{12}), (rho _4), (rho _6) and (rho _{12}) were then estimated by a weighted least square minimization using the weights previously described. The accuracy of the model was given by the weighted R2 computed with the observed and predicted values of blooming for each month. As the Fourier transform is highly flexible, it can adjust almost perfectly to the data: the median of weighted R2 was close to 1 (with a 95% confidence interval—from percentile 2.75 to 97.5—of 0.86.0 to 1).Figure 14Examples of observed time series of detections in cloud-free images (%) and their daily estimation modeled using the Fourier transform.Full size imageAfter the decomposition of the blooming signal, a daily time series of 1 year of ({widehat{bloom}}) was computed with the obtained parameters (365 values) (Fig. 14). Daily values of months with a weight of 0 were set to 0 as well as predicted negative values. Then all peaks and pits were identified in the ({widehat{bloom}}) time series. A peak or pit is an observation that is preceded and followed by, respectively, lower or higher observations62,63. For each peak, the day of start and stop were identified using the pit values. After this analysis, we were able to describe the blooming time series: that is, if there were one or more peaks and, for each peak, the days when the blooming initiated, peaked and stopped.To determine if different populations could be identified based on flowering timing, a cluster analysis was performed. A classical K-means clustering analysis was made on a dataset containing, for each pixel where Pleroma were detected, the days of start, peak, and end of blooming, the associated normalized blooming values and the xy coordinates of the pixels. If there were two peaks for a pixels, a line for each peak was created in the dataset. As the xy coordinates were in metres, they carried most of the variance in the dataset. To avoid the artefact of having clusters based only on the distance between pixels, the xy coordinates were divided by 100,000 and rounded to the nearest unit. Before the clustering analysis, all variable were scaled and centered. The number of clusters was determined based on the curve representing the total within-cluster sum of squares as a function of the number of clusters, and also to have the maximum number of clusters.To describe the association of Pleroma trees with environmental variables, we first reclassified each environmental variables into 10 classes according to the variable’s quantiles. Then a bootstrap procedure was applied. For the number N of Pleroma trees detections, N random point locations were sampled within the Atlantic Forest domain, and the value of each environmental variable at each point was extracted and stored. This operation was repeated 100 times. It enabled to compare the number of Pleroma trees in each quantile class with the mean and gave us a 95% confidence interval for the number of points obtained by random spatial sampling in each class. Using the elevation as an example, the null hypothesis of no spatial association between Pleroma trees and elevation was rejected at a level of 0.05% if the number of Pleroma trees in a quantile class of the elevation was outside the (0.025, 0.975) quantiles of the empirical distribution of elevation obtained by random location sampling in the same class. The same analysis of association with the environmental variables was made for the Handroanthus population identified by the K-means clustering analysis.All analyses were performed using R project software55. More

  • in

    A spatial analysis of lime resources and their potential for improving soil magnesium concentrations and pH in grassland areas of England and Wales

    Magnesium soil contents and relationships to livestock healthMagnesium (Mg) deficiency (hypomagnesaemia) in ruminant livestock is a serious issue for the agricultural sector and accounts for a significant number of animal deaths annually. It is caused by a diet deficient in Mg, or due to an imbalance in the supply of Mg in comparison to other mineral cations1. Hypomagnesaemia is likely to be responsible for lower-productivity and diminished well-being in more animals in a herd compared to those displaying acute symptoms, given that herds/flocks generally receive a common diet2,3. If Mg deficiency could be prevented, it would be of benefit to both animal welfare and economic productivity. Recent research has confirmed that whilst hypomagnesaemia is commonly reported by UK farmers, the reported use of preventative measures is low, and the use of pasture interventions is lower still4. Pasture interventions can include the application of Mg-rich fertiliser or lime products, or selection of sward species with a propensity to take-up elevated Mg concentrations4,5.One aspect of dietary supply is the geographic control of pasture and farm-produced fodder. It is known that total Mg and plant-available Mg concentrations in soil are controlled by geological and geographic factors6,7,8,9 and that there is little evidence for any changes in pasture soil Mg concentrations through time6,10. The magnesium content of soil relates to that of the bedrock, where it is high in the bedrock it is high in soil and vice versa. Thus, the composition of all pasture and farm-grown fodder will be influenced by this natural environmental endowment as well as pasture management decisions.Grass productivity and soil pHA key pasture management activity is that of soil pH—which here is reported as measured in water (pHw), consistent with standard agronomic laboratory practice in the UK11. Grassland mineral soil is recommended to be maintained at pHw ≥ 6.0 in Britain12,13. In Ireland, where grass-clover pasture is more widely practiced, a pHw threshold of 6.5 is recommended14. However, multiple lines of evidence exist that indicate pasture soil is frequently below these pHw recommendations. Private sector on-farm sample data summaries from the UK consistently show pH typically below recommendations in pasture soil: the most recent annual data synthesis reports 57% of grassland soil with pHw ≤ 5.99, and 27% with pHw 6.00–6.498. This is consistent with systematically collected public sector data across the north of Ireland, where 84% of pasture samples were below the clover-grass recommended threshold of pHw 6.515.Grassland production is widespread in Wales and western England (Fig. 1). Two environmental factors jointly contributing to the lower pH in these regions are (1) geological—these areas are most often on soils which are developed over rocks with low concentrations of base cations (Figs. 1 and 2); and, (2) these areas are also often upland areas, associated with typically higher rainfall16 which will further leach base cations. Added to these environmental factors are the application of nitrogen fertilisers which have an acidifying effect17. Thus, many pastoral areas require treatment using agricultural lime in order to optimize soil pH for grass growth18.The use of liming materialsThe opportunity to improve grazing livestock Mg nutrition through use of Mg-rich lime is identified in guidance available to farmers12. This can have the dual benefits of maintaining soil pH for grass growth and ensuring Mg levels in livestock feed is at sufficient levels19,20,21. The combination of soil treatment for pH and Mg would therefore appear to be an efficient solution to solve issues surrounding Mg deficiency22. Conversely, for many soils with existing high Mg levels it may be important to treat with low Mg liming materials to ensure an optimal Ca–Mg balance to preserve the soil structure23.The use of Mg lime is only one of many methods of controlling Mg levels in livestock feed. Other methods, such as direct additions to feeds, salt licks, pelletised fertilisers products are also effective in reducing incidences of hypomagnesaemia and need to be considered as part of holistic review of a individual farms requirements, this is discussed in Kumssa et al.5.Maintaining optimal soil pH will directly affect the productivity of grass used for grazing, and will increase fertiliser use efficiency14. However, in some cases, for example upland sheep farming, a low investment—low return approach, with minimum interventions such as liming, may be entirety sensible and appropriate to the farm business and local landscape24,25.The extent to which agricultural lime is used in Britain is captured through the annual British Survey of Fertiliser Practice (BSFP) and can also be inferred from commodity production statistics. Production can be regarded as a good proxy for consumption since due to its high bulk and low price it is not exported in significant quantities.The BSFP, is an annual Department for Environment, Food and Rural Affairs (DEFRA) survey26, which representatively samples fertiliser and lime use across the British farming sector. This captures information on lime use in three geological material categories as used in arable and pastoral systems, as well as use of sugar beet lime and ‘other’ options. Sugar beet lime use is very low on grassland (generally unrecorded on ‘permanent’ pasture); ‘other’ categories are generally on a par with Mg-lime, but more detailed liming characteristics are not reported. Figure 3 shows a clear trend in decreasing production of agricultural liming material over the last 40 years. Lime use in the UK peaked in the late 1950s and mid 1960s likely due to a subsidy for agricultural lime in place at the time, this ended in 1978, causing prices to increase and subsequently lime use to decrease27. The use of agricultural lime has continued on a declining, or flat trend, likely due to reluctance to engage in soil treatments that are seen to be costly and a lack of knowledge over its potential benefits.Figure 1Classes of land cover considered by this study, based on Land Cover Map 201533. Some features of this map are based on digital spatial data licensed from the UK Centre for Ecology and Hydrology. Created using ArcMap 10.7.1, ESRI, 2019.Full size imageFigure 4 shows the use of geological lime products to be low at present in respect of the proportion of fields to which lime is reported to be applied, and that this is particularly pervasive on permanent grassland, with a 10-year average to 2019 of 2.9% (range 2.0–4.1%). Of this, a 10-year average of 0.4% of fields had Mg-lime applied, with 1.8% of fields having limestone applied. The limited use of chalk (0.1% of fields) probably reflects the distance between the majority of pasture and the outcrop of the chalk. Recent grassland ( More

  • in

    Deforestation is the turning point for the spreading of a weedy epiphyte: an IBM approach

    1.de Wet, J. M. J. & Harlan, J. R. Weeds and domesticates: Evolution in the man-made habitat. Econ. Bot. 29(2), 99–108. https://doi.org/10.1007/BF02863309 (1975).Article 

    Google Scholar 
    2.Ceballos, G. et al. Accelerated modern human—Induced species losses: Entering the sixth mass extinction. Sci. Adv. 1(June), 1–6. https://doi.org/10.1126/sciadv.1400253 (2015).Article 

    Google Scholar 
    3.Wilcove, D. S. Nest predation in forest tracts and the decline of migratory songbirds. Ecology 66(4), 1211–1214 (1985).Article 

    Google Scholar 
    4.Airoldi, L. & Bulleri, F. Anthropogenic disturbance can determine the magnitude of opportunistic species responses on marine urban infrastructures. PLoS ONE https://doi.org/10.1371/journal.pone.0022985 (2011).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    5.Baker, H. G. The evolution of weeds. Annu. Rev. Ecol. Syst. 5, 1–24. https://doi.org/10.2307/2096877 (1974).ADS 
    Article 

    Google Scholar 
    6.Richardson, D. M. et al. Naturalization and invasion of alien plants: Concepts and definitions. Divers. Distrib. 6, 93–107 (2008).Article 

    Google Scholar 
    7.van Etten, M. L., Conner, J. K., Chang, S. M. & Baucom, R. S. Not all weeds are created equal: A database approach uncovers differences in the sexual system of native and introduced weeds. Ecol. Evol. 7(8), 2636–2642. https://doi.org/10.1002/ece3.2820 (2017).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    8.Booth, B. D. & Swanton, C. J. Assembly theory applied to weed communities 50th Anniversary—Invited Article Assembly theory applied to weed communities. Weed Sci. 50(3), 2–13. https://doi.org/10.1614/0043-1745(2002)050 (2002).CAS 
    Article 

    Google Scholar 
    9.Kuester, A., Conner, J. K., Culley, T. & Baucom, R. S. How weeds emerge: A taxonomic and trait-based examination using United States data. New Phytol. 202(3), 1055–1068. https://doi.org/10.1111/nph.12698 (2014).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    10.van Kleunen, M. et al. The ecology and evolution of alien plants. Annu. Rev. Ecol. Evol. Syst. https://doi.org/10.1146/annurev-ecolsys-110617-062654 (2018).Article 

    Google Scholar 
    11.de Bona, S. et al. Spatio-temporal dynamics of density-dependent dispersal during a population colonisation. Ecol. Lett. 22, 634–644 (2019).PubMed 
    Article 

    Google Scholar 
    12.Baker, H. G. Self-compatibility and establishment after “long-distance” dispersal. Evolution 9(3), 347. https://doi.org/10.2307/2405656 (1955).Article 

    Google Scholar 
    13.Razanajatovo, M. et al. Plants capable of selfing are more likely to become naturalized. Nat. Commun. 7, 13313. https://doi.org/10.1038/ncomms13313 (2016).ADS 
    CAS 
    Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    14.Vallejo-Marín, M., Dorken, M. E. & Barrett, S. C. H. The ecological and evolutionary consequences of clonality for plant mating. Annu. Rev. Ecol. Evol. Syst. 41(1), 193–213. https://doi.org/10.1146/annurev.ecolsys.110308.120258 (2010).Article 

    Google Scholar 
    15.Rodger, J. G., Van Kleunen, M. & Johnson, S. D. Pollinators, mates and Allee effects: The importance of self-pollination for fecundity in an invasive lily. Funct. Ecol. 27(4), 1023–1033. https://doi.org/10.1111/1365-2435.12093 (2013).Article 

    Google Scholar 
    16.Barrett, S. C. H. & Harder, L. D. The ecology of mating and its evolutionary consequences in seed plants. Annu. Rev. Ecol. Evol. Syst. https://doi.org/10.1146/annurev-ecolsys-110316-023021 (2017).Article 

    Google Scholar 
    17.Klimeš, L., Klimešová, J., Hendriks, R. & van Groenendael, J. Clonal plant architecture: A comparative analysis of form and function. In The Ecology and Evolution of Clonal Plants (eds De Kroon, H. & Van Groenendael, J. M.) 1–29 (Backhuys, 1997).
    Google Scholar 
    18.Barrett, S. C. H. Influences of clonality on plant sexual reproduction. Proc. Natl. Acad. Sci. 112(29), 8859–8866. https://doi.org/10.1073/pnas.1501712112 (2015).ADS 
    CAS 
    Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    19.Heywood, J. S. Spatial analysis of genetic variation in plant populations. Annu. Rev. Ecol. Syst. 22, 335–355 (1991).Article 

    Google Scholar 
    20.Barrett, S. C. H. Evolution of mating systems: Outcrossing versus selfing. In The Princeton Guide to Evolution (ed. Losos, J. B.) 356–362 (Princeton University Press, 2013).
    Google Scholar 
    21.Barrett, S. C. H., Arunkumar, R. & Wright, S. I. The demography and population genomics of evolutionary transitions to self-fertilization in plants. Philos. Trans. R. Soc. B Biol. Sci. 369(1648), 20130344 (2014).Article 

    Google Scholar 
    22.Picó, F. X., Quintana-Ascencio, P. F., Mildén, M., Ehrlén, J. & Pfingsten, I. Modelling the effects of genetics and habitat on the demography of a grassland herb. Basic Appl. Ecol. 10(2), 122–130. https://doi.org/10.1016/j.baae.2008.02.006 (2009).Article 

    Google Scholar 
    23.Ellstrand, N. C. & Roose, M. L. Patterns of genotypic diversity in clonal plant species. Am. J. Bot. 74, 123–131 (1987).Article 

    Google Scholar 
    24.Loh, R., Scarano, F. R., Alves-Ferreira, M. & Salgueiro, F. Implications of clonality to population genetic structure of the nurse species Aechmea nuducaulis (L.) Griseb. (Bromeliaceae). Bot. J. Linn. Soc. 178, 329–341 (2015).Article 

    Google Scholar 
    25.Hedrick, P. W. Purging inbreeding depression and the probability of extinction: Full-sib mating. Heredity 73, 363–372. https://doi.org/10.1038/hdy.1994.183 (1994).Article 
    PubMed 

    Google Scholar 
    26.Arunkumar, R., Ness, R. W., Wright, S. I. & Barrett, S. C. H. The evolution of selfing is accompanied by reduced efficacy of selection and purging of deleterious mutations. Genetics 199(3), 817–829. https://doi.org/10.1534/genetics.114.172809 (2015).Article 
    PubMed 

    Google Scholar 
    27.Pannell, J. R. & Barrett, S. C. H. Baker’s law revisited: Reproductive assurance in a metapopulation. Evolution 52(3), 657–668. https://doi.org/10.2307/2411261 (1998).Article 
    PubMed 

    Google Scholar 
    28.Hamrick, J. L. & Trapnell, D. W. Using population genetic analyses to understand seed dispersal patterns. Acta Oecologica 37, 641–649 (2011).ADS 
    Article 

    Google Scholar 
    29.Côrtes, M. C. et al. Low plant density enhances gene dispersal in the Amazonian understory herb Heliconia acuminata. Mol. Ecol. 22, 5716–5729 (2013).PubMed 
    Article 
    CAS 

    Google Scholar 
    30.Trapnell, D. W., Hamrick, J. L., Ishibashi, C. D. & Kartzinel, T. R. Genetic inference of epiphytic orchid colonization; it may only take one. Mol. Ecol. 22, 3680–3692. https://doi.org/10.1111/mec.12338 (2013).Article 
    PubMed 

    Google Scholar 
    31.Chung, M. Y. et al. Fine-scale genetic structure in populations of the spring ephemeral herb Megaleranthis saniculifolia (Ranunculaceae). Flora Morphol. Distrib. Funct. Ecol. Plants 240, 16–24 (2018).
    Google Scholar 
    32.Roberts, N. R., Dalton, P. J. & Jordan, G. J. Epiphytic ferns and bryophytes of Tasmanian tree-ferns: A comparison of diversity and composition between two host species. Austral Ecol. 30(2), 146–154. https://doi.org/10.1111/j.1442-9993.2005.01440.x (2005).Article 

    Google Scholar 
    33.Cardelús, C. L. & Chazdon, R. L. Inner-crown microenvironments of two emergent tree species in a lowland wet forest. Biotropica 37(2), 238–244. https://doi.org/10.1111/j.1744-7429.2005.00032.x (2005).Article 

    Google Scholar 
    34.Quaresma, A. C., Piedade, M. T. F., Wittmann, F. & ter Steege, H. Species richness, composition, and spatial distribution of vascular epiphytes in Amazonian black-water floodplain forests. Biodivers. Conserv. 27(8), 1981–2002. https://doi.org/10.1007/s10531-018-1520-3 (2018).Article 

    Google Scholar 
    35.Claver, F. K., Alaniz, J. R. & Caldíz, D. O. Tillandsia spp.: Epiphytic weeds of trees and bushes. For. Ecol. Manag. 6(4), 367–372. https://doi.org/10.1016/0378-1127(83)90044-0 (1983).Article 

    Google Scholar 
    36.Bartoli, C. G., Beltrano, J., Fernández, L. V. & Caldíz, D. O. Control of the epiphytic weeds Tillandsia recurvata and Tillandsia aeranthos with different herbicides. For. Ecol. Manage. 59, 289–294 (1993).Article 

    Google Scholar 
    37.Flores-Palacios, A., García-Franco, J. G. & Capistrán-Barradas, A. Biomass, phorophyte specificity and distribution of Tillandsia recurvata in a tropical semi-desert environment (Chihuahuan Desert, Mexico). Plant Ecology and Evolution 148(1), 68–75 (2015).Article 

    Google Scholar 
    38.Birge, W. I. The anatomy and some biological aspects of the “ball moss”, Tillandsia recurvata, 1–24. L. Bull. Univ. Tex. 194(20) (1911).39.Smith, L. B. & Downs, R. J. Tillandsioideae (Bromeliaceae). In Flora Neotropica Monograph 14(2), 663–1492 (1977).40.Hewitt, G. M. (1996). Some genetic consequences of ice ages, and their role in speciation. Biological Journal of the Linnaean Society, 58(July), 247–276. Retrieved from papers3://publication/uuid/B9DB7D5E-D6AE-404C-BFFC-9F813345329441.McWilliams, E. Chronology of the Natural Range Expansion of Tillandsia recurvata (Bromeliaceae) in Texas. Contributions to Botany 15(2), 343–346 (1992).42.Flores-Palacios, A., Barbosa-Duchateau, C. L., Valencia-Díaz, S., Capistrán-Barradas, A. & García-Franco, J. G. Direct and indirect effects of Tillandsia recurvata on Prosopis laevigata in the Chihuahua desert scrubland of San Luis Potosi, Mexico. J. Arid Environ. 104, 88–95. https://doi.org/10.1016/j.jaridenv.2014.02.010 (2014).ADS 
    Article 

    Google Scholar 
    43.Benzing, D. H. Bromeliaceae: Profile of an Adaptive Radiation (Cambridge University Press, 2000).Book 

    Google Scholar 
    44.Benzing, D. H. Air Plants: Epiphytes and Aerial Gardens (Cornell University Press, 2012).Book 

    Google Scholar 
    45.Foster, M. D. Blueprint of the jungle as depicted by the altitude of growth of the Bromeliadswith notes on the culture of certain tropical epiphytes. Bull. N. Y. Bot. Garden 46, 9–16 (1945).
    Google Scholar 
    46.Soltis, D. E., Gilmartin, A. J., Rieseberg, L. & Gardner, S. Genetic variation in the epiphytes Tillandsia ionantha and T. recurvata (Bromeliaceae). Am. J. Bot. 74(4), 531–537 (1987).CAS 
    Article 

    Google Scholar 
    47.Orozco-Ibarrola, O. A., Flores-Hernández, P. S., Victoriano-Romero, E., Corona-López, A. M. & Flores-Palacios, A. Are breeding system and florivory associated with the abundance of Tillandsia species (Bromeliaceae)?. Bot. J. Linn. Soc. 177(1), 50–65. https://doi.org/10.1111/boj.12225 (2015).Article 

    Google Scholar 
    48.Chilpa-Galván, N. et al. Seed traits favouring dispersal and establishment of six epiphytic Tillandsia (Bromeliaceae) species. Seed Sci. Res. https://doi.org/10.1017/S0960258518000247 (2018).Article 

    Google Scholar 
    49.Southwood, T. & Kennedy, C. Trees as islands. Oikos 41(3), 359–371. https://doi.org/10.2307/3544094 (1983).Article 

    Google Scholar 
    50.Burns, K. C. Network properties of an epiphyte metacommunity. J. Ecol. 95(5), 1142–1151 (2007).Article 

    Google Scholar 
    51.Trapnell, D. W., Hamrick, J. L. & Nason, J. D. Three-dimensional fine-scale genetic structure of the neotropical epiphytic orchid, Laelia rubescens. Mol. Ecol. 13, 1111–1118 (2004).CAS 
    PubMed 
    Article 

    Google Scholar 
    52.Torres, E., Riofrío, M.-L. & Iriondo, J. M. Complex fine-scale spatial genetic structure in Epidendrum rhopalostele: an epiphytic orchid. Heredity https://doi.org/10.1038/s41437-018-0139-1 (2018).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    53.Victoriano-Romero, E., Valencia-Díaz, A., Toledo-Hernández, V. H. & Flores-Palacios, A. Dispersal limitation of Tillandsia species correlates with rain and host structure in a central Mexican tropical dry forest. PLoS ONE 12(2), e0171614 (2017).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    54.Martins, S. E. (2009). Flora fanerogâmica do estado de São Paulo. FAPESP: Instituto de Botânica.55.Chaves, C. J. N., Dyonisio, J. C. J. C. & Rossatto, D. R. D. R. Host trait combinations drive abundance and canopy distribution of atmospheric bromeliad assemblages. AoB Plants 8(October 2015), plw010. https://doi.org/10.1093/aobpla/plw010 (2016).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    56.Epps, C. W. & Keyghobadi, N. Landscape genetics in a changing world: Disentangling historical and contemporary influences and inferring change. Mol. Ecol. 24(24), 6021–6040. https://doi.org/10.1111/mec.13454 (2015).Article 
    PubMed 

    Google Scholar 
    57.Cushman, S. A., Shirk, A. & Landguth, E. L. Separating the effects of habitat area, fragmentation and matrix resistance on genetic differentiation in complex landscapes. Landsc. Ecol. 27(3), 369–380. https://doi.org/10.1007/s10980-011-9693-0 (2012).Article 

    Google Scholar 
    58.Jackson, N. D. & Fahrig, L. Habitat amount, not habitat configuration, best predicts population genetic structure in fragmented landscapes. Landsc. Ecol. 31(5), 951–968. https://doi.org/10.1007/s10980-015-0313-2 (2016).Article 

    Google Scholar 
    59.Grimm, V. & Railsback, S. F. Individual-Based Modelling and Ecology (Princeton University Press, 2005).MATH 
    Book 

    Google Scholar 
    60.Csilléry, K., Blum, M. G. B., Gaggiotti, O. E. & François, O. Approximate Bayesian Computation (ABC) in practice. Trends Ecol. Evol. 25(7), 410–418. https://doi.org/10.1016/j.tree.2010.04.001 (2010).Article 
    PubMed 

    Google Scholar 
    61.Udupa, S. M. & Baum, M. High mutation rate and mutational bias at (TAA)n microsatellite loci in chickpea (Cicer arietinum L.). Mol. Genet. Genom. 265(6), 1097–1103. https://doi.org/10.1007/s004380100508 (2001).CAS 
    Article 

    Google Scholar 
    62.Anmarkrud, J. A., Kleven, O., Bachmann, L. & Lifjeld, J. T. Microsatellite evolution: Mutations, sequence variation, and homoplasy in the hypervariable avian microsatellite locus HrU10. BMC Evol. Biol. 8(1), 1–10. https://doi.org/10.1186/1471-2148-8-138 (2008).CAS 
    Article 

    Google Scholar 
    63.Marriage, T. N. et al. Direct estimation of the mutation rate at dinucleotide microsatellite loci in Arabidopsis thaliana (Brassicaceae). Heredity 103(4), 310–317. https://doi.org/10.1038/hdy.2009.67 (2009).CAS 
    Article 
    PubMed 

    Google Scholar 
    64.Bernal, R., Valverde, T. & Hernández-Rosas, L. Habitat preference of the epiphyte Tillandsia recurvata (Bromeliaceae) in a semi-desert environment in Central Mexico. Can. J. Bot. 83(10), 1238–1247 (2005).Article 

    Google Scholar 
    65.Chaves, C. J. & Rossatto, D. R. Unravelling intricate interactions among atmospheric bromeliads with highly overlapping niches in seasonal systems. Plant Biol. 22(2), 243–251 (2020).CAS 
    PubMed 
    Article 

    Google Scholar 
    66.Vekemans, X. & Hardy, O. J. New insights from fine-scale spatial genetic structure analyses in plant populations. Mol. Ecol. 13(4), 921–935. https://doi.org/10.1046/j.1365-294X.2004.02076.x (2004).CAS 
    Article 
    PubMed 

    Google Scholar 
    67.Ward, S. Genetic analysis of invasive plant populations at different spatial scales. Biol. Invasions 8(3), 541–552. https://doi.org/10.1007/s10530-005-6443-8 (2006).Article 

    Google Scholar 
    68.Pettengill, J. B., Briscoe Runquist, R. D. & Moeller, D. A. Mating system divergence affects the distribution of sequence diversity within and among populations of recently diverged subspecies of Clarkia xantiana (Onagraceae). Am. J. Bot. 103(1), 99–109. https://doi.org/10.3732/ajb.1500147 (2016).CAS 
    Article 
    PubMed 

    Google Scholar 
    69.Atwater, D. Z., Fletcher, R. A., Dickinson, C. C., Paterson, A. H. & Barney, J. N. Evidence for fine-scale habitat specialization in an invasive weed. J. Plant Ecol. 11(2), 189–199. https://doi.org/10.1093/jpe/rtw124 (2018).Article 

    Google Scholar 
    70.Li, J. & Dong, M. Fine-scale clonal structure and diversity of invasive plant Mikania micrantha H.B.K. and its plant parasite Cuscuta campestris Yunker. Biol. Invasions 11(3), 687–695. https://doi.org/10.1007/s10530-008-9283-5 (2009).MathSciNet 
    Article 

    Google Scholar 
    71.Ren, M. X., Cafasso, D., Cozzolino, S. & Pinheiro, F. Extensive genetic differentiation at a small geographical scale: Reduced seed dispersal in a narrow endemic marsh orchid, Anacamptis robusta. Bot. J. Linn. Soc. 183(3), 429–438. https://doi.org/10.1093/botlinnean/bow017 (2017).Article 

    Google Scholar 
    72.Barluenga, M. et al. Fine-scale spatial genetic structure and gene dispersal in Silene latifolia. Heredity 106(1), 13–24. https://doi.org/10.1038/hdy.2010.38 (2011).CAS 
    Article 
    PubMed 

    Google Scholar 
    73.Charbonneau, A. et al. Weed evolution: Genetic differentiation among wild, weedy, and crop radish. Evol. Appl. https://doi.org/10.1111/eva.12699 (2018).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    74.Sagnard, F., Oddou-Muratorio, S., Pichot, C., Vendramin, G. G. & Fady, B. Effects of seed dispersal, adult tree and seedling density on the spatial genetic structure of regeneration at fine temporal and spatial scales. Tree Genet. Genomes 7(1), 37–48. https://doi.org/10.1007/s11295-010-0313-y (2011).Article 

    Google Scholar 
    75.Counsens, R. & Mortimer, M. Dynamics of Weed Populations (Cambridge University Press, 1995).Book 

    Google Scholar 
    76.Loreau, M. et al. Unifying sources and sinks in ecology and Earth sciences. Biol. Rev. 88, 365–379 (2013).PubMed 
    Article 

    Google Scholar 
    77.dos Santos, L. S. et al. Generalized Allee effect model. Theory Biosci. 133, 117–124 (2014).PubMed 

    Google Scholar 
    78.Spruch, L. et al. Modeling community assembly on growing habitat “islands”: A case study on trees and their vascular epiphyte communities. Theor. Ecol. 12, 1–17 (2019).Article 

    Google Scholar 
    79.Einzmann, H. J. R. & Zotz, G. “No signs of saturation”: long-term dynamics of vascular epiphyte communities in a human-modified landscape. Biodivers. Conserv. 26, 1393–1410 (2017).Article 

    Google Scholar 
    80.Belinchón, R., Harrison, P. J., Mair, L., Várkonyi, G. & Snäll, T. Local epiphyte establishment and future metapopulation dynamics in landscapes with different spatiotemporal properties. Ecology 98(3), 741–750. https://doi.org/10.1002/ecy.1686 (2017).Article 
    PubMed 

    Google Scholar 
    81.Vergara-Torres, C. A., Pacheco-Álvarez, M. C. & Flores-Palacios, A. Host preference and host limitation of vascular epiphytes in a tropical dry forest of central Mexico. J. Trop. Ecol. 26(6), 563–570. https://doi.org/10.1017/S0266467410000349 (2010).Article 

    Google Scholar 
    82.Barrett, S. C. H. & Kohn, J. R. Genetic and evolutionary consequences of small population size in plants: Implications for conservation. In Genetics and Conservation of Rare Plants (eds Falk, D. A. & Holsinge, K. E.) 3–30 (Oxford University Press, 1991).
    Google Scholar 
    83.Nathan, R., Horn, H. S., Chave, J. & Levin, S. A. Mechanistic models for tree seed dispersal by wind in dense forests and open landscapes. In Seed Dispersal and Frugivory-Ecologie, Evolution, Conservation 69–82 (2002). https://doi.org/10.1079/9780851995250.006984.Cousens, R. et al. Dispersal in Plants. A Population Perspective (Oxford University Press, 2008).Book 

    Google Scholar 
    85.Snäll, T., Ehrlén, J. & Rydin, H. Colonization-extinction dynamics of an epiphyte metapopulation in a dynamic landscape. Ecology 86(1), 106–115 (2005).Article 

    Google Scholar 
    86.Ruiz-Cordova, J. P., Toledo-Hernández, V. H. & Flores-Palacios, A. The effect of substrate abundance in the vertical stratification of bromeliad epiphytes in a tropical dry forest (Mexico). Flora Morphol. Distrib. Funct. Ecol. Plants 209(8), 375–384. https://doi.org/10.1016/j.flora.2014.06.003 (2014).Article 

    Google Scholar 
    87.Flores-Palacios, A., Bustamante-Molina, A. B., Corona-López, A. M. & Valencia-Díaz, S. Seed number, germination and longevity in wild dry forest Tillandsia species of horticultural value. Scientia Hortic. 187, 72–79 (2015).Article 

    Google Scholar 
    88.Goodman, R., & Herold, M. (2014). Why maintaining tropical forests is essential and urgent for a stable climate. Center for Global Development Working Paper, (385).89.Seymour, F. & Busch, J. Why Forests? Why Now?: The Science, Economics, and Politics of Tropical Forests and Climate Change (Brookings Institution Press, 2016).
    Google Scholar 
    90.Stephenson, N. L. et al. Rate of tree carbon accumulation increases continuously with tree size. Nature 507(7490), 90–93 (2014).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    91.Tel-Zur, N., Abbo, S., Myslabodsky, D. & Mizrahi, Y. Modified CTAB procedure for DNA isolation from epiphytic cacti of genera Hylocereus and Selenicereus (Cactaceae). Plant Mol. Biol. Rep. 17, 249–254 (1999).CAS 
    Article 

    Google Scholar 
    92.Chaves, C. J. N., Aoki-Gonçalves, F., Leal, B. S. S., Rossatto, D. R. & Palma-Silva, C. Transferability of nuclear microsatellite markers to the atmospheric bromeliads Tillandsia recurvata and T. aeranthos (Bromeliaceae). Braz. J. Bot. 41, 931–935. https://doi.org/10.1007/s40415-018-0494-4 (2018).Article 

    Google Scholar 
    93.Keenan, K., Mcginnity, P., Cross, T. F., Crozier, W. W. & Prodöhl, P. A. DiveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors. Methods Ecol. Evol. https://doi.org/10.1111/2041-210X.12067 (2013).Article 

    Google Scholar 
    94.Slatkin, M. A measure of population subdivision based on microsatellite allele frequencies. Genetics 139, 457–462 (1995).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    95.Nei, M. Genetic distances between populations. Am. Nat. 106, 283–292 (1972).Article 

    Google Scholar 
    96.Edwards, A. W. F. Distance between populations on the basis of gene frequencies. Biometrics 27, 873–881 (1971).CAS 
    PubMed 
    Article 

    Google Scholar 
    97.Reynolds, J. B., Weir, B. S. & Cockerham, C. C. Estimation of the coancestry coefficient: Basis for a short-term genetic distance. Genetics 105, 767–779 (1983).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    98.Kamvar, Z. N., Tabima, J. F. & Grünwald, N. J. Poppr: An R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ https://doi.org/10.7717/peerj.281 (2014).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    99.Paradis, E. pegas: an R package for population genetics with an integrated-modular approach. Bioinformatics 26, 419–420 (2010).CAS 
    PubMed 
    Article 

    Google Scholar 
    100.Excoffier, L., Smouse, P. E. & Quattro, J. M. Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data. Genetics 131, 479–491. https://doi.org/10.1007/s00424-009-0730-7 (1992).CAS 
    Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    101.Loiselle, B. A., Sork, V. L., Nason, J. & Graham, C. Spatial genetic structure of a tropical understory shrub, Psychotria officinalis (Rubiaceae). Am. J. Bot. 82(11), 1420–1425 (1995).Article 

    Google Scholar 
    102.Bailleul, D., Stoeckel, S. & Arnaud-Haond, S. RClone: A package to identify MultiLocus Clonal Lineages and handle clonal data sets in r. Methods Ecol. Evol. 7(8), 966–970. https://doi.org/10.1111/2041-210X.12550 (2016).Article 

    Google Scholar 
    103.Harrison, S. et al. Beta diversity on geographic gradients in Britain. J. Anim. Ecol. 61(1), 151–158 (1992).Article 

    Google Scholar 
    104.Jost, L. Partitioning diversity into independent alpha and beta components. Ecology 88(10), 2427–2439. https://doi.org/10.1890/07-1861.1 (2007).Article 
    PubMed 

    Google Scholar 
    105.Charney, N. & Record, S. Vegetarian: Jost diversity measures for community data. https://cran.r-project.org/web/packages/vegetarian/index.html (2012). Accessed Jul 2018.106.Wilensky, U. NetLogo (Northwestern University, Center for Connected Learning and Computer-Based Modeling, 1999).
    Google Scholar 
    107.Grimm, V. et al. A standard protocol for describing individual-based and agent-based models. Ecol. Model. 198(1–2), 115–126. https://doi.org/10.1016/j.ecolmodel.2006.04.023 (2006).Article 

    Google Scholar 
    108.Grimm, V. et al. The ODD protocol: A review and first update. Ecol. Model. 221(23), 2760–2768. https://doi.org/10.1016/j.ecolmodel.2010.08.019 (2010).Article 

    Google Scholar 
    109.Kooijman, B. & Kooijman, S. A. L. M. Dynamic Energy Budget Theory for Metabolic Organisation (Cambridge University Press, 2010).
    Google Scholar 
    110.Sibly, R. M. et al. Representing the acquisition and use of energy by individuals in agent-based models of animal populations. Methods Ecol. Evol. 4(2), 151–161 (2013).Article 

    Google Scholar 
    111.Johnston, A. S. A., Hodson, M. E., Thorbek, P., Alvarez, T. & Sibly, R. M. An energy budget agent-based model of earthworm populations and its application to study the effects of pesticides. Ecol. Model. 280, 5–17 (2014).CAS 
    Article 

    Google Scholar 
    112.van der Vaart, E., Johnston, A. S. A. & Sibly, R. M. Predicting how many animals will be where: How to build, calibrate and evaluate individual-based models. Ecol. Model. 326, 113–123 (2016).Article 

    Google Scholar 
    113.Garza, J. C. & Williamson, E. G. Detection of reduction in population size using data from microsatellite loci. Mol. Ecol. 10, 305–318 (2001).CAS 
    PubMed 
    Article 

    Google Scholar 
    114.Excoffier, L., Laval, G. & Schneider, S. Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evol. Bioinform. Online 1, 47–50 (2005).CAS 
    Article 

    Google Scholar 
    115.Csilléry, K., François, O. & Blum, M. G. abc: An R package for approximate Bayesian computation (ABC). Methods Ecol. Evol. 3(3), 475–479 (2012).Article 

    Google Scholar 
    116.Pastur, G. M., Lencinas, M. V., Cellini, J. M. & Mundo, I. Diameter growth: Can live trees decrease?. Forestry 80(1), 83–88. https://doi.org/10.1093/forestry/cpl047 (2007).Article 

    Google Scholar  More

  • in

    Shared patterns in body size declines among crinoids during the Palaeozoic extinction events

    1.Smith, F. A. et al. Body size evolution across the Geozoic. Annu. Rev. Earth. Planet. Sci. 44, 523–553 (2016).ADS 
    CAS 
    Article 

    Google Scholar 
    2.Sallan, L. & Galimberti, A. K. Body-size reduction in vertebrates following the end-Devonian mass extinction. Science 350, 812–815 (2015).ADS 
    CAS 
    Article 

    Google Scholar 
    3.Heim, N. A., Knope, M. L., Schaal, E. K., Wang, S. C. & Payne, J. L. Cope’srule in the evolution of marine animals. Science 347, 867–870 (2015).ADS 
    CAS 
    Article 

    Google Scholar 
    4.Kammer, T. W. & Ausich, W. I. The, “Age of Crinoids”: A Mississippian biodiversity spike coincident with wide spread carbonate ramps. Palaios 21, 238–248 (2006).ADS 
    Article 

    Google Scholar 
    5.Wright, D. F. Phenotypic innovation and adaptive constraints in the evolutionary radiation of Palaeozoic crinoids. Sci. Rep. 7, 13745 (2017).ADS 
    Article 

    Google Scholar 
    6.Segessenman, D. C. & Kammer, T. W. Testing reduced evolutionary rates during the Late Palaeozoic Ice Age using the crinoid fossil record. Lethaia 51, 330–343 (2018).Article 

    Google Scholar 
    7.Cole, S. R. & Hopkins, M. J. Selectivity and the effect of mass extinctions on disparity and functional ecology. Sci. Adv. 7, eabf4072 (2021).ADS 
    Article 

    Google Scholar 
    8.Baumiller, T. K. Echinoderms Through Time (Echinoderms Dijon) (eds. David, B., Guille, A., Féral, J. P. & Roux, M.) 193–198 (Balkema, 1994).9.Foote, M. Ecological controls on the evolutionary recovery of post-Paleozoic crinoids. Science 274, 1492–1495 (1996).ADS 
    CAS 
    Article 

    Google Scholar 
    10.Sallan, L. C., Kammer, T. W., Ausich, W. I. & Cook, L. A. Persistent predator-prey dynamics revealed by mass extinction. Proc. Natl. Acad. Sci. U.S.A. 108, 8335–8338 (2011).ADS 
    CAS 
    Article 

    Google Scholar 
    11.Ausich, W. I. & Kammer, T. W. Mississippian crinoid biodiversity, biogeography and macroevolution. Palaeontology 56, 727–740 (2013).Article 

    Google Scholar 
    12.Brom, K. R., Salamon, M. A. & Gorzelak, P. Body-size increase in crinoids following the end-Devonian mass extinction. Sci. Rep. 8, 9606 (2018).ADS 
    Article 

    Google Scholar 
    13.Borths, M. R. & Ausich, W. I. Ordovician-Silurian Lilliput crinoids during the end-Ordovician biotic crisis. Swiss J. Palaeontol. 130, 7–18 (2011).Article 

    Google Scholar 
    14.Payne, J. L., Jost, A. B., Wang, S. C. & Skotheim, J. M. A shift in the long-term mode of foraminiferan size evolution caused by the end-Permian mass extinction. Evolution 67, 816–827 (2013).CAS 
    Article 

    Google Scholar 
    15.Finnegan, S., Heim, N. A., Peters, S. E. & Fischer, W. W. Climate change and the selective signature of the Late Ordovician mass extinction. Proc. Natl. Acad. Sci. U.S.A. 109, 6829–6834 (2012).ADS 
    CAS 
    Article 

    Google Scholar 
    16.Peters, S. E. & Ausich, W. I. A sampling-adjusted macroevolutionary history for Ordovician-Early Silurian crinoids. Paleobiology 43, 104–116 (2008).Article 

    Google Scholar 
    17.Ausich, W. I. & Deline, B. Macroevolutionary transition in crinoids following the Late Ordovician extinction event (Ordovician–Early Silurian). Palaeogeogr. Palaeoclimatol. Palaeoecol. 361, 38–48 (2012).Article 

    Google Scholar 
    18.Huttenlocker, A. K. Body size reductions in non mammalian eutheriodont therapsids (synapsida) during the end-permian mass extinction. PLoS ONE 9, e87553 (2014).ADS 
    Article 

    Google Scholar 
    19.Urbanek, A. Biotic crises in the history of upper Silurian graptoloids: A palaeobiological model. Hist. Biol. 7, 29–50 (1993).Article 

    Google Scholar 
    20.Alroy, J. et al. Phanerozoic trends in the global diversity of marine invertebrates. Science 321, 97–100 (2008).ADS 
    CAS 
    Article 

    Google Scholar 
    21.Stigall, A. L. Speciation collapse and invasive species dynamics during the Late Devonian “Mass Extinction”. GSA Today 22, 4–9 (2012).Article 

    Google Scholar 
    22.Hone, D. W. E., Keesey, T. M., Pisani, D. & Purvis, A. Macroevolutionary trends in the Dinosauria: Cope’srule. J. Evol. Biol. 18, 587–595 (2005).CAS 
    Article 

    Google Scholar 
    23.Hunt, G. & Roy, K. Fittings and comparing models of phyletic evolution: Random walks and beyond. Paleobiology 32, 578–601 (2006).Article 

    Google Scholar 
    24.Hunt, G., Hopkins, M. J. & Lidgard, S. Simple versus complex models of trait evolution and stasis as are sponse to environmental change. Proc. Natl. Acad. Sci. U.S.A. 112, 4885–4890 (2015).ADS 
    CAS 
    Article 

    Google Scholar 
    25.Syverson, V. J. & Baumiller, T. K. Temporal trends of predation resistance in Paleozoic crinoid arm branching morphologies. Paleobiology 40, 417–427 (2014).Article 

    Google Scholar 
    26.Sheridan, J. A. & Bickford, D. Shrinking body size as an ecological response to climate change. Nat. Clim. Change. 1, 401–406 (2011).ADS 
    Article 

    Google Scholar 
    27.Ebert, T. A. Negative growth and longevity in the purple sea urchin Strongylocentrotus purpuratus (Stimpson). Science 157, 557–558 (1967).ADS 
    CAS 
    Article 

    Google Scholar 
    28.Sato, K. N. et al. Response of sea urchin fitness traits to environmental gradients across the southern California oxygen minimum zone. Front. Mar. Sci. 5, 258 (2018).Article 

    Google Scholar 
    29.Brom, K. R. Body-size trends of cyrtocrinids (Crinoidea, Cyrtocrinida). Ann. Paleontol. 105, 109–118 (2019).ADS 
    Article 

    Google Scholar 
    30.Webster, G. D. & Webster, D. W. Bibliography and Index of Paleozoic Crinoids, Coronates, and Hemistreptocrinoids. 1758–2012. http://crinoids.azurewebsites.net (2014)31.Hunt, G., & Paleo, T.S. Analyze Paleontological Time Series. R Package Version 0.5.2. http://cran.r-project.org/web/packages/paleoTS/ (2019).32.RStudio Team. RStudio: Integrated Development for R. (RStudio, PBC, 2020). http://www.rstudio.com/.33.Bergman, N. M., Lenton, T. M. & Watson, A. J. COPSE: A new model of biogeochemical cycling over Phanerozoic time. Am. J. Sci. 304, 397–437 (2004).ADS 
    CAS 
    Article 

    Google Scholar 
    34.Royer, D. L., Berner, R. A., Montañez, I. P., Tabor, N. J. & Beerling, D. J. CO2 as a primary driver of Phanerozoic climate. GSA Today 14, 4–10 (2004).Article 

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

    Google Scholar 
    36.Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D. & R Core Team. Linear and Nonlinear MixedEffectsModels. R Package Version 3.1-143. http://cran.r-project.org/package=nlme (2019).37.Sookias, R. B., Benson, R. B. & Butler, R. J. Biology, not environment, drives major patterns in maximum tetrapod body size through time. Biol. Lett. 8, 674–677 (2012).Article 

    Google Scholar 
    38.Rego, B. L., Wang, S. C., Altiner, D. & Payne, J. L. Within- and among-genus components of size evolution during mass extinction, recovery, and background intervals: A case study of Late Permian through Late Triassic foraminifera. Paleobiology 38, 627–643 (2012).Article 

    Google Scholar 
    39.Zhang, Z., Augustin, M. & Payne, J. L. Phanerozoic trends in brachiopod body-size from synoptic data. Paleobiology 41, 491–501 (2015).Article 

    Google Scholar  More