Latitudinal patterns in phylogenetic and functional diversity of Diptera in temperate bogs

Abstract

Conservation of biodiversity is growing in interest, and wetlands are disappearing at an alarming rate, so understanding how communities are assembled and how interactions among species and ecosystems influence evolution is critical to the management of threatened habitats. We compared diversity and assemblages of peatland Diptera within and between ecoregions in Québec, Canada. We then determined the phylogenetic structure of peatland Diptera communities and how the structure differed with spatial scale (trap, site, ecoregion). Finally, we tested alpha and beta diversity along environmental and spatial gradients to determine which processes influence Diptera communities and diversity. Bogs across the three ecoregions support similar abundance, species richness, and functional diversity. We found that the major forces structuring Diptera assemblages in bogs across Québec are stochastic processes such as dispersal limitations. However, those random patterns change to clustering when anthropogenic disturbances modify the landscape. Assembly rules are mostly dictated by patch and landscape parameters specific to each ecoregion affecting dispersal and establishment between sites. Conservation of mobile organisms in habitats such as bogs will depend on conservation plans focusing on both patch quality and surrounding landscape, and that different conservation strategies need to be applied in different ecoregions.

Introduction

Two evolving and complementary fields have made progress in explaining how local communities are assembled from the regional species pool: metacommunity ecology and community phylogenetics. Metacommunity ecology aims to predict species distribution, abundance, and interactions on a variety of spatial and temporal scales (Leibold et al. 2004; Gravel et al. 2006). Alternatively, community phylogenetics measures how species traits and relatedness explain patterns of community structure in an evolutionary framework (Webb 2000; Cavender-Bares et al. 2009; Peres-Neto et al. 2012). Two theoretical approaches have primarily been used to explain the distribution of species in communities: niche (Diamond 1975; Webb 2000) and neutral (Bell 2001; Hubbell 2001) models. The evolution of species traits, environmental filtering, and competitive exclusions are processes that can structure a community in a phylogenetic context, with constituent species forming a phylogenetically clustered or overdispersed assemblage, supporting the niche model. In contrast, a random pattern of phylogenetic relationships in a community can indicate that these processes are not dominant in structuring the community and that stochastic processes such as chance colonisation, random extinction, and ecological drift (random changes in the relative abundance of species) are generating the patterns observed, thus supporting the neutral model (Hubbell 2001; Hardy and Senterre 2007). The distribution and abundance of organisms are influenced by different forces at different spatial and temporal scales (Graham and Fine 2008; Kraft and Ackerly 2010). Taxonomic, functional, and phylogenetic diversity are complementary metrics that can be used to quantify different components of biodiversity (Petchey and Gaston 2002; Magurran 2004; Cadotte et al. 2010). Consequently, using a multimetric approach is essential to understand the basic principles of community ecology, diversity, and composition to predict new combinations of species and environments, functioning of ecosystems, and the benefits and problems of these changes to society (Jackson and Blois 2015).

Peatlands are an important habitat that fulfils important ecosystem functions, such as improving water quality, protecting shorelines, stabilising water supplies, moderating the effects of floods, and recharging groundwater aquifers (Keiper et al. 2002). They are ecosystems saturated in acidic water, nutrient-poor, low in calcium and magnesium, and covered with peat, which create anaerobic conditions and characteristic vegetation composition (Gorham and Janssens 1992; National Wetlands Working Group 1997). Yet they are threatened systems due to anthropogenic activities, because in Canada more than 60% of wetlands have been transformed through agriculture, industrial/urban development and contamination, climate change, and non-native invasive species (Environment Canada 1993). Given the strong environmental constraints and specificity of the habitat, flora and fauna assemblages have been found to be characteristic (Batzer and Boix 2016). Insect assemblages, flies in particular, although poorly known, are diverse in peatlands and are influenced by both small-scale resource specificity (such as soil acidity) and large-scale factors (such as anthropogenic disturbance of the surrounding matrix) (Beaulieu and Wheeler 2001; Savage et al. 2011; Grégoire Taillefer and Wheeler 2017). Because of this, peatlands and their related fly assemblages represent a good study system where an assemblage is physically constrained with strong environmental pressures.

Flies (Diptera) are an ideal taxon for studying taxonomic and phylogenetic community composition in peatlands at broad spatial scales. Diptera are abundant, species-rich, and ecologically and trophically diverse in peatlands (Blades and Marshall 1994; Pollet 2001; Keiper et al. 2002; Grégoire Taillefer and Wheeler 2012). Diptera are associated with a variety of terrestrial, semi-aquatic, or aquatic peatland habitats including mud shores, vascular plants, decaying organic matter, emergent vegetation, macrophytes, and algal mats (Ferrar 1987; Keiper et al. 2002). Some Diptera species are peatland specialists, because their larval or adult food resource or breeding media is either restricted to or characteristic of peatlands (Teskey and Burger 1976; Farkas and Brust 1986; Marshall 1994).

In this study, we use a multimetric approach (taxonomic, phylogenetic, and functional) to assess the composition and diversity of Diptera communities in peatlands across multiple ecoregions in Québec, Canada. Our model group was 21 families within Schizophora that encompasses a large range of life history and biological traits, and which are abundant and species-rich in peatlands. Consequently, we examined the taxonomic, functional, and phylogenetic species composition and diversity patterns of peatland Diptera in three Québec ecoregions, arranged broadly from south to north, and within and between ecoregions.

