Interaction network structure reflects the ecological mechanisms acting within biological communities, which are affected by environmental conditions. In tropical forests, higher precipitation usually increases fruit production, which may lead frugivores to increase specialization, resulting in more modular and less nested animal-plant networks. In these ecosystems, El Niño is a major driver of precipitation, however, we still lack knowledge of how species interactions change under this influence. To understand bat-plant network structure during an extreme ENSO event, we determined the links between frugivorous bat species and the plants they consume by DNA barcoding seeds and pulp in bat faeces. These interactions were recorded in the dry forest and rainforest of Costa Rica, during the dry and the wet seasons of an extreme El Niño year. From these we constructed seasonal and whole-year bat-plant networks and analyzed their structures and dissimilarities. In general, networks had low nestedness, high modularity, and were dominated by one large compartment which included most species and interactions. Contrary to our expectations, networks were less nested and more modular in drier conditions, both in the comparison between forest types and between seasons. We suggest that increased competition, when resources are scarce during drier seasons and habitats, lead to higher resource partitioning among bats and thus higher modularity. Moreover, we have found similar network structures between dry and rainforests during El Niño and non El Niño years. Finally, most interaction dissimilarity among networks occurred due to interaction rewiring among species, potentially driven by seasonal changes in resource availability.

Materials and Methods

Study sites

Fieldwork was conducted at two forest sites in Costa Rica that show contrasting seasonality and precipitation: the Atlantic rainforest of La Selva Biological Station (10°25′19” N, 84°00′54” W) and the Pacific dry forest at Sector Santa Rosa of Área de Conservación Guanacaste (ACG) (10°48’53” N, 85˚36’54” W) (Figure 1). La Selva Biological Station covers 1,611 ha of lowland wet tropical forest between 35 to 137 m on the Caribbean slope of the Cordillera Central mountain range. It has a mean annual temperature of 25˚ C with a mean annual precipitation of 3,962 mm (Sigel, Sherry, & Young, 2006). La Selva Biological Station has a mild dry season that ranges from January to April and a wet season that goes from May until December (Sanford et al. 1994). Sector Santa Rosa (of ACG) covers >38,000 ha of tropical dry forest ranging from 0 m to 300 m, and is part of Área de Conservación Guanacaste (Asensio, Schaffner, & Aureli, 2015). Sector Santa Rosa (of ACG) has a mean annual temperature of 25˚C with a mean annual precipitation of 1,575 mm. The dry season ranges from December to April, while the wet season extends from May to November, where 85%-97% of the precipitation falls (Castro et al. 2018) (Figure 1). The annual precipitation is higher at La Selva Biological Station (range 2,809-6,164 mm) than at Área de Conservación Guanacaste (range 880-3,030 mm, six-month dry season) (Gillespie, Grijalva, & Farris, 2009).

Bat sampling

We captured bats using four to six mist nets (6m – 12m) opened along trails and near watercourses in the study area from 18h – 22h. In addition, a canopy net and harp trap were used in 2009 but these had low capture rates and so were not used in 2015. Sampling took place in the wet season May-July (Santa Rosa of ACG) (2009), and in the wet season July-August (Sector Santa Rosa of ACG) and September-October (La Selva) (2015), and in the dry season during January-February (Sector Santa Rosa of ACG) and March-April (La Selva) (2015). Sampling and bat identification during the non El Niño year was conducted as described in Clare et al. (2019). Sampling effort using mist nets was equal to approximately 2,250 m² hours within each season during the El Niño year, and approximately the same during the non El Niño year. We collected wing punches for another study and these also served to avoid recaptures in order to maintain the independence of the data in our analysis. We measured the forearm length with callipers (0.1 mm precision) and identified species following Reid (1997), Timm & LaVal (1998), and LaVal & Rodriguez-Herrera (2002). Bats were held in cloth bags for a maximum of two hours for the collection of faecal samples. All samples were frozen after collection (-20° C).

DNA extraction, PCR and sequencing

