Drivers of avian species richness and community structure in urban courtyard gardens

Increasing global urbanisation has steered research towards understanding biodiversity in urban areas. Old city spaces throughout Europe have a proliferation of urban court gardens, which can create a mosaic of habitat pockets in an urban area. This article examines the patterns and drivers of avian species richness and community structure in 20 gardens of the constituent colleges of the University of Oxford. We conducted morning surveys across 7weeks in May and June 2017 and used an information-theoretic approach and model averaging to identify important habitat predictors of species richness. We also studied community structure with Sorensen indices and non-metric multi-dimensional analysis. A total of 43 avian species were observed across all sites. Our sites generally differed in their avian assemblages, with greater species turnover than nestedness between sites. Site area was the strongest predictor of site species richness and surrounding habitat composition was the dominant driver of community structure. Thus, the largest gardens were the most species rich, but species composition among gardens differed based on the habitats in which they were embedded. We support using island biogeography theory to understand the avian species assemblages of urban ecosystems and stress the suitability of our study sites for future urban ecosystem research and generating wildlife awareness.


Introduction
The world's urban population is expected to increase to 6.4 billion by 2050 (Gaston 2010). As natural environments continue to shrink in the face of urban growth, conservationists have increasingly recognised the potential for urban conservation (Dunn et al. 2006;Dearborn and Kark 2010;Gaston 2010). Practical biodiversity conservation in populated areas can be informed by studies of the basic ecology of cities (Gaston 2010;Lepczyk et al. 2017b). Birds are a group of organisms that features prominently in these studies (Murgui and Hedblom 2017).
Their conspicuousness, ease to study and role as a taxonomic indicator (Lepczyk et al. 2017b) make them a model taxon for identifying patterns of urban biodiversity and informing conservation decisions in urban environments.
The species assemblages and environmental characteristics of the small city of Oxford have been the subject of research in recent years (e.g. Dale 2015). However, a cohesive investigation into the patterns and drivers of avian species assemblages of central Oxford has not been conducted. This study focuses on avian species richness, community structure and their respective drivers in a sample of 20 University of Oxford Colleges collected during 7 weeks of the 2017 breeding season (from 30 April to 16 June). We also investigated the gardening regimes enacted in these colleges (referred to interchangeably in this article as 'sites'). We hope our study highlights the potential of these sites for future urban ecology research.
We chose to examine Oxford colleges because: (1) college gardens are examples of courtyard gardens that feature in many urban centres in the UK and Europe, (2) each college's defined boundary allows for species-area relationships to be reliably analysed and (3) each college is maintained by a different college gardener and thus likely subjected to different gardening regimes. This last factor creates an interesting sample of urban greenspace where every college is a unique environment-a strong contrast to studying a set of city parks maintained in a uniform manner by local government. This spread of different management practices allows for correlations and regression relationships to be investigated thoroughly.
We also consider the gardening practices implemented within each site. A gardener's actions determine a site's plant assemblage, and this is a particularly important aspect to note, as the extent and composition of a plant assemblage has significant implications for avian richness (Mills, Dunning, and Bates 1989;Ferná ndez-Juricic 2000;Chamberlain et al. 2007;Grimm et al. 2008;Gaston 2010). This investigation into garden maintenance also sheds light on how human-controlled processes and human activities are key contributors towards urban species richness levels (Grimm et al. 2008;Gaston 2010;Faeth, Bang, and Saari 2011;Hanmer 2017), emphasising the importance of human involvement in urban wildlife conservation.
We hypothesise that patterns of avian species richness and community structure will differ among sites. We hypothesise that this will be due to differences in the characteristics of each site, such as site area and site habitat composition. We suspect two sets of factors: (1) those with positive effects, such as patch size (Chamberlain et al. 2007;Gaston 2010), vegetation cover (Ferná ndez-Juricic 2000; Oneal and Rotenberry 2009;Gaston 2010), supplementary feeding by humans (Recher 1972), the presence of water bodies (Chamberlain, Cannon, and Toms 2004;Chamberlain et al. 2007) and a greater proportion of adjacent habitats (Chamberlain et al. 2007;Zhou, Fung, and Chu 2012); and (2) factors with negative effects such as urbanisation (Chace and Walsh 2006;Grimm et al. 2008;Gaston 2010;Faeth, Bang, and Saari 2011). Urbanisation has a homogenising effect on species richness, this effect being defined as the replacement of numerous specialist species with a few generalist species (Devictor et al. 2007). Similarly, Chace and Walsh (2006) highlight how the effects of urbanisation on avian richness depend on the type of species under investigation (e.g. native versus non-native).