We tested four hypotheses on Schizophora (Diptera) diversity patterns between the three ecoregions: (1) Species richness will increase from north to south and differences in traits and species composition will be mostly related to climate. (2) Smaller and more disturbed wetlands will be characterised by lower overall trait diversity and dominated by species possessing traits related to higher resilience (Townsend et al. 1997; Vieira et al. 2004): smaller body size, shorter life cycles, and good dispersal ability. (3) More anthropogenically transformed landscape matrices will support increasing phylogenetic clustering. If traits susceptible to fragmentation are phylogenetically conserved, entire clades could be lost and remaining individuals would tend to co-occur with close relatives, reducing the phylogenetic distance between them (Webb et al. 2002; Cavender-Bares et al. 2009). (4) Species and phylogenetic beta diversity at the landscape scale will be related to local characteristics such as soil cover composition and surrounding characteristics such as community isolation, wetland size, and/or percentage of surrounding forest and wetlands (Patitucci et al. 2011; Savage et al. 2011). At the regional scale, climate will explain most of the variation in species and phylogenetic beta diversity (Økland et al. 2005).

Materials and methods

Study sites, sampling, and processing

The three selected ecoregions (Ricketts and Imhoff 2003) differed by both climatic conditions and human occupancy, which influence the distribution of vegetation and animals (Fig. 1): Eastern Canadian Forest (Montréal Metropolitan area and surroundings, 45.2°N, 73.9°W; region code: Montréal, MTL); Eastern Great Lakes Lowland Forest (Saguenay area, 48.8°N, 72.2°W; region code: Saguenay, SAG); and Central Canadian Shield (southern James Bay area, 49.83°N, 77.27°W; region code: James Bay, BJM). To standardise the peatland type across ecoregions, we selected five open-raised bogs per ecoregion with pH < 4.5 dominated by Sphagnum Linnaeus (Sphagnaceae) mosses, Ericaceae shrubs, a tree cover < 30%, and no open pools of water. Sites in the Eastern Canadian Forest ecoregion were chosen from high-resolution digital aerial photographs available from Ducks Unlimited Canada (Beaulieu et al. 2010) and through communication with experts. Sites in the Eastern Great Lakes Lowland Forest ecoregion were chosen from peatlands inventoried by Calmé et al. (2002). Sites in the Central Canadian Shield ecoregion were chosen in the field by driving along the James Bay Road (Route de la Baie-James).

Fig. 1. Map of the 15 study peatland sites across Québec, Canada. Map is separated into ecoregions: Montréal – Eastern Great Lakes Lowland Forest; Saguenay – Eastern Canadian Forest; James Bay – Central Canadian Shield (Ricketts and Imhoff 2003). Map created with SimpleMappr (www.simplemappr.net).

Sampling was conducted weekly for five weeks during the period of highest Diptera species richness and activity. Montréal sites were sampled from 28 June to 31 July 2013, Saguenay sites from 2 July to 8 August 2014, and James Bay sites from 29 June to 29 July 2015. In each site, samples were collected using yellow pan traps (4 cm deep and 12 cm of diameter) and sweeping in an area of 30 × 30 m, at least 30 m from the edge. Sweep samples were collected weekly on three random transects of 20 m with 20 sweeps on each transect in suitable weather conditions. Three transects of four pan traps each were placed 10 m apart in the 30 × 30 m plot. Pan traps were placed in the soil with their upper rim flush with the ground surface and filled with a 50% solution of propylene glycol and water, with a drop of liquid detergent as a wetting agent. Pan traps were serviced every 6–8 days.

Insects were preserved in 95% ethanol for subsequent DNA extraction. Small flies were dried using hexamethyldisilazane and then mounted; larger flies were transferred into ethyl acetate, then pinned, and air-dried. Specimens were deposited in the Lyman Entomological Museum (McGill University, Ste-Anne-de-Bellevue, Québec). Taxa were identified by morphological characters to named species or morphospecies. Target Schizophora families were selected based on their abundance, diversity of functional traits, and availability of identification keys and expert taxonomists. Phoridae, Tachinidae, Muscidae, and Anthomyiidae were excluded as these families are complex to identify using morphology, or experts were not available for identification. Finally, each species in the 21 families selected was assigned to a category in each of eight biological traits: feeding habits, size, specialisation, habitat preference for oviposition, preferred substrate, voltinism, overwintering stage, temperature range, and wetland specialisation. Also, three indirect rarity measures were used: frequency of occurrence, range, rarity (Table 1). Frequency was the number of sites in which the species occurred, and rarity was the number of individuals in the total collection across the 15 bog sites. Trait values were determined from direct measurements/counts, published literature, and consultation with Diptera specialists.

Table 1. Life history and ecological traits considered in this study including categories, abbreviations, and trait determination notes.

Traits were assigned at the species level and adult stage except for feeding habits and specialisation (larval stage).

Habitat and vegetation variables

At the local scale, the following ground cover attributes were surveyed at week three in five 1 × 1 m quadrats placed at each corner and in the middle of the 30 × 30 m area at each site: percent plant species cover, litter cover, canopy cover, open water cover, and bare soil cover. Fluctuation of water table depth (cm) was estimated with the polyvinyl chloride tape discoloration method (Belyea 1999). Polyvinyl chloride tape was mounted along 1-m-long bamboo sticks, one of which was inserted vertically in each site with 15 cm left above the surface and left for the five weeks of sampling. Chemical parameters of the substrate (pH, conductivity, and temperature) were measured with Hanna pocket combo pH / conductivity / TDS tester in the field at week three at two random locations within the 30 × 30 m area.

