Aim: To assess the range size patterns of ferns and lycophytes along elevational gradients at different latitudes in an ecographical transition zone and search for predictors of range size from a set of environmental factors. Location: Mexico, from 15° to 23° N. Taxon: Ferns and lycophytes. Methods: All terrestrial and epiphytic species were recorded in 658 plots of 400 m2 along eight elevational gradients. To test if the range size within assemblages increases with elevation and latitude, we calculated the latitudinal range using the northern and southern limits of each species and averaged the latitudinal range of all species within assemblages weighted by their abundances. We related climatic factors and the changes with latitude and elevation with range size using linear mixed-effects models. Results: Species richness per plot increased with elevation up to about 1500-2000 m, with strong differences in overall species richness between transects and a reduction with increasing latitude. The mean weighted range size of species within assemblages declined with elevation, and increased with latitude, as predicted by theory. However, we also found marked differences between the Atlantic and Pacific slopes of Mexico, as well as low range size in humid regions. The best models described about 76-80 % of the variability in range size and included the seasonality in both temperature and precipitation, and annual cloud cover. Main conclusion: Latitudinal and elevational patterns of range size in fern assemblages are driven by an interplay of factors favouring wide-ranging species (higher latitudes with increasing temperature seasonality; dryer habitat conditions) and those favouring species with restricted ranges (higher elevations; humid habitat conditions), with additional variation introduced by the specific conditions of individual mountain ranges. Climatically stable, humid habitats apparently provide favourable conditions for small-ranged fern species, and should accordingly be given high priority in regional conservation planning. Methods Study area: The Mexican transition zone is the complex area where the Neotropical and Nearctic biotas overlap, and in a strict sense corresponds to the mountain highlands of Mexico, Guatemala, El Salvador and Nicaragua (Halffter & Morrone, 2017). We here present data from eight elevational gradients at a range of 0 m to 3500 m elevation at 15-23° latitude N on both the Pacific and Atlantic (Gulf of Mexico) sides of Mexico (Figure 1, Table S4 in appendices). Three transects have been considered in previous studies: Los Tuxtlas (Krömer et al., 2013; Acebey, Krömer, & Kessler, 2017), Perote (Carvajal-Hernández & Krömer, 2015; Carvajal-Hernández, Krömer, López‐Acosta, Gómez‐Díaz, & Kessler, 2017), and Oaxaca (Hernández-Rojas et al., 2018). Los Tuxtlas including abundances was not published before (“Los Tuxtlas a”). Both transects from Los Tuxtlas were combined for the majority of the analysis. Fern sampling: On each gradient, we sampled the fern assemblages at regular elevational intervals of 100-300 m (every 500 m at Perote), depending on accessibility. At each elevation, depending on the suitability of the slope, 4-8 plots of 20 x 20 m (400 m2) were sampled with a consistent, standardized methodology (Kessler & Bach, 1999; Karger et al., 2014). The plots were established in natural zonal forest, avoiding special structural features like canopy gaps, ridges, ravines, riparian areas, tree fall gaps, landslides, and other disturbed areas whenever possible, which all change microenvironmental conditions and have special fern assemblages. In each plot, all fern species and their abundances were recorded for terrestrial (soil, rocks, and dead wood) and for epiphytic substrates. Species with long creeping rhizomes were counted as patches. Epiphytes were sampled up to heights of 8 m with trimming poles and recorded at greater heights by using binoculars, climbing lower parts of trees, and searching recently fallen trees and branches within and adjacent to the plots (Gradstein, Nadkarni, Krömer, Holz, & Nöske, 2003; Sarmento Cabral et al., 2015).