Choosing survey sites
We randomly selected 20 University of Oxford colleges to survey from a list of 36. The spatial distribution of the 20 sites across the City of Oxford is shown in Fig. 1.

Organising bird surveys
We conducted bird surveys from 1 h before dawn to 1 h after. Multiple sources state it is optimal to listen to birds in urban settings during the 'dawn chorus' (Mayntz 2017; BTO-British Trust for Ornithology, n.d.), when avian song activity is at its highest and surrounding ambient noise and human activity are at their lowest. A.P.B. surveyed four sites every morning of the working week. Each 2-h morning slot was divided into quarters, so that each site received an approximately 30-min visit per week.
The 20 sites were divided into five groups due to the timing constraints imposed by distance between sites on efficient bird surveying (see Supplementary Section S1). The order in which blocks were surveyed during each working week was randomly allocated. The particular timing of a college bird survey was decided by randomly allocating sites to a single 30-min slot of each working week.

Conducting bird surveys
A.P.B. recorded bird species numbers along pre-established transects within each site. These transects were designed through reconnaissance visits in the 3 weeks leading up to the start of the survey and followed the sites' public footpaths and walkways. Avian species were mainly recorded via song, but also visually, whilst walking very slowly along the transect. Birds that flew above the surveyed site but did not land were not recorded.
Absolute species abundance was not measured in this study. Species identification took place in the field and was verified by reviewing recordings from a Sennheiser ME62 omnidirectional microphone and Marantz PMD661 digital recorder. All species recorded were identified by A.P.B. or A.G.Z. Nord University (Hemb 2018) and Xeno-Canto (Vellinga 2019) websites were used for training purposes and subsequent reference.

Calculation of independent variables
We calculated five independent variables: (1) site area, (2) site land use, (3) water cover, (4) vegetation cover and (5) land use surrounding sites. Site area (or the surveyed area, as was the case when the whole college grounds were too big to survey within 30 min, or where access to some quadrangles was restricted) was determined via the construction of layers of each surveyed area in Quantum Geographic Information System (QGIS version 2.14.14). Site area was calculated from these layers in hectares.
Site land use was determined by dividing site land use into six broad categories: (1) buildings and impervious cover, (2) mown grass, (3) meadow, (4) water, (5) shrubs, hedgerow and bushes and (6) tree cover, following Dale (2015). Proportions of each college habitat were ascertained via satellite imagery analysis on Google Earth Pro (2017), again following Dale (2015). A grid made of 10 Â 10 m 2 squares, generated using GE-Path Pro Software, was superimposed on a Google Earth Pro satellite image of each site. Each square was counted for the habitat type it contained and percentage cover of each site's habitat type was calculated. The determination of some habitat types such as meadow and tree cover proved to be difficult to ascertain solely from satellite imagery. Thus, on-the-ground information from bird survey visits was also incorporated when evaluating habitat type in each square.
We reduced the dimensionality of habitat type variables with principal component analysis (PCA). The PCA was run on R version 3.6.0 (R Core Team 2019) using the 'FactoMineR' package. Components with eigenvalues larger than 1 were included in the data analysis, as per the Kaiser-Guttmann criterion (Syms 2008).
The size of site water bodies was determined in the aforementioned site habitat diversity analysis. Vegetation cover was calculated by combining meadow, forest and shrub cover into one variable (in hectares).
The land use surrounding sites was determined using geospatial data surrounding each college within a 500-m radius. We chose a 500-m radius for two reasons: (1) a radius larger than 500 m meant that even sites distant from one another, such as Worcester and Magdalen, would have overlapping radii and therefore it would be difficult to discern differences in richness between these colleges due to surrounding land use; and (2) we considered a radius of less than 500 m inappropriate due to birds' abilities to fly over long distances and move easily between habitats.
The geospatial data for surrounding land use was collected from the Ordnance Survey (OS) and National Forest Inventory (NFI) data inventories (see Acknowledgements). We condensed the geospatial data using PCA in R with the 'FactoMineR' package. The variables included in the PCA were the (1) areas occupied (in hectares) by buildings [OS], (2) woodland [OS], (3) woodland [NFI], (4) greenspace [OS] and (5) water [OS] in the surrounding 500-m radius area around each college. Components with eigenvalues larger than 1 were included in the data analysis.
The geospatial data was collated into maps with QGIS software to provide a visual illustration of the surrounding land use of each site. Please see Supplementary Section S2 for the 20 maps generated.
We also interviewed college gardeners about the gardening practices applied in each site. The answers serve to record management practices driving habitat structure within sites and to document how encounters with avifauna influence garden management decision-making. These answers are not presented in the main text and can be found in Supplementary Section S3.