At the regional scale, wetland size and the nature of the surrounding matrix were quantified using QuantumGIS version 1.8.0 (Quantum GIS Development Team 2012) software with the Earth Observation for Sustainable Development of Forests land cover datasets of year 2000 representing 23 land cover classes at a spatial resolution of 25 m (Wulder et al. 2008). Land use (open water, exposed land, urban development, low vegetation, wetland, forest, and agriculture) surrounding each wetland was measured (m²) within a circle of 2 km radius. Diptera have been found to respond to the surrounding matrix at this scale (Meats and Smallridge 2007; Savage et al. 2011). Climatic variables that are known to influence Arthropoda diversity (Bowden and Buddle 2010) were also extracted from WorldClim version 1.4 (Hijmans et al. 2005) in 30-arcsecond resolution. Values at the site locations were extracted with raster, rgdal, and foreach packages in R version 3.2.2 (R Core Team 2015). The seven variables selected were annual mean temperature, maximum temperature of warmest month, minimum temperature of coldest month, mean temperature of warmest quarter, mean temperature of coldest quarter, temperature seasonality (standard deviation * 100), and annual precipitation.

Statistical analysis: taxonomic and functional composition

Diptera communities in each peatland (Supplementary Table 1) were characterised by different biodiversity metrics to account for the probable difference in sampling effort among sites; two common approaches for comparing diversity were used: rarefaction of the data to the lowest effort (Srare) and use of a non-parametric estimator (Chao1) that is robust to unequal sampling efforts and that can quantify sampling completeness. Diversity was assessed with rarefaction curves based on 1000 permutations. Extrapolated species richness was assessed using a bias-corrected Chao index (O’Hara 2005). Observed abundance (Abun) and observed species richness (S) were also calculated to assess how common or rare a species is relative to other species in a community. All measures were calculated from samples pooled by sites throughout the season. Significance of differences in abundance and species richness indices among ecoregion was determined based on analysis of variance F tests. All analyses were performed using the R vegan package (Oksanen et al. 2012).

Functional diversity, which considers the distribution and range of functions of co-occurring species (Supplementary Table 2) in a community, was measured with three multidimensional indices (Villéger et al. 2008). Functional richness (FRic) was standardised by the total FRic to constrain the values between 0 and 1. Functional evenness (FEve) and functional divergence (FDiv) were weighted by the abundance of species. First, a Gower dissimilarity matrix was computed (Podani 1999) because traits were quantitative and categorical, and lingoes corrections were applied to obtain Euclidean distance matrices. Dimensionality of the trait matrix was reduced to 40 out of 157 principal coordinate analysis (PCoA) axes due to computational power limitations. Significance of differences between sites and ecoregion was determined with analysis of variance F tests. Community-level weighted trait means (Lavorel et al. 2008) were used to determine dominant functional composition at each site, where quantitative traits are weighted by abundance and categorical traits are returned as the dominant category. All analyses above were computed with the FD package (Laliberté et al. 2014).

Redundancy analysis was used to test for the similarity in overall taxonomic and functional assemblages and for the relationships to environmental variables (Supplementary Table 3) using the vegan package. Prior to redundancy analysis, species abundances were Hellinger-transformed (Legendre and Gallagher 2001) as it gives low weights to low counts and many zeros. A PCoA of the dissimilarity matrix of community-level weighted trait means was computed; principal coordinate eigenvalues of community-level weighted trait means were used as response variables in the subsequent redundancy analysis. Principal component analysis was used to reduce the number of ground cover variables (Supplementary Table 3). Significant axes were identified with the Kaiser-Guttman criterion (Yeomans and Golder 1982), and the first five significant principle component analysis axes retained 82% of the variation explained (Supplementary Table 4). To reduce the number of environmental variables in the subsequent analyses and to avoid autocorrelation between them, remaining environmental and climatic variables were observed with scatterplots. Only soil temperature was removed. Latitudes and longitudes were transformed into corresponding coordinates in X (east–west) and Y (north–south) distances in the package SoDA (Chambers 2008). The resultant environmental matrix included local variables (PC1–5, canopy cover, open water cover, fluctuation of water table, pH, conductivity) and all landscape variables (open water, exposed land, urban development, low vegetation, wetland, forest, agriculture, area). Forward selection of explanatory variables was applied with the packfor package and retained the variables with a P-value < 0.05. Spatial (X and Y) and climate variables were excluded from redundancy analyses and incorporated further in variation partitioning. A permutation test was used for the significance of the axes eigenvalues associated with significant environmental variables. A permutational multivariate analysis of variance based on a Bray–Curtis distance matrix was used to assess the significance of differences among ecoregions for overall species assemblages and environmental variables using the function Adonis of the vegan package.

A fourth-corner analysis (Dray and Legendre 2008) using the ade4 package (Dray and Dufour 2007) was used to assess the relationship simultaneously between separate traits, species abundance, and environmental/climatic variables. Missing entries (251) of the categorical traits were replaced by predicted values with the MissMDA package (Josse and Husson 2016), which corresponded to 11% of all entries. Associations between categorical traits and quantitative environmental variables were measured with Pearson correlation coefficient. The significance was tested by a combination of the permutation model 2 and 4 with 999 permutations to obtain a correct level of Type I error.