Island biogeography explores how biodiversity in island ecosystems arises and is maintained. The topographical complexity of islands can drive speciation by providing a diversity of niches that promote adaptive radiation and speciation. However, recent studies have argued that phylogenetic niche conservatism, combined with topographical complexity and climate change, could also promote speciation if populations are episodically fragmented into climate refugia that enable allopatric speciation. Adaptive radiation and phylogenetic niche conservatism therefore both predict that topographical complexity should encourage speciation, but they differ strongly in their inferred mechanisms. Using genetic (mtDNA and SNP) and morphological data we show high species diversity (22 species) in an endemic clade of Fijian Homalictus bees, with most species restricted to highlands and frequently exhibiting narrow geographical ranges. Our results indicate that elevational niches have been conserved across most speciation events, contradicting expectations from an adaptive radiation model but concordant with phylogenetic niche conservatism. Climate cycles, topographical complexity, and niche conservatism could interact to shape island biodiversity. We argue that phylogenetic niche conservatism is an important driver of tropical bee biodiversity but that this phylogenetic inertia also leads to major extinction risks for tropical ectotherms under future warming climates. Methods Sample locations and collection methods. Collections throughout Fiji were made between 2010 and 2017 from multiple localities including the main islands of Viti Levu, Vanua Levu, Kadavu and Taveuni, as well as multiple small islands in the Lau group (SI Appendix, Fig. S5). Sampling of specimens at each location was not biased towards particular species because, for these very small bees, only H. achrostuscould be easily identified in the field due to its distinctive coloration; all other species required microscopy or DNA sequencing for species identification.

Samples were collected from 3 m to 1,324 m asl (highest elevation of Fiji) by sweep netting both native and introduced plants, and from nesting aggregations along roadsides. For each collection site, latitude, longitude and elevation were recorded using a Garmin 550 (Garmin Ltd., USA); latitude and longitude were then checked against satellite images (Google Earth) to confirm accuracy. Once collected, bees were immediately transferred into vials containing 98% ethanol. Vials were kept cool at ~5˚C and ethanol was replaced within a week of collection to reduce DNA degradation.

Maps of Fiji were produced in ArcMap (50)and a digital elevation model (DEM) of the archipelago was provided by Fiji Lands Information System (FLIS).

Sampling bias and elevational species richness. It was not possible to evenly sample bees across all geographical and elevational regions of Fiji because physical access to many regions was restricted by terrain and lack of roads. Access constraints could therefore affect sampling effort and this, in turn, could influence ability to recover true species richness in different elevational bands. Here, we quantize sampling effort as the number of DNA sequences obtained for different elevations, categorized into 200 m asl bands. Because specimens were only identified to species levels after DNA sequencing, the number of obtained sequences represents sampling effort. We examined whether this sampling effort may have influenced our estimates of species richness using multiple regression, where the number of detected species was the dependent variable and the number of sequences (sampling effort) and elevational band were the independent variables. The relative importance of sampling effort and elevation band for detected species richness can then be explored by regression bvalues and their statistical significance.

DNA extraction and sequencing. Tissue samples for DNA extraction were obtained by removing a single hind leg from each of the 764 specimens. For all samples obtained after 2014, DNA extraction and PCR amplification was completed at the South Australian Regional Facility for Molecular Ecology and Evolution (SARFMEE). DNA extraction and PCR amplification of COI prior to the 2014 samples was completed at the Canadian Centre for DNA Barcoding (CCDB) at the Biodiversity Institute of Ontario (31)and amplification used the universal primer pair LepF1 and LepR2 (31, 51). Extractions at SARFMEE followed protocols described by (52)with the subsequent DNA eluted into 75 µL of TLE buffer. PCR amplification of the 710 bp fragment of the DNA (COI) was completed using the primers LCO1490 (forward) and HCO2198 (reverse). The 25 µL PCR reactions comprised the following reagents: Sterile H2O (15.9 µL), MRT buffer (5 µL), 1 µL (5 µM) of LCO1490, 1 µL(5 µM) of HCO2198, Immolase Taq (0.1 µL) and DNA from specimen (2 µL). The thermocycling regime comprised of one cycle at 94˚C for 10 minutes, then five cycles at 94˚C for 60 seconds, 46˚C for 90 seconds, 72˚C for 75 seconds, followed by 35 cycles at 94˚C for 60 seconds, 51˚C for 90 seconds, 72˚C for 75 seconds, followed by 72˚C for 10 minutes and then 25˚C for 2 minutes.

Sequences were checked against the NCBI BLAST database to screen for non-target DNA. Forward and reverse sequences were aligned and chromatograms visually checked before creating final consensus sequences in Geneious version 10.2.2 (53). Initial alignments were trimmed to 630 bp to avoid any problems associated with missing data.