Sampling variables
We recorded the ambient temperature during each college visit with a digital thermometer. The relative time at which each college survey occurred during each morning (i.e. 1-4, corresponding to the order of site visitation) and the week during which each survey occurred (1-7) were also recorded.

Data analysis
Our data analysis had two approaches: (1) determining species richness (i.e. number) patterns and drivers via an informationtheoretic (IT) approach; and (2) investigating species composition patterns and drivers using Sorensen indices and non-metric multi-dimensional scaling (NMDS). Determining drivers of species richness across different sites From this point onwards, the terms 'species richness' refers to the total number of species that A.P.B. observed and recorded from all visits at each site. We built a global generalised linear model with a 'Poisson' distribution family and applied an IT/model-averaging approach to produce an average-weighted model following Grueber et al. (2011). The IT approach allows models to be ranked and weighted to provide a measure of relative support for competing hypotheses (Grueber et al. 2011).
Before constructing our global model to describe drivers of site species richness, we examined correlations among predictor variables to ensure no highly correlated pairs of variables (jrj > 0.8) were both included in the model. If pairs of highly correlated variables were found, then the variable that is simplest to measure was included and the remaining variable was excluded. Our global model included all variables mentioned in Section 'Calculation of independent variables' (unless correlations were found, see later). Interactions were not included due to the relatively low sample size (20 gardens). The residual plots of the global model are in Supplementary Section S4. We also checked for overdispersion in the global model and found none.
The global model was standardised [standardize() function, 'arm' package] before producing the candidate set of models via the dredge() function ('MuMIn' package). Standardising model coefficients allows for direct comparison of effect sizes between candidate models. The candidate models were limited to having a maximum of four variables per model due to the sample size.
Site area, site vegetation cover and site urbanisation were all highly correlated (see Supplementary Table S5.1). We therefore decided to remove site vegetation cover and site urbanisation from the global model and to keep site area. This is because we viewed site area as the underlying factor driving site vegetation cover and site urbanisation, as well as being the easiest variable to measure. Vegetation cover and site urbanisation were not converted to proportions-e.g. the proportion of vegetation cover-as these proportional variables were still highly correlated with each other and site area and showed the same trends with species richness as their non-proportional counterparts. The global model was 'dredged' to produce 31 candidate models, each with a maximum of four predictors. A 95% confidence set of 14 models was created from the larger set using the get.model() function ('MuMIn' package). An average-weighted model was produced from this smaller set using the zeroaveraging method, where a parameter estimate of zero is substituted into those models where the given parameter is absent, and the parameter estimate is obtained by averaging over all models in the top model set (Grueber et al. 2011). This method was chosen over the conditional average method-'where the parameter estimate for each predictor is averaged only over models in which that predictor appears and is weighted by the summed weights of these models' (Grueber et al. 2011)-due to the suggestion that the zero-averaging method is better when trying to identify which drivers have the strongest effect on the response variable, which is the aim of this study.
The summary() output of the average-weighted model and the confidence intervals of each estimate are shown in the Results. Confidence intervals were obtained via the confint() function in R, 'stats' package. Predictor importance is also included. Predictor importance is a measure computed from all models in the 95% confidence model set, using the Akaike weights summed for all models that contain a particular predictor variable (Burnham 2015). The predictor variable with the largest weight is estimated to be the most important of the predictors, whilst the variable with the smallest sum is estimated to be of least or no importance of all (Burnham 2015).
Assessing species community structure First, we calculated three Sorensen indices on a multi-site scale using the beta.multi() and beta.pair() functions. These simple statistics describe if species assemblages differ among sites. We calculated Sorensen dissimilarity, Sorensen turnover and Sorensen nestedness using presence-absence data with the 'betapart' package (Baselga and Orme 2012).
Second, we employed NMDS on relative abundance data (no. of species sightings per site divided by no. of sightings per site) as per Salako, Adebanji, and Glèlè Kakaï (2013) with the 'vegan' package in R (Oksanen et al. 2019). This analysis allowed us to investigate valuable species identity data that the model analysis omitted. We used NMDS fitted with the Morisita-Horn index to explain patterns and drivers of species composition of sites. We plotted stress values and a R 2 value for non-metric fit to assess the number of dimensions to use in our NMDS analysis. Scree plots indicated to use two dimensions, for which the stress value was 0.13 and the non-metric R 2 was 0.98. Significant drivers of observed patterns in our NMDS analysis were identified using the envfit() function ('vegan' package). The function fits environmental variables onto the ordination space and also identifies variables that have a strong correlation with the sites' coordinates in ordination space [computed by the metaMDS() function]. Significance levels for all correlations are calculated via permutations. No explanatory variables were omitted in this analysis as envfit() treats each explanatory variable separately.