Phylogeny reconstruction

Hypothesised molecular phylogeny of Diptera species from peatland sites was generated with DNA barcode (658 base pairs of the mitochondrial CO1-5P gene) used as a proxy for species-level phylogenetic relationships (Joly et al. 2014). Boyle and Adamowicz (2015) found that trees built using COI alone calculate the metrics of phylogenetic community structure as accurately as the trees obtained from a multi-gene tree. A compilation of published molecular sequences in the Barcode of Life Data Systems (Ratnasingham and Hebert 2007) and our own Diptera species representing 356 specimens and 139 species submitted to the Canadian Centre for DNA Barcoding for sequencing using C_LepFolF/C_LepFolR primers (Hernández-Triana et al. 2014) were used in a matrix of DNA sequences constructed with Mesquite 3.04 (Maddison and Maddison 2015) and aligned with MUSCLE (Edgar 2004). A Bayesian approach was used for the phylogenetic analyses using BEAST version 1.7.5 (Drummond et al. 2012), and the output was examined with Tracer version 1.6 (Rambaut et al. 2014). To determine which model of nucleotide substitution best fit the alignment, jModelTest2 according to Akaike information criterion was applied. We constrained monophyly for families, subfamilies, genera, and groups for which there were recent and well-supported multi-gene or morphological phylogenies. These included Schizophora and some higher groups of acalyptrate Diptera (Yeates and Wiegmann 2005); Calyptratae (Lambkin et al. 2013); Sarcophagidae (Kutty et al. 2010; Pape et al. 2011); Ravinia Robineau-Desvoidy (Sarcophagidae) (Giroux et al. 2010; Piwczynski et al. 2014); Boettcheria Parker (Sarcophagidae) (Piwczynski et al. 2014); Luciliinae and Polleniinae (Calliphoridae) (Kutty et al. 2010); Phytomyzinae (Agromyzidae) (Scheffer et al. 2007); Scathophagidae (Bernasconi et al. 2000); Milichiidae (Brake 2000); Sciomyzini and Tetanocerini (Sciomyzidae) (Tóthová et al. 2013); Tephritidae (Han and Ro 2009); Drosophilidae (Yassin 2013); and Chloropinae (Chloropidae) (Brake 2000). Sequences of Dolichopus brevipennis Meigen (Dolichopodidae), Sphaerophoria philanthus Meigen (Syrphidae), Chelipoda truncata MacDonald (Empididae), Hypocera ehrmanni Aldrich (Phoridae), and Lonchoptera furcata (Fallén) (Lonchopteridae) were used to root the tree (outgroup).

The general time-reversible model with a proportion of invariable sites (+I) and a rate of variation across sites (+G) was selected as the appropriate model with jModelTest2, and no partitioning by codon was used. In all analyses, strong priors were set on the age of two nodes corresponding to Schizophora (lognormal: mean = 3, SD = 0.78, offset = 70) and Chloropidae (lognormal: mean = 3, SD = 0.7, offset = 42) according to paleontological and molecular data (Nardi et al. 2010; Wiegmann et al. 2011). A branching prior was set under a Yule process model, and a relaxed molecular clock was assumed using a lognormal distribution of rates (Drummond et al. 2006). The analyses were performed twice using a random starting tree for a Markov chain Monte Carlo chain length of 40 million generations (10 000 echo states, 4000 log parameters, 20% burn-in) (Fig. 2).

Fig. 2. Hypothesised phylogenetic relationships, based on CO1, among Schizophora (Diptera) species collected in the 15 bogs. Colours represent families; branch lengths represent divergence time estimates; and numbers on nodes represent posterior probabilities. Species codes can be found in Supplementary Table S1.

Phylogenetic community structure

The phylogenetic structure was first calculated for all Diptera pooled in each trap and repeated at each higher spatial scale: for each site and in each ecoregion. To determine if Diptera communities are more or less related than expected by chance, several indices were calculated: Faith’s phylogenetic diversity metric (PD) (Faith 1992), mean pairwise distance (MPD), and the mean nearest taxon distance (MNTD) were calculated for each community (Webb 2000; Webb et al. 2008). A null model shuffling taxa label to generate null communities randomised 999 times was used to determine if the phylogenetic composition of communities differed significantly from that expected by chance. These 999 null measurements were used to calculate two measures of standardised effect size multiplied by −1: net relatedness index (NRI) and nearest taxon index (NTI) (observed value minus the mean value of the 999 null values divided by the standard deviation of the 999 null values) and probabilities. To estimate whether net relatedness index and nearest taxon index values at the trap and site scales were significantly different from zero, t-tests were performed. Fisher’s combined probability tests were performed as meta-analyses of P-values using the metap package (Dewey 2016). Null models were used at the ecoregion scale to assess if communities were significantly different from random.

Phylogenetic signals in trait diversity

The phylogenetic signal is defined as the tendency for related species to resemble each other. Trait conservatism is related to a higher degree of phylogenetic signal, meaning that close relatives share similar traits, while trait convergence is the tendency for distantly related species to resemble each other more than expected (Blomberg et al. 2003). Phylogenetic signals of continuous traits (size and temperature range) were tested with Blomberg’s K, which is rescaled by dividing by the Brownian motion expectation. This gives it the property of having an expected value of 1.0 under Brownian evolution. For significance we randomly arrayed the trait data on the community phylogeny 999 times to generate a null distribution from which a P-value could be calculated using the picante package (Kembel et al. 2010). Phylogenetic signals of categorical traits were measured using Pagel’s lambda with the fitDiscrete function in the Geiger package (Harmon et al. 2008). By comparing the likelihood ratio test of a model where the tree is transformed by lambda with one where the tree is transformed into a large polytomy (λ = 0), we could predict which one fits better by a chi-square distribution.

