Intraspecific Abundance–Occupancy–Patchiness Relations in the Intertidal Benthic Macrofauna of a Cool-Temperate North Sea Mudflat

The macroecological variable of patchiness (Lloyd’s Ip index of patchiness, P) has recently been shown to be related inter- and intraspecifically to those of abundance (numbers m−2, A), and occupancy (% occurrence in samples, O) in lower latitude seagrass macrofaunas. For the first time in higher latitudes, intraspecific relationships between three spatial variables were investigated in the intertidal mudflat macrobenthos of the Scolt Head barrier island, southern North Sea (53° N, 01° E). Sampling was conducted between early July and late September 2009–2013 using 710-µm mesh for sample processing. Strong positive interspecific A-O and negative interspecific P-O and P-A relationships were present. Two of the most numerous and widespread assemblage components, however, occurred with effectively constant occupancy (Peringia ulvae, 100%, and Tubificoides benedii, 93%) across the whole 20-ha locality and therefore could not show intraspecific relationships of occupancy with other macroecological metrics. These two apart, only one other dominant species failed to show a significant positive intraspecific A-O relationship; no species showed significant P-A relations of any form; and only two showed the negative P-O ones that have been described elsewhere. The intraspecific A-O patterns appear to contrast with those of an earlier study at another North Sea locality (the Dutch Wadden Sea), although differences are more apparent than real, but the Scolt Head fauna showed fewer intraspecific P-O and P-A relations than those characterising similar circumstances in the two lower-latitude localities previously investigated. Neither developmental mode nor variation in local abundance appears to influence these patterns. A-O-P relations therefore seem widespread but may be subject to latitudinal modification.