Species richness estimates and species accumulation curves
We also calculated the predicted richness of each site (S Estimate ) and the variation around these predicted values (S Var ) using the Chao 2 equation for replicated incidence data [see Gotelli and Colwell (2011)]. Species accumulation curves were also plotted using ggplot2. These were used in conjunction with the S Estimate and S Var values to assess how well we sampled the avian assemblage of each site during the 7 weeks of bird surveys.
We found that 12 sites had species accumulation curves with asymptotes (Magdalen, Mansfield, Merton, Balliol, Exeter, Green Templeton, St. Anne's, Somerville, Christ Church, Corpus, Oriel and Pembroke). A total of 15 sites had species richness estimates that matched the observed site species total 63 species (Lady Margaret Hall, Green Templeton, St. Anne's, Corpus, Oriel, New, Christ Church, Trinity, Magdalen, Pembroke, Somerville, Exeter, St. Edmund's Hall, Lincoln and Mansfield). Based on this information, we concluded that our sites were well sampled for avian species richness. The species accumulation curves and species richness predictions can be found in Supplementary Sections S6 and S7, respectively.
Checking for spatial dependence We tested for spatial dependence in our response variables (species richness and two NMDS axes) by calculating the Moran's I statistic using the R 'spdep' package. We determined the centroid locations for each garden site and considered sites within 750 m of one another to be neighbours for the purposes of the analysis. We used the 'moran.mc' function to perform a permutation test on the Moran's I statistic, using 10 000 iterations. In all cases, the test returned small values of Moran's I and detected no significant spatial dependence (species richness: statistic ¼ 0.04, observed rank ¼ 7478, P-value ¼ 0.25; NMDS component 1: statistic ¼ À0.01, observed rank ¼ 6275, Pvalue ¼ 0.37; NMDS component 2: statistic ¼ 0.05, observed rank ¼ 7595, P-value ¼ 0.24). These results did not change appreciably if we set the neighbour threshold to either 500 or 1000 m, instead of 750 m. Due to a lack of spatial dependence in our data, we did not take steps to account for spatial dependence.

PCA results
Quantifying site land use The first principal component (PC1) on site building, water, meadow, tree, shrub and mown grass cover explained 47.8% of the observed variation and the second explained 25.4% of the variation. PC1 and PC2 were the only components with eigenvalues >1 (2.87 and 1.52, respectively) and thus both were included in the analysis as predictor variables (Supplementary Tables S8.1 and S8.2 and Fig. S8.1).
We interpreted principal components based on component loadings (see Fig. 2B and Supplementary Section S8). PC1 captures a gradient of site impervious surface cover and consequently is referred to in the rest of this article as 'site urbanisation'. PC2 captures further variation in site habitat composition after accounting for this impervious surface cover gradient, namely an axis from sites with less tree, water and mown grass cover and more shrub and meadow cover, to sites with generally more tree, water and mown grass cover and less shrub and meadow cover. This second principal component is consequently referred to as 'site habitat composition'.

Quantifying surrounding urbanisation
The first principal component on surrounding building, water, greenspace and woodland cover (PC1 A, with the 'A' label distinguishing these PCs from the PCs in Quantifying site land use) explained 63.3% of the observed variation and the second (PC2 A ) explained 22.1% of the variation. PC1 A and PC2 A were the only components with eigenvalues >1 (3.16 and 1.11, respectively) and thus both were included in the analysis as predictor variables (see Supplementary Section S9).
We interpreted the principal components as per Section 'Quantifying site land use' (see Fig. 2B). PC1 A captures a broad rural-urban gradient and is subsequently referred to as 'surrounding urbanisation'. PC2 A captures additional variation in surrounding habitat composition after accounting for this rural-urban gradient, namely an axis from areas with less woodland and greenspace and more water, to areas with generally more woodland and greenspace and less water. This second principal component is consequently referred to as 'surrounding habitat composition'.

Species assemblage patterns across sites
We observed 43 different species across all sites (S Absolute ). Worcester College had the highest total of species (30) and Mansfield College had the lowest (7). The three most frequently observed species were Turdus merula (Blackbird, 11.1%), Columba palumbus (Woodpigeon, 10.1%) and Erithacus rubecula (Robin, 9.85%). Supplementary Figure S10.1 details the relative frequency of species sightings across all sites.

Global model averaging
The average-weighted model identified site area as the only important and significant factor for predicting species richness in urban courtyard gardens (Table 1). Site area had the largest effect, with species richness increasing by a factor of e 0.48 ¼ 1.61 for every 1 SD (standardised unit) increase of area (95% CI ¼ 1.21, 2.16, importance ¼ 1). All other variables-water cover, site habitat composition, surrounding habitat composition and surrounding urbanisation-had non-significant effects and low relative importance (see Table 2).

Assessing community structure
The Simpson dissimilarity between sites (0.50) was almost twice as large as the Sorensen nestedness (i.e. similarity) between sites (0.28), producing a high Sorensen dissimilarity value of 0.77. This suggests there was greater species turnover than nestedness between our sites. Figure 4 shows the projection of sites and species on NMDS axes, and the predictor variables that were significantly correlated with these projections (a ¼ 0.05). Surrounding habitat composition was the only variable significantly correlated with the NMDS results. Table 3 shows the correlation values of the predictor variables and the significance of these correlations.

Discussion
Our dissimilarity indices and Figs 3 and 4 confirm that sites differed in species richness and community structure. We also confirm our hypothesis that these differences were driven by factors that are known in the literature to affect species richness, specifically site area for the total number of species observed in each site (Chamberlain et al. 2007;Gaston 2010) and surrounding habitat composition for community structure (Zhou, Fung, and Chu 2012;Xie et al. 2016). The effects of site vegetation cover and site urbanisation on species richness remain unclear, whilst it is apparent that site habitat composition and surrounding urbanisation do not have a significant effect on species richness. Interestingly, surrounding habitat composition and site water cover-although non-significant in the averaged model-are in the second-and third-best models of the 95% confidence set, respectively. These two variables are also significant and near-significant drivers of community structure, respectively, suggesting that they should not be overlooked as potential drivers of avian diversity.

Drivers of species richness
Our averaged model clearly shows that site area is the most important factor driving species richness. This fits the findings of previous studies (Chamberlain et al. 2007;Gaston 2010; Camacho-Cervantes, Ojanguren, and MacGregor-Fors 2018) and reinforces the importance of island biogeography in influencing biodiversity in urban settings (Davis and Glick 1978).
However, it is difficult to discern if site area alone is the factor affecting our observations. Whilst site area could be the prevailing factor in this study, the effects of site vegetation cover and site urbanisation on species richness should not be overlooked as there are a wealth of studies providing evidence for their effects on avian assemblages (Chace and Walsh 2006;Grimm et al. 2008;Oneal and Rotenberry 2009;Gaston 2010;Faeth, Bang, and Saari 2011). Sites with conifers in their grounds were noted to attract conifer-dependent species, such as Periparus ater (Coal Tit) and Regulus regulus (Goldcrest), further highlighting the potential effect of site vegetation cover on species richness. Alternative methods of quantifying vegetation and urbanisation could reduce the correlation among these Drivers of avian species richness and community structure in urban courtyard gardens | 5 three variables and improve our ability to identify their contributions to species diversity. For vegetation, one could include vegetation structure or composition, both of which have been shown to have a significant effect on urban avian richness in previous studies (Mills, Dunning, and Bates 1989;Ferná ndez-Juricic 2000;Hennings and Edge 2003;White et al. 2005;Gaston 2010;Xie et al. 2016). Future measures of site urbanisation could include the 3D structure, age and the material of site buildings. Reasons for incorporating these elements in future studies include how birds operate in a 3D environment, and observations of Cyanistes caeruleus (Blue Tit) and Parus major (Great Tit) using old buildings for nesting in New, Worcester and Merton College.
We recommend future investigations into drivers of urban avian assemblages to include these methods and improved variables.
Our averaged model did not find an effect of site water cover on species richness. Previous studies suggest that water does have a positive effect on species richness by attracting water-dependent species (Chamberlain, Cannon, and Toms 2004;Chamberlain et al. 2007). This is indeed true for some sites in our study sample. For instance, the Worcester College lake attracted Cygnus olor (Mute Swan) and Motacilla cinerea (Grey Wagtail), species unique to Worcester College. Furthermore, three sites with the highest species richness-  Figure 2: Variable factor maps (A, B) for site land use and surrounding land use, respectively. Variables factor maps illustrate the directionality and magnitude contributed by each variable towards the first two components (named here as dimensions 1 and 2). Sites' colour illustrates the block they belonged to Worcester, LMH and St. Catherine's-all contained water bodies (>0.02 ha). These sites were the only ones to attract waterdependent species such as Gallinula chloropus (Moorhen) and Branta canadensis (Canada Goose). Figure 4 shows many waterdependent species distributed around these three sites. Although the effect of water cover in our dataset was not significant in the zero-averaged model, this variable was present in 50% of models in the 95% confidence set. Therefore, future studies should continue to examine water cover as a potential driver of species richness.
Our model averaging analysis found that neither site habitat composition, surrounding habitat composition, nor surrounding urbanisation showed a consistent effect on species richness. However, the first two were present in several models in the 95% confidence set (relative importance values of 0.33 and 0.31, respectively), suggesting that further data collection might reveal additional patterns. Previous studies have found that species richness is affected by patch management (Calder et al. 2015), adjacent habitats (Chamberlain, Cannon, and Toms 2004;Chamberlain et al. 2007;Zhou, Fung, and Chu 2012;Xie et al. 2016) and surrounding urbanisation (Chace and Walsh 2006;Grimm et al. 2008;McKinney 2008;Gaston 2010;Faeth, Bang, and Saari 2011).

Drivers of community structure
The large Sorensen dissimilarity value (0.77) suggests that species assemblages differ substantially among sites. Our NMDS analysis confirms that species assemblages of some sites differ greatly and suggests that changing species composition is driven by surrounding habitat composition. Surrounding habitat composition is a principal component explaining variation in surrounding land use after accounting for the rural-urban gradient. It correlates most closely with NFI woodland, OS Water Cover and OS Greenspace (see Supplementary Table S8.2). These 'Importance' refers to the probability that a given predictor appears in the best model, and is estimated by summing the AICc weights of each best model (in the 95% confidence set) where the variable appears (Galipaud et al. 2014). Significant predictors are in bold. All predictors have been standardised on two SD as per Gelman (2008). a Are values from the 'full/unconditional average' of the summary output. The full summary table with conditional and unconditional averages can be found in Supplementary Section S11.     Figure 4: NMDS results. Both sites (number) and species (text) are plotted. Sites are numbered with a list indicating the name of each site in the top left. Drivers that have significant (a < 0.05) correlations with the NMDS coordinates of community structure are plotted in full arrows. The different colours refer to the different blocks that each site belonged to. Please see Supplementary Tables S13.1 and S13.2 for exact site coordinates and species coordinates facts would suggest that the proportion of surrounding woodland, water cover and greenspace is a significant determinant of the community structure of Oxford colleges. There are studies that support the influence of surrounding habitat patches on avian community structure such as Zhou, Fung, and Chu (2012), Chamberlain et al. (2007) and Lepczyk et al. (2017a), although these differ in habitat type and there appears to be no study supporting the effect of the relative proportions of surrounding habitats on avian community structure. We recommend that future studies continue to investigate the effects of the relative proportions of surrounding habitat on community structure. Additionally, further sampling would clarify the importance of the near-significant effects of water cover on community structure.

Future research
We highly encourage the continued sampling of Oxford college avian assemblages. Not only do they possess heterogenous environments, different species assemblages and belong to a working urban matrix, but they also bear a strong similarity to other European urban court (or cloister) gardens that make them a broadly relevant and replicable study sample for future research aimed at improving our understanding of the patterns and processes of urban avifauna. We also advance Fontana et al.'s (2011) recommendation to utilise multiple measures when studying urban biodiversity in light of the informative results provided by the IT, NMDS and Sorensen approaches we employed.
We suggest that future investigations include different avian diet types and guilds in their model analysis. For example, Faeth, Bang, and Saari (2011) showed that in urban areas the number of granivorous species tends to increase whilst the number of insectivorous species tends to decrease. Furthermore, different species classes, such as urban exploiters and non-native species, have been shown to respond positively to surrounding urbanisation (Mills, Dunning, and Bates 1989;Chace and Walsh 2006). Factoring in these different subsets of avian species as well as their abundances in the model analysis might clarify the effects of variables such as urbanisation on species richness in urban centres.
Other factors that ought to be considered in future studies of urban courtyard avifauna are the effects of abiotic factors, such as noise and pollution, on site species assemblages. For example, pollution has a strong negative effect on species richness (Sanderfoot and Holloway 2017), whilst anthropogenic noise reduces species richness by hindering the effectiveness of acoustic communication between birds in urban areas (Noam, Jonathan, and Yoram 2005;Habib, Bayne, and Boutin 2006;Gaston 2010;Perillo et al. 2017).
Subsequent studies on urban courtyard avian assemblages should consider building models of species richness that include gardening regimes, abiotic and biotic processes. The effect of biotic processes on species richness, such as intra-specific competition and avian communication, are key (Dabelsteen and Pedersen 1985;Gaston 2010). An appreciation of the effects of these biological processes on species richness already exists in some Oxford colleges. For example, the Worcester College gardener acknowledged the importance of ecological processes as determinants of avian species richness by employing a policy that focuses at the lowest trophic level (plants and invertebrates) before enacting strategies for higher trophic levels (small birds and predators).
It has also been suggested to the authors that the effects of site avifauna on the student inhabitants could be investigated. Fuller et al. (2007) suggest that the psychological benefits of greenspace increase with biodiversity, so the effects of different species richness levels on university students could be worth studying in order to see if exposure to different biodiversity levels improves psychological health. This is pertinent to British universities, where 'the number of first-year students arriving at university who report a mental health condition is now five times what it was 10 years ago' (Bewick and Stallman 2018). Admittedly, such investigation would have to collect and consider important social science data, such as individual histories, and the financial and welfare support provided by different colleges, all of which are likely factors affecting student mental health (Bewick and Stallman 2018;NHS 2018).

Overall implications
Our findings show that both local (in-site) and landscape (exsite) factors drive urban avian assemblages and that island biogeography is important in understanding urban avian diversity. These results suggest that management schemes designed to attract birds in urban areas would require coordination at multiple spatial scales and would support the incorporation of conservation principles into urban planning [as per Lepczyk et al. (2017a)]. Our interviews with stakeholders (i.e. site gardeners) showed that many managed their gardens to either attract birds directly or indirectly. Therefore, we highly encourage efforts to study and coordinate schemes promoting urban wildlife to take advantage of this model system of urban courtyard gardens.

Supplementary data
Supplementary data are available at JUECOL online.