To evaluate the variation of alpha diversity along anthropogenically modified landscapes, a series of linear models was used to evaluate linear relationships among measures of Abun, S, Srare, Chao1, FRic, FDiv, FEve, NRI, and NTI to the proportion of anthropogenically modified land (urban development plus agriculture). Regression diagnostic plots were created in R to visually assess if linear regression assumptions are met. Dependant variables were log (x + 1)-transformed to meet the assumptions of the test.

To evaluate the phylogenetic turnover in Diptera communities within and between the three ecoregions, beta phylogenetic differences were tested using comdist (betaMPD: the average mean pairwise distance for each species in a sample to all species in another sample) and comdistnt functions (betaMNTD: the average mean nearest taxon distance for all species in a sample to the nearest neighbours in another sample) of the picante package. Then PCoA of the output dissimilarity matrix was computed, and principal coordinate eigenvalues were used as response variables in redundancy analysis to determine the effect of environment on beta phylogenetic structure. To determine if beta and phylogenetic beta diversity were positively spatially autocorrelated, the Mantel statistic was calculated with spatial distance. To assess correlations between species co-occurrence and phylogenetic distances, the Mantel statistic was also performed between Bray–Curtis dissimilarity matrices and beta phylogenetic dissimilarity matrices.

Variation partitioning was used to assess the proportion of variation in community data using adjusted R ² in redundancy analysis ordination explained by environmental, spatial, and climatic variables at the taxonomic, functional, and phylogenetic levels. The variation was partitioned into several components depending on significance according to forward selection: pure environmental component (local and surrounding variables), pure climatic component, pure spatial component, spatially structured environmental and/or climatic components, and residual variation. Analysis of variance was used to calculate the significance levels of different components. These analyses were computed with the vegan package.

Results

Species and functional alpha and beta diversity

A total of 3302 specimens belonging to 202 species in 21 families were identified from the 15 bog sites in the three ecoregions (Supplementary Table 1). Neither abundance, observed species richness, rarefied species richness nor Chao1 indices were significantly different among and within the three ecoregions (P > 0.05), although there was a tendency for a decrease in diversity from south to north. FRic, FEve, and FDiv were not significantly different among sites and between the three ecoregions (P > 0.05). FRic, the amount of niche space occupied by the species within a community, was low across all sites; FEve was moderately high and FDiv was high (Table 2).

Table 2. Abundance (Abun), observed species richness (S), Chao1 indices (Chao1), rarefied species richness (Srare), functional richness (FRic), functional evenness (FEve), and functional divergence (FDiv) for all sites and mean diversity for each ecoregion.

Site codes: MTL, Montréal (Eastern Canadian Forest); SAG, Saguenay (Eastern Great Lakes Lowland Forest), and BJM, James Bay (Central Canadian Shield).

The community-weighted means showed that species inhabiting Montréal and Saguenay sites possessed similar traits, predominantly saprophagous generalist species preferring moist substrates for oviposition (Table 3). The functional composition of James Bay sites was clearly separated from the other two ecoregions as communities were dominated by saprophagous and predacious species specialised on insects or arthropods that use terrestrial conditions for oviposition. Nearly all communities were dominated by widespread wetland facultative species found in at least six sites and in the three ecoregions. Several significant relationships among species functional traits and environment and climate were found with fourth corner analysis (Fig. 3). The saprophagous generalist species inhabiting Montréal and Saguenay bogs were positively influenced by the coverage of wetlands surrounding the sites, although the warm temperature negatively impacted them. Insect predators in James Bay bogs were positively influenced by the fluctuation in the level of water, longitude, forests surrounding the sites, and tree cover within the sites. Habitat preference for oviposition, voltinism, size, frequency of occurrence, and rarity were not significantly correlated with any measured environmental parameters.

Table 3. Dominant functional composition of each site, represented by community-level weighted means.

For continuous traits, the mean trait value is weighted by abundance, while for categorical traits, the dominant class is returned. For description of traits and associated codes, refer to Table 1.

Fig. 3. Relationship between the functional traits and environmental/climatic variables of Diptera species revealed by fourth corner analysis. Blue cells correspond to significant negative relationships; red cells correspond to significant positive relationships. Refer to Table 1 for an explanation of each trait and code.