For this study, we focused on nectar and fruit eating species, which produced faecal samples consisting of either seeds or digested fruit pulp. For DNA extraction, PCR and sequencing of the samples we followed standard protocols for plants and all work was conducted by the Canadian Centre for DNA barcoding (CCDB) following these procedures (Ivanova, Kuzmina, & Fazekas, 2011). In brief, dried plant material from faeces (fruit pulp or seed) was placed in a sterile strip-tube with pre-aliquoted sterile stainless steel beads and the tissue was ground using a Tissue Lyser (Qiagen, USA). The ground material was incubated with 2x CTAB buffer at 65°C for 1 hour and DNA extraction was performed using a semi-automated glass fiber filtration method (Ivanova, Fazekas, & Hebert, 2008; Fazekas et al. 2012). Following established methods, we amplified a 552 bp fragment of the 5’ end of the large subunit of RuBisCO (rbcL) and the ~350 bp second nuclear encoded internal transcribed spacer (ITS2) flanking by the partial 5.6S and 26S ribosomal genes. Sanger sequencing was performed using a ABI 3730xl capillary sequencer (Ivanova, DeWaard, Hajibabaei, & Hebert, 2005; Ivanova & Grainger, 2006; Kuzmina & Ivanova, 2011a; Kuzmina & Ivanova, 2011b; Fazekas et al. 2012). Although plant DNA barcoding yields lower species resolution compared to fungi and animals (Hollingsworth, Graham, & Little, 2011), generally it provides robust results for identification of vascular plants at the genus level (Kress et al. 2009; Parmentier et al. 2013; Braukmann, Kuzmina, Sills, Zakharov, & Hebert, 2017). For the samples from the non-El Niño year, plant specimens collected on the sites were identified using rbcL and matK and the supplementary non-coding plastid region trnH-psbA (see Clare et al. 2019 for full methods). 


Identification of plant DNA sequences from bat faecal samples

We initially filtered all sequences for quality and excluded low quality sequences where the PHRED score was <30 as indexed on the Barcode of Life Data Management System (BOLD) (Ratnasingham & Hebert, 2007). At least one genetic marker was recovered for each faecal sample, thus no samples were excluded at this stage. We compared the obtained rbcL and ITS2 sequences with the reference libraries of GenBank and BOLD using the BLAST algorithm with default search parameters (Altschul, Gish, Miller, Myers, & Limpman, 1990) in GenBank and the combined BLAST and Hidden Markov Model methods implemented by the BOLD server (Ratnasingham & Hebert, 2007). For each reference database (BOLD, GenBank), we assigned query sequences to taxa based on highest percentage similarity, and considered a threshold of ≥97% to be a reliable assignment (Lamb, Winsley, Piper, Freidrich, & Siciliano, 2016). When there was an agreement between species-level matches for both markers (rbcL and ITS2) in both databases, with at least one match >97%, we assigned a species name. In cases where the query matched with equal similarity to multiple taxa in the same genus, we assigned the taxon to the level of genus only, and similarly we used the same approach to assign query sequences to the level of the family. Where rbcL and ITS2 sequences matched different species from different genera, both at >97%, we concluded that two taxa were present in the sample and assigned to both genera. Query sequences that did not show significant similarity to a reference were excluded from the analysis, with only one faecal sample removed due to this procedure.

To corroborate species assignments, for each candidate genus match, we reconstructed a gene phylogeny in which we included our query sequences together with all available reference sequences from species of the same genus present in BOLD that are also known to occur in Costa Rica. Sequences from rbcL and ITS2 of each plant genus were aligned with ClustalW (Larkin et al. 2007) in BioEdit v7.2.5 (Hall, 1999). For each alignment we ran a model selection test to check which would be the best method to build the phylogenetic tree based on the lowest BIC value. We ran model selection and built the phylogenetic trees using MEGA 6.06 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013). These phylogenies (not shown) recovered paraphyletic groupings for some species, perhaps through a lack of reference material, and therefore such species assignments were considered unreliable. To address this impact of taxonomic resolution, we reduced all data to genus-level designations and repeated our analyses for both species and genera data to check for consistency of results (see Supporting Information).

The identification of plant DNA sequences from bat faecal samples during the non El Niño year relied on GenBank and BOLD, with the exception of the trnH-psbA region which was not searchable within BOLD (see Clare et al. 2019 for more details) for our purposes we used the assignments as given in Clare et al. (2019).

Network matrices

We compiled the inferred interactions into weighted matrices where each cell value represented the number of observed interactions between each bat-plant taxon pair. We considered one realized interaction when the DNA of a plant taxon was detected in the faeces of one individual bat. We constructed matrices for (1) each forest site in which we pooled data from both seasons during the El Niño year (‘La Selva’ and ‘Santa Rosa’), and (2) for each forest site in which we separated the data collected for dry and wet season during the non El Niño and El Niño years (‘wet’ versus ‘dry’ for each site). As species usually distribute their interactions unevenly among partners (Ings et al 2009), weighted analysis might uncover preferences, and even modules, that are invisible to binary analysis. Thus, we additionally applied weighted metrics to characterize network topologies. All statistical analysis and network drawings were performed using R, version 3.3.2 (R Development Core Team, 2017).