Phylogenetic, elevational and species analyses. The full COI alignment consisted of 630 bp for 764 specimens. Partition finder version 2 was employed using BIC and a greedy algorithm to find the best partition schemes and DNA substitution models from widely-used (i.e. MrBayes) models (54-56). The first and second codon positions were combined into a single partition with an HKY + I substitution model. A GTR substitution model was applied to third codon position. The BEAST file and parameters for phylogenetic analyses were set using BEAUti version 1.8.4 (57). Because of the small numbers of substitutions on each branch, a strict clock was used to avoid overparameterization. To infer changes in elevation across the tree we included elevation as a continuous trait using a strict or relaxed Brownian motion model (confirmed as adequate given our λ, κ, and δestimates;Table 1). Phylogenetic analyses were implemented in BEAST version 1.10(34)with 200 million iterations sampled every 50,000thiteration. Resulting log files were analyzed in Tracer version 1.6 (58)and a burnin of 2.5x107iterations was employed, which was always after stationarity had been achieved. Maximum clade credibility trees and posterior probability support values were obtained using TreeAnnotator Version 1.8.4 (57). Each run was performed four times for each analysis to ensure consistent results and stationarity. Post-burnin log and tree files for each run were then combined using LogCombiner version 2.5.2 (59)for the final analysis.

To infer the evolutionary mode and phylogenetic signal in the elevation data, we used BayesTraits version 3.0 (29). The tree-transformation models employed in BayesTraits assume that each terminal taxon is a species, hence we repeated the BEAST analysis using only one DNA sample from each species, and elevation data as either the median or minimum for all samples of that species. The (reduced) BEAST analysis used 100 million iterations, sampling every 50,000thiteration; stationarity and burnin was checked as above. The resulting consensus tree was run in BayesTraits using the median and the minimum elevational value for each terminal taxon to estimate λ (degree of phylogenetic signal), κ (degree of punctuated evolution), and δ (degree of early burst, adaptive radiation). The model of best fit for each estimate was chosen using Akaike’s Information Criterion with 100 bootstrap replicates in Tracer (60). Analyses in BayesTraits used 500 million iterations sampled every 50,000thiteration. Each run was performed four times for each model at each elevation to ensure consistent results. BayesTraits log files for each run were then combined using LogCombiner version 2.5.2 (59)for the final analysis.

We attempted to co-estimate phylogeny and elevational niche evolution, but these analyses repeatedly failed to converge. Thus, to infer elevational changes across the full phylogeny, we mapped elevation across all post-burnin trees sampled in the full COI analysis. This was done using BEAST, under a standard rate-constant Brownian motion model, as well as a rate-variable Brownian motion model, which assumes rates vary across branches according to an uncorrelated relaxed clock (35). Stationarity and burnin were confirmed as above. Both models gave very similar ancestral state reconstructions, but the latter model fitted better and is shown in Fig. 1.

Genetic analyses of bee clades were explored using Arlequin version 3.11 (61). For each species with multiple haplotypes and a sample size of more than 10 specimens we calculated haploytype diversity (h) and pairwise FSTvalues.

SNP quality filtering and analyses. The thorax and front legs were taken from 19 individuals from H. fijiensis, H. tuiwawae, H. ostridorsum, H. groomi and H.sp.S, respectively. To perform Restriction-site Associated DNA sequencing (RAD seq), the solid state method Diversity Arrays Technology (DArT) was used (62). The restriction enzymes used werea combination of PstI and HpaIIenzymes. Only female specimens were used to avoid the impact of male haploidy on SNP diversity. Post filtering, missing data was capped at 1.16%.

A total of 62,426 SNP loci were called across all species. Using the R package ‘DArTR’ version 1.0.5 low quality loci were removed at a threshold of 0.85% removing loci with 15% or more missing values (63), leading to

Aim: Environmental and spatial factors are broadly recognized as important predictors of beta diversity patterns. However, the scale at which beta diversity patterns are evaluated will affect the outcoming results. For example, studies at larger scales will usually find spatial processes as the main predictor of beta diversity patterns. In this study we evaluate how beta diversity patterns change when analyses are conducted at different scales by reducing the scale of analysis in a hierarchical manner.

Taxon: Chiroptera.

Location: Atlantic Forest biome.