The three ecoregions were significantly different (Adonis R² = 0.47, P = 0.001) in species assemblages. Climate was significantly different between the three ecoregions (Adonis R² = 0.98, P = 0.001). Local and landscape variables were also significantly different between ecoregions (Adonis R² = 0.77, P = 0.001). The total proportion of variance explained in species redundancy analysis by four explanatory variables (surrounding forests, pH, PC1, PC2) was 53% (adjusted R ² = 0.34) (Fig. 4A) for all sites combined. Montréal species composition in redundancy analysis was mostly affected by higher pH, PC2, and PC1, while Saguenay species were not affected by precise environmental variables. James Bay species composition was primarily affected by surrounding forests and by higher coverage of bog specialist plant species Maianthemum trifolium (Linnaeus) (Asparagaceae) and Chamaedaphne calyculata Baillon (Ericaceae) (PC1). Montréal and Saguenay ecoregions were not significantly different in functional composition (Adonis R ² = 0.98, P = 0.18). Four different environmental variables (surrounding forests, area, PC2, PC4) emerged in the redundancy analysis as having a significant effect on the functional composition for all sites combined explaining 51% (adjusted R ² = 0.39) of variance (Fig. 4B). Large area predominantly affected the functional composition of Montréal (sites 3–4–5) and Saguenay bog sites, in being composed mainly of saprophagous generalist species preferring moist substrates for oviposition. Two Montréal sites that were dominated by specialist phytophagous species (community-level weighted trait means) were mostly affected by a high coverage of Carex oligosperma Michaux (Cyperaceae) and low litter cover (PC2), and high coverage of Polytrichum Hedwig (Polytrichaceae) species and low Eriophorum vaginatum Linnaeus (Cyperaceae) (PC4). James Bay communities were mostly affected by a high percentage of surrounding forest.

Fig. 4. Redundancy analyses of species (A) and functional composition (B). Arrows represent vectors of significant variables (P < 005) explaining community structure, R ² values in parentheses from the forward selected models examining the effect of environmental factors.

Phylogenetic alpha and beta composition

Montréal communities at the three scales were significantly different from random for net relatedness index and nearest taxon index, and Faith’s phylogenetic diversity metric was also significantly different at the site scale (Table 4). At the trap scale for Saguenay, the three indices were significantly different from random, while at the site scale Faith’s phylogenetic diversity metric and nearest taxon index were significant and at the ecoregion scale only Faith’s phylogenetic diversity metric was different from random. James Bay communities were not different from random at every scale. Mantel correlation between beta and phylogenetic beta diversity was significantly correlated with spatial distance for species (r = 0.68, P = 0.001), mean nearest taxon distance (r = 0.53, P = 0.002), and traits (r = 0.43, P = 0.003), although not mean pairwise distance (r = 0.001, P = 0.45). Species co-occurrence was significantly correlated (P < 0.05) with phylogenetic distances and functional distances.

Table 4. Phylogenetic diversity (PD), net relatedness index (NRI), and nearest taxon index (NTI) for Diptera communities in 15 sites across three ecoregions at three spatial scales (trap, site, ecoregion).

Fisher’s combined probability tests were performed to assess if overall P-values were significant at the trap and site scales. Null models were used at the ecoregion scale to assess if communities were significantly different from random. Significant P-values are in bold font.