Introduction
The interspecific (i.e. between-species) macroecological spatial abundance-occupancy (A-O) relationship is much researched and well-attested across a wide range of habitat types (Verberk et al. 2010;Frisk et al. 2011;Dallas et al. 2019;Miranda and Killgore 2019). Indeed, it has been argued to be one of the most fundamental patterns in ecology (Gaston et al. 2000). The macrobenthic faunas of estuarine and other sheltered coastal mudflats and seagrass beds, for example, seem characterised by such relationships (Foggo et al. 2003;Barnes 2020aBarnes , 2021a. Values of occupancy associated with the macrofaunal elements of intertidal seagrass beds in the subtropical Moreton Bay, Australia, and in the warm-temperate Knysna estuarine bay, South Africa, can be predicted with great accuracy solely from local abundance data (Barnes 2021a) using the general A-O expression proposed by He and Gaston (2003). It has also been demonstrated in those same southern-hemisphere seagrass meadows that the dispersion pattern of the invertebrate macrofaunal species, as assessed by the I p patchiness index of Lloyd (1967) (P), is also negatively related to both their interspecific abundance and occupancy (Barnes 2021a), so that the three macroecological metrics appear interlinked, at least in a seagrass habitat.
The comparable situation with respect to such relationships intraspecifically (i.e. within species), however, appears less straightforward, although it has generally been considered that most A-O patterns are similar intra-and Communicated by Judy Grassle.

3
interspecifically. Confirmatory evidence for the intraspecific A-O relationship has come largely from studies of vertebrates, especially fish and birds (Gaston et al. 1999;Zuckerberg et al. 2009;Frisk et al. 2011;Faulks et al. 2015), but although its occurrence is also regarded as widespread in other taxa (Gaston et al. 2000), various exceptions have been noted (Blackburn et al. 2006). All the dominant members of the Moreton Bay and Knysna seagrass systems mentioned above that displayed strong interspecific A-O relationships also displayed clear intraspecific ones (Barnes 2021b), but most members of similar macrobenthic soft-sediment assemblages of the cool-temperate Wadden Sea have been recorded as not doing so (Bijleveld et al. 2018). All intraspecific P-O relations within the South African and Australian seagrass species were negative, but only significantly in two-thirds of the species. Whilst most P-A ones were also negative, less than half of them were significant (Barnes 2021b). One third of the component seagrass species at those localities showed no significant P-O or P-A correlations. Blackburn et al. (2006) considered that the strongest A-O relationships are likely to be found in marine habitats, but they have relatively rarely been investigated there. To date, the only sites that have supported studies of both interand intraspecific A-O-P relationships are the Moreton Bay and Knysna seagrass beds, whilst the only site with studies of the magnitude of equivalent intraspecific A-O relations amongst the coastal macrobenthos is the Dutch Wadden Sea (Bijleveld et al. 2018). The present work therefore sought to extend this very limited range of intraspecific coastal studies by investigating a northern-hemisphere site similar to the Wadden Sea site: the north Norfolk coast in the UK which also comprises (or before reclamation comprised) a chain of sand-dune barrier islands semi-isolating a small 'Wadden Sea' in their lee that is largely exposed during low tide (Steers 1960). It is known that the macrofauna of a seagrass system there shows the typical strongly positive interspecific A-O relationship, together with weaker negative interspecific P-O and P-A correlations (Barnes 2021a). Hence, it is possible further to explore the generality of A-O-P relations by testing at that general site three aspects of the intraspecific patterns: (1) whether Bijleveld et al.'s (2018) 'atypical' A-O results in the Dutch Wadden Sea are replicated in another southern North Sea locality, (2) whether P-O and P-A macroecological patterns are similar in the cool-temperate northern hemisphere to those in the lower latitude, southern hemisphere sites, and (3) whether those individual species that do not seem to conform to the apparent norm of strong positive A-O and less-strong negative P-O patterns also show a common suite of systematic or ecological features.

Methods
Twelve historical datasets of the intertidal macrobenthic fauna of the 20-ha Cockle Bight mudflat on the sheltered lee side of the Scolt Head Island National Nature Reserve, southern North Sea (52° 59′ 05″ N, 00° 40′ 35″ E) (Fig. 1), mostly collected between the summer months of early July to late September in 2009-2013, were available for analysis. Since the mudflat's Zostera noltei seagrass beds and the adjacent unvegetated Arenicola marina dominated muddy sand show essentially the same fauna with only minor differentiation (Barnes 2014), in this study, the two habitats were not distinguished, and the data analysed included both those from the Z. noltei beds and unvegetated flats of Barnes and Ellwood (2011) and Barnes (2014). This faunal similarity, of course, runs contrary to the general 'vegetated versus unvegetated' paradigm of differential abundance and species richness across the two habitat types (Hemminga and Duarte 2000), but it is typical of the southern North Sea (Polte et al. 2005;Bouma et al. 2009), possibly because of the short-leaved and open nature of the Z. noltei beds there. The qualifications for inclusion of a species in the analyses were its presence (1) in densities of > 1 ind. in any single core sample in at least six individual datasets, (2) in overall abundances of > 125 m −2 , and (3) with overall occupancies of > 20%. Ten such species were present, together comprising 99% of the total assemblage of individuals and 83% of the total species occurrences. Analysis of intraspecific patterns was therefore confined to the more widespread and abundant species, not only because confidence in the nature of macroecological metrics assesssed from relatively few samples is low (e.g. Green 1966) but also because A-O relationships at least may be different in rare versus widespread species (Borregaard and Rahbek 2006;Bijleveld et al. 2018), regression slopes being shallower and coefficients of determination (R 2 ) weaker in the rare category (Freckleton et al. 2006;Buckley and Freckleton 2010).
All samples from the intertidal Cockle Bight macrobenthic assemblage were obtained using the same methodology, involving series of cores each of nominal 60-mm diameter (0.0027-m 2 area) and 50-mm depth, taken at 1-m intervals across each site, with a minimum of 20 core samples per site and a total of 320 samples. Sampling strategy is always a compromise between number of replicates and number of organisms present in unit sample, and the core crosssectional area was chosen as one that would yield > 100 individuals per sample (see Barnes and Ellwood 2011). Core depth was selected because most (> 98%) individual macrofauna in such circumstances are known to occur in the uppermost 3 cm (see, e.g. Klumpp and Kwak 2005;Ysebaert et al. 2005) and expected numbers of the more deeply burrowing species present in the area (Arenicola marina and Heteromastus filiformis) were in any event unlikely to approach 1 core −1 . Samples were collected during low tide and were gently sieved ('puddled') through 710-µm mesh on site. This procedure collects the smaller (mostly < 5 mm) and more numerous members of the macrofauna that constitute the large majority of invertebrate biodiversity (Bouchet et al. 2002;Albano et al. 2011), though not the meiofauna nor sessile species attached to the seagrass leaves. Warwick et al. (2006) have shown that different patterning rules may apply to meiofauna and macrofauna, and likewise Davidson et al. (2004) and Leopardas et al. (2014) to sessile species. Retained material from each core was transported immediately to a local laboratory on the island and was there placed in a 30 × 25-cm white tray in which the living fauna was located by visual examination until no further invertebrates could be observed. Individuals were identified to species level under a binocular microscope, with all organismal nomenclature here being as listed in the World Register of Marine Species (www. marin espec ies. org) (Accessed April 2021).
All abundance data are presented as numbers m −2 and for individual species calculation of mean densities at a given site included unoccupied samples (i.e. zero values). Occupancies are percentages of the total samples at a given site in which a species was present. In conformity with the comparable inter-(Barnes 2021a) and intraspecific (Barnes 2021b) data presented earlier for equivalent macrobenthic assemblages in the warm-temperate Knysna estuarine bay and subtropical Moreton Bay (Barnes 2021a,b), dispersion pattern was ascertained by spatial point pattern analysis of count data using Lloyd's I p index of patchiness (Lloyd 1967), I p = [1 + 1/k], where k is the dispersion parameter of the negative binomial distribution, i.e. = [1 + (v-m)/m 2 ], 'm' is the mean abundance across samples and 'v' is the associated variance. This index has been demonstrated to yield equivalent results to those of the spatially explicit Moran's spatial auto-correlation index for equivalent intertidal macrobenthos (Barnes and Hamylton 2019). It is also independent of sample size over a wide range of area, provided that the individuals are randomly distributed within a patch and that the patches are large relative to sample size (Lloyd 1967;Myers 1978). Granted that the core area was 0.0027 m 2 , it seems unlikely that macrofaunal patches were smaller than that. Again, to enable direct comparison with the earlier localities mentioned above, occupancy data are also presented in (log) percentage occurrence form, rather than as logit transformations of proportional data: transformation of occupancy values has no effect on rank A-O or P-O correlations or on the pattern of relative position of individual data points across logarithmic space (and values of R 2 of power-law relationships are not significantly different: see Barnes 2021b).
The Wadden Sea data of Bijleveld et al. (2018) are available in the form of one value for each of A and O per species per year per tidal basin at http:// doi. org/ 10. 5281/ zenodo. 11203 47, and these were re-analysed for correlations between intraspecific abundance and occupancy.
All calculations were carried out in Microsoft Excel for Mac 16.50 with the StatPlus:mac Pro 7.5.0 add-on, or via PAST 4.04 (Hammer et al. 2020); curves were fitted using KaleidaGraph 5.0; and correlations between metrics are given as Spearman's rank correlation coefficient, C. Occupancy data were, however, logit-transformed to test homogeneity of slopes of the P-O relationship by ANCOVA. Powerlaw scaling coefficients (exponents) are abbreviated below to β and the associated coefficients of determination of the power-law relationships to R 2 . Statistical probability of departures of values of Lloyd's I p from 1 were determined by Monte Carlo testing with 9999 replicates each with a random distribution of points across samples. To avoid occasions where there might be possible confusion between patchiness (P) and probability, the latter is abbreviated to 'Prob.'

Results
A total of 38 species was encountered of which 12 were singletons, the remaining 26, totalling 45,666 individuals, occurring at a mean density of 52,530 m −2 (range 31,455-78,400 m −2 ). This assemblage generally conformed to the sheltered and often estuarine southern North Sea fauna described in detail by Wolff (1973) and Ysebaert et al. (2003). Together, the 26 species displayed both highly significant interspecific A-O and interspecific P-O and P-A relationships (Fig. 2). The overall interspecific A-O relationship ( Fig Values of neither assemblage abundance nor its patchiness showed any significant correlations across the 12 sites with other macroecological metrics (i.e. N 0 species density, N 2 species diversity, or J evenness) (all C < 0.43, regardless of sign; all Prob > 0.17).
Of these almost 46,000 individuals, 79% were of the one species, Peringia ulvae, which occurred with an occupancy of 100% (see Electronic Supplementary Material) at a minimum observed density of 2,200 m −2 (≅ 1 per 0.00045 m 2 ) and a maximum of > 130,000 m −2 . Since the value of occupancy for P. ulvae was spatially invariant, intraspecific A-O or P-O relationships were not possible. Occurring at such overwhelming densities, P. ulvae could conceivably negatively influence the spatial pattern of its sympatric fauna, but the only negative correlation of its numbers with those of any other species was with Ampharete acutifrons (Prob = 0.001) and with a C = − 0.18 in that respect the correlation was very weak. No negative correlation in the entire dataset was other than in the 'none to very weak' category of Moore et al. (2018), as were all positive correlations with the exception of a weak though highly significant one (C = 0.43; Prob << 0.00001) between the two canalipalpatan polychaetes Pygospio elegans and Fabricia stellaris.
Intraspecifically, only two of the 10 most abundant and widespread species failed to display a highly significant A-O relationship in addition to the Peringia ulvae case above. One of these was Tubificoides benedii which showed almost invariant values of occupancy (occurring in a mean 93% of samples, and in 100% of samples at six of the twelve sites) and yet nevertheless achieved a close to significant A-O relationship at Prob = 0.07 (Table 1, Fig. 3). Slopes of the relationship differed significantly between the seven species that displayed a significant A-O correlation (ANCOVA F = 3.29, Prob(same) = 0.008). No species displayed a significant P-A correlation (Table 2, Fig. 4). Only two species, T. benedii and A. acutifrons, showed a significant intraspecific P-O relationship (Table 2, Fig. 4), both being negative and with a moderate level of correlation, although values of R 2 were relatively low (< 0.5). Slopes of the P-O relationships in the two were not significantly different (ANCOVA F = 4.27, Prob(same) = 0.053). There were no evident differences between the biology of T. benedii and A. acutifrons on the one hand and all the other species on the other that might account for their differential P-O responses. Neither were there any significant correlations within any species between the degree of variation in abundance across sites (as assessed by the relevant coefficient of variance) and the magnitude of the A-O correlation (C = 0.12, Prob = 0.8) or of the associated value of R 2 (C = 0.35, Prob = 0.4) or magnitude of the P-O correlation (C = − 0.02, Prob = 0.9).
The Wadden Sea dataset of Bijleveld et al. (2018) is not directly comparable with the present ones, but re-analysis of their data for each of the 16 dominant macrofaunal species at that locality shows a significant positive intraspecific correlation between A and O in 14 of the 16 species, highly so in 11 of them (C = 0.50-0.86, Prob << 0.00001) including Peringia ulvae. Only in Hediste diversicolor and Maculomangus tenuis was there no statistically significant intraspecific A-O correlation (C < 0.25, Prob > 0.15). Comparisons between the frequencies of P-A-O relations observed across the four localities thus far investigated in any of these respects (Scolt Head, Wadden Sea, Knysna and Moreton Bay) are shown in Table 3.

Discussion
These results suggest that relationships between macrobenthic patchiness, abundance, and occupancy continue into lower latitudes but are less marked than those demonstrated previously in the warm-temperate and subtropical southern hemisphere. They also have a bearing on the conclusions of Bijleveld et al. (2018) on abundance-occupancy relations in the equivalent macrobenthic assemblage in the Wadden Sea.
The Dutch Wadden Sea macrofauna investigated by Bijleveld et al. (2018) is similar although not identical to that at the southern North Sea site investigated here. Both have Peringia ulvae as by far the most abundant species, but only three others of the 16 individual Wadden Sea species they studied (Cerastoderma edule, Abra tenuis and Limecola balthica) were amongst the 10 dominant Scolt Head ones, although another six were also present. At face value, the Fig. 2 Interspecific A-O-P relationships of macrobenthic species in Cockle Bight mudflat (mean abundance m −2 vs mean % occupancy of each species), with logarithmic axes and fitted power-laws: A A-O pattern; B those of P-A and P-O (C = Spearman's rank correlation coefficient, and P = probability of C value; β = scaling coefficient of fitted power law, and R 2 = coefficient of determination) results on intraspecific A-O relationships at these two sites appear very different, however. The situation in which two members of the Scolt Head fauna, P. ulvae and Tubificoides benedii, were present in all or virtually all of the 320 individual 0.0027-m 2 samples and therefore could not participate in relations involving variation in occupancy is clearly very unusual, and excluding those two, only one species (that present in the fewest qualifying datasets) failed to display a significant positive A-O relationship on Scolt Head. In the Wadden Sea, on the other hand, Bijleveld et al. (2018) concluded that although positive spatial A-O relationships were found, overall these were weak (slopes of the relationships being all < 0.33) and only manifested at the spatial scale of the entire 15,000 km 2 Wadden Sea. Even at the largest spatial scale, they considered that only seven of the 16 species showed reasonably strong A-O relationships. Changes in such relations across scales have been rarely investigated elsewhere, but in the Knysna estuarine bay, similar though interspecific A-O relationships have certainly been observed in intertidal macrobenthos across different spatial scales, i.e. at landscape (entire estuarine bay), regional (significantly different faunal regions of the bay) and local (individual sites), at least within the bay's total area of 16 km 2 (Barnes 2019). The entire scale of the studied Knysna bay, however, would fall within the 'local' scales of Blackburn et al. (2006) and Bijleveld et al. (2018), and data on variation in intraspecific A-O correlations across scales are currently only available for the Wadden Sea.
Differences between the present intraspecific results and those of Bijleveld et al. (2018), however, may be more apparent than real. The research on Scolt Head (together with that at Knysna and Moreton Bay referred to above) sought whether there were significant correlations between observed levels of A and O at a given locality, and if so the form of these relationships. Bijleveld et al. (2018), however, were examining the extent to which the magnitude of one of these variables explained the other, via assessing associated values of R 2 and differences between observed and predicted values. Statistically significant correlation does not necessarily indicate that most of the variation is thereby explained, although a mean value of 0.90 for the correlation between A and O for each of the seven species showing significantly correlated values on Scolt Head is not only high but as high as those known from the southern hemisphere. Re-analysis of Bijleveld et al. data indicated that although positive correlations were less marked than those on Scolt Head, values for those species common to both localities showed a highly significant positive correlation (with a mean value of 0.63), and, overall, 14 of the Wadden Sea species did indeed show a significant, though moderate, positive correlation of A and O, even though the associated values of R 2 indicated that one accounted for only 30-53% of the variation in the other.
Spatial variation in abundance is a commonly described feature of intertidal soft-sediment macrofaunal species and assemblages (Bassindale and Clark 1960;Azovsky et al. 2000;Dittman 2002;Zajac 2004;Schenone and Thrush 2020), even across apparently uniform stretches of habitat (Kraan et al. 2009;Barnes and Laurie 2018) and over small distances (Reise 1979;Barnes and Hamylton 2019). Nevertheless, a relationship of such patchiness to A-O patterns has only recently been proposed (Barnes 2019) and to date P-A-O correlations have only been investigated at three widely separated localities. The increasing latitudinal sequence of these, i.e. Moreton Bay → Knysna → North Sea, shows decreasing levels of intraspecific P-A and P-O relations, as well as increasing assemblage abundance (< 3000 to > 50,000 m −2 ), decreasing overall species density (> 200 to < 40), and increasing mean occupancy (< 12 to > 20%)  (Barnes 2021a and present paper). The level of overall assemblage patchiness, however, was the same across all three systems, at Lloyd's I p values of 1.15-1.18 (Barnes 2020a(Barnes , 2021a. Although fewer species on Scolt Head showed significant intraspecific A-O correlations than at the lower latitude sites, the magnitude of those that were significant was comparable, i.e. C = 0.90 and R 2 = 0.69 on Scolt Head versus mean values of C = 0.90 and R 2 = 0.81 across the Moreton Bay and Knysna localities (Barnes 2021a). The equivalent values for the significant P-O relationships at Moreton Bay and Knysna were lower than for A-O at mean values of C = (−)0.70 and R 2 = 0.49 (Barnes 2021a), and the only two datapoints available for P-O from Scolt Head were both lower still.
Overall assemblage abundance appears remarkably constant at the Moreton Bay locality, measured at the same time each year from 2011 to 2018 having a range of only 2300-2600 m −2 (Barnes and Hamylton 2019); the present data from Scolt Head, in contrast, had a mean of 52,530 and a range of 31,455-78,340 m −2 . Spatial and/or temporal variation in abundance correlates with equivalent variation in occupancy under positive A-O relationships (He and Gaston 2003), and such variation has in turn been correlated with fluctuations in local recruitment, for example in the Wadden Sea (Bijleveld et al. 2018). Hence, attempts have been made to correlate this intraspecific spatial variation with juvenile dispersal rates (Holt et al. 2002;Foggo et al. 2007;Webb et al. 2009). Abra tenuis, for example, deposits its eggs in the sediment (Holmes et al. 2004) and shows poorer dispersal than the otherwise ecologically similar tellinoid Limecola balthica that spawns into the overlying water (Honkoop and van der Meer 1997). Correspondingly, A. tenuis showed the weaker A-O relationship in the Wadden Sea (Bijleveld et al. 2018). On Scolt Head, however, both species had highly significant A-O relationships and L. balthica displayed the lower value of R 2 . The only species not to show a significant intraspecific A-O relationship on Scolt Head (other than the omnipresent P. ulvae and T. benedii) was indeed the non-planktonic brooder Fabricia stellaris (Rouse 1995). But Tubifex costatus, Pygospio elegans and Ampharete acutifrons are all species with very limited dispersal (Birtwell and Arthur 1980;Price and Warwick 1980;Kesäniemi et al. 2012;Thonig et al. 2017), and they all displayed the same relationship as those with dispersive planktonic larvae. Further, the available database for F. stellaris was smaller than those of all the other species and this may have been responsible for failure to achieve statistical significance. Thus, there was no clear separation of A-O response based on reproduction mode among the Scolt Head species, nor indeed among those in Moreton Bay or Knysna (Barnes 2021b).
Neither was there a correlation between degree of variation in local population density (as assessed by the coefficients of variation in abundance) and P-O relationship across the Scolt Head sites, nor across the species at the three contrasting localities so far investigated (C = − 0.02, Prob = 0.9). There was also no evidence that members of the same group of invertebrates with, as far as is known similar ecologies, reacted in a similar fashion across sites A mean abundance m −2 and O mean % occupancy, across qualifying datasets C rank correlation coefficient, Prob probability of C value (with probabilities < 0.05 shown in bold), β scaling coefficient (power law exponent) where a significant correlation exists, R 2 coefficient of determination of the power-law fit, n number of qualifying datasets All values of I p (P) indicate significantly patchy dispersion (all Prob < 0.004)  Bijleveld et al. (2018) pointed out, that although the interspecific A-O relationship may be sufficiently strong to permit accurate estimation of the level of O solely from that of A, that is not the case with the intraspecific one; and P-O and P-A relationships are weaker still. Nevertheless, within intertidal estuarine and coastal mud-and sandflats, the macrofauna may display very close interspecific A-O relationships, and many if not all the dominant components of these assemblages also display corresponding intraspecific A-O and P-O patterns. The fact that significant intraspecific A-O-P relationships do occur, even if they are of lesser magnitude than interspecific ones, therefore requires explanation. So, perhaps to an even greater degree, do exceptions to the apparent generalities, since a number of the dominant species do not show some, or potentially any, such patterns. An explanation for this is not yet evident, and neither is whether the potential latitudinal gradient in P-O and P-A relationships shown in Table 3 is real. Further investigation of these exceptions, of the relationships between patchiness of the macrofauna and of their microphytobenthic food base (see, e.g. Daggers et al. 2020), and of more localities in total are obviously needed before explanation of the nature of intraspecific A-O and P-O responses is possible.
Author Contribution RSKB conceived, designed and executed this study and wrote the manuscript. No other person is entitled to authorship.
Data Availability All data generated or analysed during this study are included in this published article and its supplementary information files.

Conflict of interest
The author declares no competing interests.
Ethical Approval All applicable international, national and/or institutional guidelines for sampling, care and experimental use of organisms for the study were followed and all necessary permissions and approvals were obtained in respect of the original collections of the data.

Consent for Publication
The author consents to publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.