Methods: Information on the occurrence of 59 bat species were obtained from the Atlantic Bats and Species Link database. We partitioned beta diversity into its two components (nestedness and turnover), and calculated these indexes hierarchically: the biome in its entirety (all ecoregions); between larger regions (north, central and south); and between ecoregions within each region. We performed a Generalized Dissimilarity Model (GDM) to identify and predict the turnover of bat species in the Atlantic Forest based on geo-climatic predictors. We obtained 19 geo-climatic data from AMBDATA, an environmental dataset based on different data sources commonly used in species distribution modeling.

Results: We found that turnover was the main component influencing a latitudinal gradient when the biome was analysed in its entirety. However, when the scale of the analysis was reduced, we found that species loss (nestedness component) had a large effect in determining beta diversity dissimilarity. We also found that nestedness was the main pattern explaining beta diversity dissimilarity along a longitudinal gradient.

Main conclusions: Beta diversity patterns changed with the scale of analysis, which indicates that bat species composition does not follow the same pattern throughout the Atlantic Forest. This corroborates the importance of analysing beta diversity patterns at different scales in order to understand how environmental dissimilarity across geographic space can influence species distribution patterns.

Methods Occurrence and geo-climatic data

    Our study was based on 2,626 bat occurrence data points for 525 sites (coordinates) within the Atlantic Forest. Our data came from two sources. We extracted 1,795 occurrence data points for 160 sites from Muylaert et al. (2017). This dataset was compiled by bat specialists who also reviewed the taxonomy of the species and the coordinates of sampling sites. We then used the species list provided in Muylaert et al. (2017) to search for other occurrence records in speciesLink (data downloaded from http://splink.cria.org.br/). We obtained 830 occurrence records of bat species for 364 sites. We reviewed the dataset obtained from speciesLink according to reliability of information regarding: i) coordinates and site correspondence (we used google maps to check if the coordinates were referring to the places indicated), ii) correct taxonomy (we excluded species with “sp”, “ssp”, “cf” and “aff”), and iii) voucher specimens (we only considered records with specimens that were deposited in a museum). We also included a single occurrence record of Natalus macrourus (Trajano, 1984) from Parque Estadual Turístico do Alto Ribeira (PETAR), which was not considered by Muylaert et al. (2017) or SpeciesLink. Occurrence records belonging to the bat families Molossidae, Vespertilionidae and Embalonuridae were not included in our study because they are seldomly captured in mist-nets (Nogueira, Pol &, Peracchi, 1999; Nogueira, Pol, Monteiro &, Peracchi, 2008), which was the predominant method used for sampling bat species represented in our data sources.

    We obtained geo-climatic data from AMBDATA (available at http://www.dpi.inpe.br/Ambdata/index.php). The AMBDATA is an environmental dataset systematized from different data sources and commonly used in species distribution modelling. It consists of 19 bioclimatic variables at 30 arc-sec resolution (approx. 1 km). These are: 1) annual mean temperature (ºC); 2) mean diurnal range (ºC); 3) isothermality (mean diurnal range divided by annual temperature range, and multiplied by 100); 4) temperature seasonality (standard deviation *100); 5) maximum temperature of warmest month (ºC); 6) minimum temperature of coldest month (ºC); 7) temperature annual range (ºC); 8) mean temperature of wettest quarter (ºC); 9) mean temperature of driest quarter (ºC); 10) mean temperature of warmest quarter (ºC); 11) mean temperature of coldest quarter (ºC); 12) annual precipitation (mm); 13) precipitation of wettest month (mm); 14) precipitation of driest month (mm); 15) precipitation seasonality (coefficient of variation); 16) precipitation of wettest quarter (mm); 17) precipitation of driest quarter (mm); 18) precipitation of warmest quarter (mm); and 19) precipitation of coldest quarter (mm). We also included three non-climatic environmental variables from AMBDATA: 1) tree cover at a 500 m resolution (percentage); 2) elevation (m) at 3 arc-sec horizontal resolution (about 90 m) and a vertical resolution of 1 m and lastly, 3) declivity (degrees) generated from the elevation grid.

Beta diversity and the Generalized Dissimilarity Model (GDM)

    There are various dissimilarity indices to measure changes in species composition between assemblages. We used the Sorensen index (βsØr) as implemented in ‘Betapart package’ – ‘R-project’ (Baselga & Orme, 2012). The input data table consists of the presence and absence of bat species for each study site (latitude and longitude). The package computes the total dissimilarity across all sites, and calculates turnover (Simpson’s index, βsim) and nestedness (the difference between the Sorensen and Simpson index, βsne) components. ‘Betapart’ returns cluster and dissimilarity matrices (between pairwise sites, and pairwise matrices of shared and non-shared species between sites) of turnover and nestedness.

First, we computed total beta diversity and its two components, nestedness and turnover, among the ten Atlantic Forest ecoregions proposed by Olson et al. (2001). Then we split the Atlantic Forest into three larger regions (southern, central and northern). Lastly, nestedness and turnover were calculated among the ecoregions making up each of the three regions. Each region was treated separately.

    We used a species presence data frame with the coordinates of the occurrence sites to perform Generalized Dissimilarity Modeling (GDM), which analyses spatial patterns of pairwise dissimilarity in species composition between sites, using a nonlinear regression matrix. GDM quantifies dissimilarity using the Soresen Index (total beta diversity), then associates the turnover component (βsim) with biological distance (predictor variables) between sites (Fitzpatrick et al., 2013). The GDM procedure was used to predict bat species turnover across the Atlantic Forest based on environmental data. We used the ‘R packageGDM’ (Fitzpatrick & Lisk, 2016) to fit a GDM with the 22 environmental variables and the geographical distance (decimal degrees) between occurrence sites. The latter was calculated using the option ‘geo=T’ in the function ‘gdm’ of the GDM package. We used the parameter ‘weightType= richness’ to weight sites relative to the number of species to minimize sampling bias. We chose not to exclude sites with few species (i.e. less than five) because in our database over 360 occurrence sites had five or fewer species, and 90 sites had 10 species or fewer. Less than 11 sites recorded 50% or more of the total number of species (59), so that excluding sites with few species would lead to a significant loss of data. In addition, the average number of species per site was equal to five, and sites with few species are evenly distributed throughout the biome. Therefore, maintaining all sites while correcting for species richness, even with a low number of species, does not weaken our model (see Fig. S1.2). Patterns of species turnover can be visualized on a raster with RGB colour standards; areas with similar colours contain similar assemblages. The GDM matrix regression used was I-spline with three basic functions, meaning that we used three points (the minimum) to form the I-spline curve (Fitzpatrick & Lisk, 2016). I-Splines can be visualized in a graph showing the relationship of predicted biological distance versus observed biological distance, providing an indication of how species composition changes along each environmental gradient (Fitzpatrick & Lisk, 2016). The selection of the best subset of predictors for our model followed Williams, Belbin, Austin, Stein, & Ferrier (2012): the initial model included all predictors; variables that contributed less than 2% to model explanation were iteratively removed. Variable removal was done on a stepwise basis beginning with the elimination of the variable that contributed the least to model explanation. Variables were reassessed regarding their importance and significance during each step of model reduction (i.e., backward elimination). Our model started with 23 predictor variables and ended with 11.

To enable reproduction of analyses in the linked journal article, we provide a table of habitats, woody linear features and elevation of all pollinator survey transect sections. This table can be linked to the publicly available data, based on square ID and section ID, to reproduce the analyses (except for one discontinued square with ID marked ‘NA’ in the publicly available data).

Journal article abstract:

1. Pollinating insects provide economic value by improving crop yield. They are also functionally and culturally important across ecosystems outside of cropland. To understand landscape-level drivers of pollinator declines, and guide policy and intervention to reverse declines, studies must cover (1) multiple insect and plant taxa and (2) a range of agricultural and semi-natural land uses. Furthermore, in an era of woodland restoration initiatives and rewilding ideologies, the contribution of woodland and woody linear features (WLFs; e.g. hedgerows) to pollinator abundance demands further investigation. 2. We demonstrate fine-scale analysis of high-quality, co-located measurements from a national environmental survey. We relate pollinator transect counts to ground-truth habitat and WLF maps across 300 1km squares in Wales, UK. We look at effects of habitat type, flower cover, WLF density and habitat diversity on summer abundance (July and August) of eight insect groups, representing three insect orders (Lepidoptera, Hymenoptera and Diptera). 3. Compared with improved grassland (the dominant habitat in Wales), pollinator abundance is consistently higher in cropland and woodland - especially broadleaved woodland. For solitary bees and two hoverfly groups, abundance is predicted to be at least 1.5× higher in woodland ecosystems than elsewhere. Furthermore, we estimate contributions of WLFs to abundance in agriculturally improved habitats to be up to 14% for honeybees and up to 21% for hoverflies. 4. The abundance of all insect groups increases with flower cover, which is a key mechanism through which woodland, cropland and grassland support pollinators. Importantly, we observe diminishing returns of increasing flower cover for abundance of non-Apis pollinator groups, expecting roughly twice the increase in abundance per % flower cover from 0-5%, as compared with 10-15%. However, the shape of the relationship was inverted for honeybees, which showed steeper increases in abundance at higher flower cover. 5. Synthesis & applications: We provide a holistic view of the drivers of pollinator abundance in Wales, in which flower cover, woodland, WLFs and cropland are critical. We propose a key role for woodland creation, hedge-laying and farmland heterogeneity within future land management incentive schemes. Finally, we suggest targeting of interventions to maximise benefits for non-Apis pollinators. Specifically, increasing floral provision in areas where existing flower cover is low – e.g. in flower-poor improved grasslands - could effectively increase pollinator abundance and diversity, while prioritising wild over managed species.

Methods Further information on habitat and landscape feature survey methodology is available in supporting documentation for the published datasets (Maskell et al., 2020a, 2020b). This dataset represents an intersection of spatial versions of the published datasets with locations of a set of 200m pollinator transect sections. The main linked paper (Alison et al. 2021) also provides lots of additional information.

Habitat, woody linear features (WLFs) and elevation variables were extracted for each transect section using ArcGIS Desktop 10.6 (ESRI, Redlands, California). To compare transect sections in different habitats, we classified the underlying broad habitat for 200m pollinator transect sections through intersection with habitat polygons. For each section, the broad habitat accounting for the greatest proportion of its length was assigned as the dominant habitat type. During this process we sought to maximise sample size, while avoiding unrepresentative classifications related to missing or ambiguous habitat data. As such, we made no dominant habitat assignment if (1) the dominant broad habitat accounted for <100m of the section, (2) a section had incomplete (<90%) overlap with habitat survey data or (3) the dominant broad habitat was recorded as “Mosaic”. We also calculated the Shannon index of habitat diversity for each transect section, taking into account the total number of broad habitats and the dominance among them. To relate flower cover and insect counts to habitat classifications at the finest possible scale and resolution, our habitat intersections were length-based, and not area- or buffer-based. However, a Euclidean buffer was necessary to extract the total density of WLFs (m ha-1) within 10m of each transect section, including managed WLFs (hedges), unmanaged WLFs (lines of trees), and forestry linear features (belts of trees or scrub). Outside of GMEP survey data, a 5m resolution raster of elevation was provided by Welsh Government (the Nextmap Britain DTM by Intermap Technologies). The elevation of each transect section was taken as the mean elevation of all vertices in the digitised transect section. The final modelled dataset included pollinator counts from 4,449 section-visits across 295 1km squares.

Maskell, L. C., Astbury, S., Burden, A., Emmett, B. A., Garbutt, A., Goodwin, A., Henrys, P., Jarvis, S., Norton, L. R., Owen, A., Sharps, K., Smart, S. M., Williams, B., Wood, C. M., & Wright, S. M. (2020a). Landscape and habitat area data from the Glastir Monitoring and Evaluation Programme, Wales 2013-2016. NERC Environmental Information Data Centre. (Dataset). https://doi.org/10.5285/82c63533-529e-47b9-8e78-51b27028cc7f

Maskell, L. C., Astbury, S., Burden, A., Emmett, B. A., Garbutt, A., Goodwin, A., Henrys, P., Jarvis, S., Norton, L. R., Owen, A., Sharps, K., Smart, S. M., Williams, B., Wood, C. M., & Wright, S. M. (2020b). Landscape linear feature data from the Glastir Monitoring and Evaluation Programme, Wales 2013-2016. NERC Environmental Information Data Centre. (Dataset). https://doi.org/10.5285/f481c6bf-5774-4df8-8776-c4d7bf059d40