The three ecoregions were significantly different in phylogenetic beta diversity (betaMNTD Adonis R ² = 0.68, P = 0.001; betaMPD Adonis R ² = 0.38, P = 0.001). The total variance explained by the redundancy analysis in betaMNTD for all sites combined by four explanatory variables was 31% (adjusted R ² = 0.22) (Fig. 5A). Phylogenetic structure at the Montréal sites at the tip of the tree was mostly influenced by substrate coverage (PC1), while Saguenay communities were mainly affected by the large bog areas and substrate coverage (PC2) in this ecoregion. James Bay communities were affected by low vegetation surrounding the sites. The total variance explained by the redundancy analysis in betaMPD for all sites combined by two explanatory variables was 19% (adjusted R ² = 0.06) (Fig. 5b). Surrounding exposed lands were affecting the majority of Montréal communities, while surrounding forests were affecting James Bay communities. Phylogenetic turnover was higher deeper in the phylogeny (betaMPD) for Montréal communities than within the other two ecoregions (Fig. 5B), as sites were more spread apart. Forward selection analyses showed that overall betaTD, betaFD, and betaMNTD were significantly associated with geographical distances, while betaMPD was only significantly associated with east–west coordinates. Variation partitioning analyses (Fig. 6) explained 45% of the total variation of betaTD with the shared spatial, environmental, and climatic fraction explaining the highest proportion of variation (24%, P = 0.001). Pure environmental variables accounted for 8% (P = 0.001) of the variation, while climate accounted for 2% (mean temperature of warmest quarter (TempWarm), maximum temperature of warmest month (MaxTemp, P = 0.001). The shared component between environmental, spatial, and climatic variables accounted for 18% (P = 0.003), the pure environmental fraction accounted for 14% (P = 0.002), and climate accounted for 9% (TempWarm, MaxTemp; P = 0.001) of the total variation (50%) of betaFD. The shared fraction of spatial, environmental, and climatic (annual mean temperature (AnnuTemp, MaxTemp) variables explained 16% (P = 0.001) of the total variation (26%) of betaMNTD. The shared spatial, environmental, and climatic fraction explained the highest proportion of variation (4%, P = 0.001) of the total variation (11%, P = 0.01) in betaMPD. The pure climatic variables (TempWarm, MaxTemp) explained only 0.7% (P = 0.001) of the variation.

Fig. 5. Redundancy analyses of beta phylogenetic composition for mean nearest taxon distance (A) and mean pair-wise distance (B). Arrows represent vectors of significant variables (P < 005) explaining community structure, R ² values in parentheses from the forward selected models examining the effect of environmental factors.

Fig. 6. Variation partitioning diagrams of taxonomic (betaTD), functional (betaFD), and phylogenetic (betaMNTD, betaMPD) composition. Circles represent variation explained by unique and shared fractions of significant (based on forward selection) environmental (Env), climatic, and spatial variables, while numbers correspond to the adjusted coefficient of determination (adjusted R ²) associated with each circle. Only significant and positive adjusted R ² are shown.

Linear models indicated that species richness was significantly positively related to agriculture (df = 12, F = 3.29, adjusted R ² = 0.25, P = 0.04). Urban development and agriculture had an additive effect on FDiv with a negative and positive association respectively (df = 12, F = 8.73, adjusted R ² = 0.52, P = 0.005). Phylogenetic diversity had a similar association with agriculture as species richness (df = 12, F = 4.04, adjusted R ² = 0.30, P = 0.03). Nearest taxon index was positively related to urban development (df = 12, F = 5.02, adjusted R ² = 036, P = 0.012); therefore, with increasing urban coverage, clustering is increasing.

Tests of phylogenetic signal showed that traits showing a signal are phylogenetically labile, with close relatives less similar than expected under a Brownian model of evolution. The K statistic was significantly different than a signal expected by chance but less similar (K < 1) than expected under a Brownian model for size (K = 0.68, P = 0.001) and temperature range (K = 0.29, P = 0.01). The variance in size and temperature range found within clades could be caused by either departure from Brownian motion evolution, such as adaptive evolution that is uncorrelated with the phylogeny, or measurement error. Pagel’s lambda was not significantly different (P > 0.05) from a giant polytomy for all other categorical traits.

Discussion

Only a few studies have examined the diversity of wetland Diptera (e.g., Blades and Marshall 1994; Barták and Roháček 1999; Grégoire Taillefer and Wheeler 2010, 2012, 2013, 2017), and even fewer have looked at their phylogenetic community structure (Pfenninger et al. 2007; Espíndola et al. 2012; Silver et al. 2012). Phylogenetic relatedness and functional structure complemented the information provided by classical metrics. Different dimensions of Diptera diversity in bogs across different ecoregions are more evident with this multimetric approach and help to predict if communities are mainly maintained by competition-driven differentiation, environmentally filtered coexistence, or stochasticity. Bogs across the three ecoregions support similar abundance, species richness, and functional diversity. Our results suggest that major forces structuring Diptera assemblages in bogs across Québec are stochastic processes such as dispersal limitations. However, those random patterns change to clustering when anthropogenic disturbances modify the landscape. Assembly rules are mostly dictated by patch and landscape parameters specific to each ecoregion affecting dispersal and establishment between sites.

We did not find evidence of limiting similarity as a mechanism of Diptera community assembly in bogs at any spatial scale. A random pattern can result when opposing patterns of clustering and overdispersion are found within the same community (Hardy and Senterre 2007; Mayfield and Levine 2010; Pangjanda and Pramual 2016). However, biotic interactions have rarely been considered important mechanisms to control invertebrate species biodiversity in wet environments (Batzer and Wissinger 1996; Vinson and Hawkins 1998), and this still holds with the advent of phylogenetic analyses (Ruhí et al. 2013; Grégoire Taillefer and Wheeler 2017). Our results suggest that dispersal limitation of small-sized species, coupled with environmental filtering caused by the strength of the environmental gradients, is driving the spatial patterns observed in these bog Diptera communities. The most common species collected had small sizes, and many of these species may be dispersal-limited because several acalyptrates do not appear to fly readily or for long distances (T.A.W., personal observations). Therefore, they likely benefited from dispersal by wind (Kovats et al. 1996; Heino 2013) over a long period of time, which explains their occurrence in the three ecoregions. A historical imprint caused by dispersal limitations over time likely created these phylogenetic random patterns (Hubbell 2001). In two of our other studies on wetland Diptera communities in Québec, we found that inter-site distance mostly affected beta taxonomic and phylogenetic diversity of small and common species (Grégoire Taillefer and Wheeler 2017), and that bog disturbance had the greatest effect on the distribution of small Brachycera and Acalyptrate flies (Grégoire Taillefer and Wheeler 2012). Abundant, small-bodied, multivoltine insect species were affected by stochastic processes in tropical stream metacommunities (Saito et al. 2015a). These functional traits are similar to our dominant traits in every bog communities. As local variations in species diversity and composition are dependent on historical diversification and dispersion, a larger latitudinal gradient would have been necessary to detect a change in species diversity and an effect of species pools (Martin 2016).

However, with recent drastic human modification of the landscape, less suitable patches and more barriers to dispersal are found neighbouring those bogs, so agriculture and urban development act as filters for the small proportion of species in the species pool that can disperse in these conditions. Contrary to our hypothesis, pure climatic variables were not the major determinants of community structure at the regional scale. In Montréal, assembly rules are mostly dictated by patch and landscape parameters affecting dispersal and establishment between sites in an urban matrix. The high turnover of species in Montréal sites is mostly affected by surrounding exposed lands (e.g., non-vegetated, non-developed, sediments, beaches, burnt area) where Diptera species probably have environmental requirements suitable for bog habitats and local conditions. As species emigrate from the surrounding habitats, they are filtered by anthropogenic barriers and the unique local vegetation. A selective filtering role of anthropogenic disturbances has been found for several other taxa (Brunbjerg et al. 2012; Ding et al. 2012; Concepción et al. 2016; Hausberger and Korb 2016; Mykrä et al. 2016). Increase environmental stresses in aquatic Diptera wetland communities have also been found to increase phylogenetic clustering in Alberta, Canada (Silver et al. 2012).

Bog communities are distributed along broad environmental gradients in surrounding land use and anthropogenic disturbance. Montréal bogs can be seen as more isolated systems, where exchange with the anthropogenically disturbed surrounding environment is difficult. On average 28% of the proportion of the 2-km-radius area around each bog was anthropogenically modified. Isolation is also illustrated by the higher variation in betaMPD observed compared with the other two ecoregions, which means that species composing Montréal sites were found in different families. Our results are similar to Angold et al. (2006) who found that anthropogenic activities and urban environmental conditions influence ground beetle (Coleoptera: Carabidae) species diversity and composition. Vegetation structure, level of disturbance affecting the dominance of generalists, and woodland ground beetle specialists are the most significant influence on assemblages in urban habitat patches. Saguenay bogs are more intermediate in surrounding land use to Montréal and James Bay. They have less anthropogenic disturbances than Montréal and more wetland coverage than James Bay (approximately 5% of anthropogenic land use, 36% forest, 43% wetlands). James Bay bogs are more open systems, where exchange can happen with the surrounding undisturbed environment (approximately 0.03% of anthropogenic modification). Predacious insect species in James Bay immigrate from the surrounding forests as the bogs in this region are connected to a complex of forests (approximately 80% land cover) and wetlands (approximately 9% land cover). Consequently, community assembly within regions with high fragmentation of habitat, characteristic local conditions, and low dispersal rates usually converge in phylogenetic clustering of assemblages; whereas in regions where dispersal is prevalent, assembly patterns should be largely stochastic and dominated by species dispersal capacity (Emerson and Gillespie 2008; Weiher et al. 2011).

High levels of functional divergences are associated with a high degree of niche differentiation among species within communities: the most abundant species are very dissimilar and likely compete weakly. As shown in the community-level weighted trait means analysis, common species were facultative wetland inhabitants, and FDiv was increasing with the coverage of agricultural land but decreasing with urban development. Consequently, the recent species turnover in Montréal must originate from species inhabiting surrounding agricultural land that possess very dissimilar traits to peatland inhabiting species. And as these peatlands possess unique features in an urban matrix, clustering at the tip of the phylogeny (NTI) is probably due to in situ environmental filtering (PC1 and PC2) and low dispersal rates of generalist saprophagous species from surrounding wetlands and exposed lands. Another of our study encompassing more environmental heterogeneity in Montréal region (Grégoire Taillefer and Wheeler 2017) showed that bogs have harsher environmental conditions compared with other wetland types, which also act as selective pressure. Therefore, heterogenisation of communities, FRic, and redundancy levels in bogs are systematically limited despite differences in species assemblages.

The large proportion of unexplained variation in beta diversity for each facet of biodiversity considered suggests that other factors are determining community patterns. These factors are probably spatially structured environmental factors, such as a combination of local and surrounding conditions with dispersal limitations. Moreover, taxonomic turnover was high among ecoregions, which increased phylogenetic turnover; and beta diversity was correlated with spatial distance except betaMPD, which was only affected by longitude. Soininen et al. (2007) found in a meta-analysis that shifts in community composition were happening at the scale of 1000 km in line with large-scale variability in environmental conditions and that passive versus active dispersal did not have an effect on species turnover with increasing spatial scale. Here, betaMPD might be correlated more to Pleistocene glaciations and Holocene post-glacial dispersal, as these had a profound impact on Nearctic biotas in a diversity of habitats (Lafontaine and Wood 1988; Harris and Taylor 2010; Solecki et al. 2016). Long-term dispersal limitation in an east–west direction, from coastal to more continental conditions, seems to shape the original beta diversity, as explained variation in betaMPD was mostly driven by pure spatial factors. Arnan et al. (2015) found similar results to ours for the three diversity measures calculated in ant assemblages. Spatial factors thus dispersal limitation played a major role in shaping ant communities in different climatic and human disturbed environments across Europe.

Conclusions

The major forces shown to structure Diptera assemblages in this study in bogs across Québec are stochastic processes such as dispersal limitations of abundant, small, multivoltine species. The random patterns observed change to environmental filtering with increased anthropogenic disturbances modifying the landscape. The massive Diptera radiation that began in the early Cenozoic and the different historical disturbances are probably at the origin of the functional and phylogenetic structure observed for peatland Diptera. Using only traditional metrics, it would not have been apparent that Montréal bog communities are impacted by land use changes and that these impacts change the species pool capable of inhabiting these isolated unique habitats (Saito et al. 2015b). This suggests that the three levels of diversity studied are complementary and give a different picture of macroecological patterns in those threatened habitats. Our results show that conservation of mobile organisms in bogs will depend on conservation plans focussing on both patch quality and surrounding landscape. Different conservation strategies need to be applied in different ecoregions. Montréal should be prioritised for biodiversity conservation through the re-establishment of habitat connectivity (Tewksbury et al. 2002; Staddon et al. 2010), and unique roles of wetlands should be key factors to maintain high species richness. In Saguenay, the large bog areas should be maintained as this diminishes the surrounding landscape pressure; and with the continuing increase of human-induced changes, this will become more important. Finally, protecting clusters of forest and wetlands in James Bay is a more effective conservation strategy than preserving islands of peatlands, as mobile organisms may see clusters of habitat as continuous potential niches instead of unreachable islands (Hatfield and LeBuhn 2007). Developing a standard approach and efficient techniques for large-scale biodiversity monitoring, such as high-throughput DNA sequencing (Gibson et al. 2015) to assess the efficacy of mitigation measures, would be highly valuable.