First network analysis of interspeci ﬁ c associations of abyssal benthic megafauna reveals potential vulnerability of abyssal hill community

The distribution of organisms is related to both environmental factors and interactions between organisms. However, such associations between organisms across an abyssal megafaunal community have not previously been investigated at landscape scale because of a lack of positional data on specimens over such scales. We quanti ﬁ ed spatial distributions and investigated interspeci ﬁ c associations in benthic megafaunal communities in three contrasting habitats on the Porcupine Abyssal Plain, two on the abyssal plain and one on the ﬂ ank of a modest abyssal hill (~50 m above the plain). We used a Bayesian Network Inference Algorithm approach, which considers the ecosystem as a network, facilitated by robust positioning of specimens determined through seabed photography captured with an autonomous underwater vehicle. We found non-random intraspeci ﬁ c distributions of most morphotypes in all areas. The organisms in two interspeci ﬁ c networks on the abyssal plain had high connectance and link density, while the network at the Hill site was notable in the lack of inter-dependencies and highly dependent on one morphotype – Ophiuroidea. The reduced connectance of the hill community suggests that it is operating under a di ﬀ erent regime and potentially more vulnerable to perturbation than those on the plain. Interspeci ﬁ c dependencies on the abyssal plain occurred across broad taxonomic groupings, and were thought to be a result of similar relationships between pairs of organisms and the substrate, and competition for detrital resource. In addition, some intraspeci ﬁ c pairs changed dependency direction at di ﬀ erent scales. Our results suggest that the scales of inter- and intraspeci ﬁ c aggregation will be important considerations in the design of community assessments, and in spatial planning for their conservation.


Organism distribution in abyssal seabed habitats
The spatial distribution of organisms is related to environmental factors and biotic interactions between individual organisms in a habitat, and can influence ecosystem dynamics (Massol et al., 2011). Key to understanding the susceptibility or vulnerability of an ecosystem to change is the network of these interactions and associations between taxa (May, 1977), with communities of highly connected organisms (complex networks) more resistant to change than simpler networks (Rooney and McCann, 2012). These ecological networks could consist of only trophic relationships (food webs), but may also include habitat associations and inter-specific interactions, or a combination of all or some of these (Ings et al., 2009;Olff et al., 2009). The structure of these networks is crucial for understanding ecosystem dynamics, specifically individuals within and between taxa, with some resulting in aggregation (e.g. predation, symbiosis, high competition for common resource). Thus, the spatial heterogeneity of the deep-sea benthos likely influences the distribution of organisms in communities there.
The influence of organism interactions on the distribution of deepsea benthic megafauna is less well understood, with many dispersion studies limited by small sample sizes. Intraspecific aggregations have been observed at scales of 1-100 m, such as 'herds' of the holothurian Kolga hyalina (Billett and Hansen, 1982), aggregations of holothurians Peniagone sp. and Elpidia minutissima (Lauerman and Kaufmann, 1998), Paroriza prouhoi (Tyler et al., 1992b), echinoids Echinus affinis and Phormosoma (Grassle et al., 1975), xenophyophores (Lauerman et al., 1996), and ophiuroids (Lauerman and Kaufmann, 1998;Smith and Hamilton, 1983). Such aggregations were thought to facilitate the maximal use of deposited food and reproduction (Tyler et al., 1992a). Fewer studies examine interspecific interactions between abyssal megafaunal taxa. Serendipitous observations of direct associations between abyssal taxa have been made using seabed photography, such as predation on a polychaete worm by the anemone Iosactis vagabunda (Durden et al., 2015b). Photography has also captured interactions between organisms accessing the same food source, such as the hosting of the anemone Kadosactis commensalis by the holothurian Paroriza prouhoi (Bronsdon et al., 1993), augmented with data from collected specimens. Similar methods have also revealed interactions between organisms with a common method of exploiting environmental conditions, such as the attachment of anemones to glass sponge spicules (e.g. Amphianthus bathybium and Daontesia porcupina) where both organisms extend into the water column above the seabed to exploit suspended material in areas of enhanced currents. Despite these observational efforts, no studies have examined interspecific associations between organisms across a megafaunal community on abyssal plains and considered the impact on ecosystem dynamics.

Bayesian network inference -A tool for examining organism associations
One approach to understanding ecosystem dynamics is to consider the ecosystem as a network (e.g. Mitchell and Neutel, 2012). Ecological networks can be computed using discrete Bayesian Network Inference Algorithms (BNIAs) with morphotypes considered "nodes" and dependencies between them (where the abundance of one morphotype depends upon that of another) described as "edges" which link correlating morphotypes together (see Section 2.3) (Heckerman et al., 1995). BNIAs can infer network structures and non-linear interactions, and have been used extensively to reveal gene regulatory networks (Yu et al., 2002(Yu et al., , 2004, neural information flow networks and ecological networks (Smith et al., 2006, Milns et al., 2010 and more recently palaeontological communities (Mitchell and Butterfield, 2018). The structure produced by the BNIA reflects the associations caused by co-localizations rather than a specific interaction. For example, a negative dependency between two morphotypes could correspond to competitive exclusion or exclusive niches and positive dependencies could reflect trophic interactions, habitat associations or facilitations. While it is not possible to infer the underlying processes resulting in dependencies or connectance from our BNI analyses alone, the most likely processes can sometimes be inferred using biological observations. By using BNIA, direct dependencies between morphotypes can be detected, minimising auto-correlation between two variables. For example, if A depends on B which depends on C, there could be a correlation between A and C. However, this correlation would not represent an interaction or association between A and C, merely the two correlations between A and B and B and C, so the BNIA would report the edges A to B and B to C, and not report A to C. The BNIA approach enables only the realised dependencies to be found, ensuring only actual interactions and associations between morphotypes are found.

Aim
We quantify spatial distributions and investigate interspecific associations in benthic megafaunal communities on an abyssal plain using the BNIA approach. The approach is facilitated by robust positioning of specimens determined through seabed photography captured with an autonomous underwater vehicle, covering an extensive seabed area (Morris et al., 2014). We investigate the benthic invertebrate megafaunal communities at hectare scale in three areas at the well-studied Porcupine Abyssal Plain (PAP) Sustained Observatory (PAP-SO), a longterm time series site at 4850 m water depth in the northeast Atlantic (Hartman et al., 2012), two on the abyssal plain and one on the flank of an abyssal hill (mean water depth 50 m above the abyssal plain). The sediment conditions and benthic community on the hill were recently found to be significantly and substantially different to those in the areas on the plain; sediments in this area were coarser but organic carbon content was lower, and total numerical and biomass density, diversity, and suspension feeder density were higher on the hills than on the abyssal plain (Durden et al., 2020a). Differences between the two areas on the abyssal plain were more subtle; differences in sediment particle size were not detected, but the compositions of the two communities were significantly different, related to differences in the densities of the two most common morphotypes (Ophiuroidea and the anemone Iosactis vagabunda), and in terms of biomass related to differences in contributions from large holothurians. These three areas provide contrasting habitats in which to examine abyssal ecosystem dynamics in terms of intra-and interspecific associations of benthic megafauna. We consider these associations in terms of ecosystem dynamics, including organism interaction with the seabed substrate and feeding modes.

Megafauna in seabed photographs
Individual megafaunal specimens (> 1 cm in size) were located and identified as one of 75 morphotypes in randomized downward-facing seabed images captured by an autonomous underwater vehicle. The seabed images were captured along transects spaced at 100 m intervals forming a 'fine scale grid' pattern in each of three locations: PAP Central (the location of sampling for the 30-year time series), on the flank of an abyssal hill ('Hill' grid) and on the abyssal plain to the northwest of this hill ('North Plain' grid), as described in Durden et al. (2020a). Full details of image capture and processing are provided in Morris et al. (2014); in brief, image processing involved some colour correction, the removal of overlap between adjacent images and combining~10 successive images into a 'tile' (seabed area~14 m 2 ), and georeferencing the location of each tile based on the location of the vehicle (latitude/longitude) during image capture. A total of 3374 tiles (representing 44255 m 2 seabed) across the three grids were used in this analysis. Some morphotypes consisted of groups of species (i.e. genus, family-or order-level identification), because species-level identification from images is not possible or practicable (e.g. Ophiuroidea, Elpidiidae spp., Porifera).
Locations of identified megafaunal specimens were derived from the respective tile locations, then transformed to a two-dimensional flat surface using the R package geosphere (Hijmans, 2019). The area of each 1 km × 1 km fine-scale grid was split into cells, so that each contained the same length of photographic transect (i.e. the same sampling intensity); cell sizes were 100 m × 100 m (1 ha). The positions of the megafauna within the tiles was at the sub-10 m scale, so negligible within each hectare cell. Morphotype densities were calculated within each cell, so that the associations found describe interactions and associations that occur between the cells; these density data were used in the computation of both intra-specific distributions and the inter-specific association networks.

Analysis of non-random intra-specific distributions
In order to investigate the variation of intra-specific distributions in each photographic grid, chi-squared analyses of the counts per cell were performed in R (R Core Team, 2017). In this analysis, detectability of non-randomness is a function of abundance and strength of the pattern. That is, strong correlations can be detected with low abundances, and high abundances provide the power to find weaker correlations. Significance was reported at α = 5%. In order to check that relative abundance was not driving the non-random intra-specific distributions we performed a Kruskal-Wallis test between the abundances and the chi-squared p-values.

Bayesian network analysis of inter-specific associations
The Bayesian network (BN) consists of the best-fitting set of edges between nodes (i.e. the variables) which are found through a heuristic search mechanism Banjo (Smith et al., 2006) because finding the bestfitting BN is computationally intractable (Chickering, 1996). Assessment of a given BN is performed by calculating the Bayesian Scoring Metric, which calculates the probability that the network encodes of the statistical dependences of the observed data (Milns et al., 2010 and references within Appendix B). For ecological networks, "greedy" searches are an optimal technique for determining the network structure (Milns et al., 2010). Within Banjo greedy searches generate a random network, score how well the network fits the observed data and then randomly add edges to the network. If the score increases, then the edge is kept and another randomly added, and if the score decreases then it is removed. This process is repeated until there are no possible edge additions that increase the score and so corresponds to a local maxima within the search space. Because this best-fit BN is only a local maxima, the search needs to be repeated multiple times to ensure the global maxima is found. See Milns et al. (2010) and references within Appendix B for more details of the scoring and searching algorithms for ecological networks.
Bayesian network inference was performed in Banjo (Smith et al., 2006), with data preparation for Banjo (grouping and discretization) and statistical analysis in R (R Core Team, 2019). Further analysis of Banjo outputs, when required, used the functional language Haskell (Jones, 2003); scripts are available on Github (https://github.com/ egmitchell/bootstrap).

Data preparation -Discretization and contingency test filtering
The BNIA used requires discrete data, which ensures that correlations are not over-fitted to data noise and only the relative densities of each morphotype are important (Milns et al., 2010;Yu et al., 2004). Morphotype abundances by cell were split into three bins: zero counts, and counts below and above the median for the morphotype by site. For ecological species abundance datasets, three different bins have been deemed to be a good balance between maintaining information in the dataset (achieved with many bins) and increasing statistical power (achieved with fewer bins; Yu, 2005).
In order to avoid Type I errors introduced by high numbers of morphotypes with zero counts, we excluded morphotypes which were found in fewer than 33% of the cells. Generally, the result is that low density (less common) morphotypes are excluded; however, this method of exclusion could potentially mean that a morphotype with high abundance in a very limited area is excluded from analyses. This latter type of exclusion happened in the Hill grid, where two taxa with sufficient overall abundance to be included (Actinauge abyssorum and Peniagone spp.), were excluded from the network analyses because of the aggregation of these taxa over a very limited area. While there may be subtle inter-specific correlations that we are unable to detect, it is unlikely that these less common morphotypes would have a significant influence on community dynamics.
To further guard against Type I errors, we also used a method of contingency test filtering that removed from consideration an edge between two variables whose joint distribution showed no evidence of deviation from the distribution expected from their combined marginal distributions (chi-squared tests, p > 0.25; Milns et al., 2010). This threshold was used to ensure no chance of removing truly dependent dependencies, so that only artefacts, such as those found between high zero counts, were removed from consideration. These links were provided to the BNIA to exclude from consideration.

Bayesian network inference
The BNIA software used was Banjo v2.0.0, a publicly available Java based algorithm (Bernard and Hartemink, 2006;Smith et al., 2006). For details of the algorithm please see Smith et al. (2006) and Milns et al. (2010). The discretized data was input into Banjo which then generated a random network based on the input variables. A 'greedy search' was repeated 10 million times for each set of input data and the most probable network was then output. The maximum number of edges leading to a node was set to 3 to balance predictive power with overfitting limit artefacts (Milns et al., 2010;Yu, 2005).
To minimise bias from outliers, we bootstrapped at 95% level (Magurran, 2013) by randomly-selecting 95% of the total number of grids cells for each subsample and then finding the subsample network using Banjo. For each edge calculated, the probability of occurrence was calculated, and the resultant distributions analysed find the number of Gaussian sub-distributions using normal mixture models (Fraley et al., 2012). This probability distribution was bimodal for each dataset, which suggests that there were two distributions of edges, those with low probability of occurrence, and those highly probable edges. The final network for each area was taken to be those edges which were highly probable. The threshold for being labelled 'highly probable' depended on the network (as determined by the normal mixture modelling analyses): 55% for PAP Central, 49% for Hill and 52% for North Plain. The magnitude of the occurrence rate is indicated in the network by the width of the line depicting the edge.
The direction of the edge between nodes in the network indicates which node (morphotype) has a dependency on the other node (morphotype); this direction is indicated in the network by an arrowhead. For each edge, the directionality was taken to be the direction which occurred in the majority of bootstrapped networks. Where there was no majority (directional edges have a probability between 0.4 and 0.6) the edge was said to have bi-directionality, or indicated a mutual dependency; these are shown without arrows).
The influence score (IS) can be used to gauge the type and strength of the interaction between two nodes. If the IS = 1, this corresponds with a positive correlation: that is, a high density of morphotype 1 corresponds to a high density of morphotype 2, and the dependency arrow would point from morphotype 2 to morphotype 1. An IS of −1 corresponds to a negative correlation: a high density of morphotype 1 corresponds to a low density of morphotype 2. An IS = 0 does not mean there is no correlation between the two nodes, but rather that the correlation is non-monotonic so that the interaction is positive at low densities and negative at high densities (or vice versa). The mean IS for each edge was calculated for each site.

Intraspecific spatial distributions
In the grids on the abyssal plain, 14 and 23% of all morphotypes were included in the analysis (at PAP Central and North Plain, respectively), while 28% of morphotypes in the Hill grid were included ( Table 1). As these communities are dominated by few morphotypes (which were included), the included morphotypes accounted for the vast majority of specimens (92.6%; 92.5% and 85.3%, respectively). Seven included morphotypes were common to all grids; of these, three were mobile deposit feeders, one was a mobile predator/scavenger and three were sessile suspension feeders.
Total megafaunal density across all three study areas was nonrandom (Kolmogorov-Smirnov tests all p < 0.0001). The distribution of all fourteen morphotypes tested in the PAP Central fine grid were non-random (Table 2), while five of eight morphotypes in the North Plain grid were non-randomly distributed, and fifteen of nineteen morphotypes on the Hill randomly distributed. The three most common morphotypes, Iosactis vagabunda, Elpidiidae spp. and Ophiuroidea, were non-randomly distributed in all grids (Fig. 1). Morphotypes included in analyses of PAP Central and Hill grids with non-random distributions included Stalked tunicate, Oneirophanta mutabilis and Stalked crinoid. Cnidaria sp.9, Aphroditid and Porifera were also found to have nonrandom distributions at these sites, but random distributions in the North Plain grid. Amphianthus bathybium and Indeterminate (Indet.): Table 1 Total megafaunal morphotypes for the areas studied, morphotypes included in the analysis (those that occurred in > 33% of grid cells of 100 m × 100 m) and the number of connected morphotypes are those included in the Bayesian Network.  hydroid were distributed non-randomly at PAP Central, but randomly in the Hill grid.

Bayesian networks of interspecific associations
On the abyssal plain (PAP Central and North Plain grids), all but one of the included morphotypes in each grid were connected to others (92% and 88%, respectively), while in the Hill grid 77% of included morphotypes were connected (Table 1). These results were not reflections of relative abundancewe performed a Kruskal-Wallis test to compare the total abundances for each taxa to the p-values of the intra-specific chi-squared tests (from Table 2) between the spatial quadrats and found that there was not a significant correlation (p > 0.05).
Our Bayesian network inference found that all three areas had multiple dependencies between morphotypes (Fig. 2). The majority of the dependencies were positive (Table 3), with only positive dependencies for North Plain grid (6 of 6), 13 positive dependencies for Hill grid (of 15) and 17 positive dependencies (of 19) for PAP Central. There was only one negative dependency, in PAP Central. Where the interaction strength (IS) was zero, this indicates different behaviour at different densities, of which there were three in PAP Central and one in Hill grid. The highest mean IS was in the Hill grid network (0.6267) and the lowest in PAP Central (0.2989), with a mean IS for North Plain grid of 0.4960. Most dependencies were directional, with no mutual dependencies on Hill grid, two out of six mutual dependencies on North Plain grid and one mutual dependency on PAP Central.
The most complex area was PAP Central, exhibiting high connectance and link density (Table 3), highest dependency chain length of 11 morphotypes. The Hill grid was notable in the lack of inter-dependencies (Fig. 2, Table 3), with half the connectance of the grids on the abyssal plain. Both the PAP Central and North Plain grids had similar mean chain length, maximum length of dependency chain of 5 and connectance, but the North Plain grid differed in having fewer chains and lower link density than PAP Central grid.

Contrasting interspecific network structures in contrasting habitats
We found substantial differences in the interspecific associations of the benthic megafauna in the three sites studied. The networks of fauna in both grids on the abyssal plain are highly connected, while connectance in the network at the Hill is low. The dependency on a single morphotype (Ophiuroidea) within the Hill grid suggests that this network has little redundancy, so changes to the density of Ophiuroidea may spur substantial variation in the dynamics of this community. Connectance correlates with increased ecosystem stability and robustness independently of species richness (Dunne et al., 2002;Gardner and Ashby, 1970), through the increase of stabilizing processes such as of stabilizing feedback loops cf (Mitchell and Neutel, 2012) so that the significantly lower connectance of the Hill site suggests that it is less robust to perturbations that the plain sites.
The starkly different network structure at the Hill grid suggests different community dynamics from the abyssal plain community. Megafaunal communities at the Hill and plain sites have different community structures (Durden et al., 2020a), which are likely related to increased hard substrate availability and coarser sediments on the hill, a result of hydrodynamic conditions (Turnewitsch et al., 2015). In addition, the difference in the abundance and biomass k-dominance plots (a measure of the successional state of a community; Clarke, 1990;Warwick and Clarke, 1994) was similar for the two grids on the abyssal plain, but significantly different on the hill, suggesting that the community there was in a different successional state, or more recently disturbed (Durden et al., 2020a). The difference in successional state may be related to the differences in connectance and dependence found in the networks; specifically, the lower connectance of the Hill network and the suggestion that it is less robust to perturbations than the communities on the abyssal plain may be related to recent disturbance there.
It is also likely that the communities of the Hill and abyssal plains are subject to contrasting organizing forces with respect to available detrital resource. Deposit-feeding communities generally exist in trophically stable environments, which may be resource limited, while suspension-feeding communities exist in trophically more variable and resource-independent environments (Levinton and Kelaher, 2004). These generalities appear to apply to the PAP-SO megabenthos, where the community on abyssal hills is known to be more trophically diverse than on the plain (Durden et al., 2015;Durden et al., 2020a), and nearseabed detrital inputs in these two habitats have been estimated to be different in both quantity and its partitioning in the community . Furthermore, both top-down and bottom-up processes influence patch size in deposit feeding communities, such as the abyssal plain, with some local-scale patchiness likely introduced by localized inputs of detritus. These localized inputs of detritus may in turn influence the distribution of organisms with less mobility, while highly mobile deposit feeders alter the patchiness through grazing at distances greater than the existing patch size (Levinton and Kelaher, 2004). Unfortunately, detrital patch sizes at the scales examined in this Table 3 Properties of interspecific Bayesian network analysis of megafaunal morphotypes in each of the three photographic grids at the Porcupine Abyssal Plain. The 'mean IS' is the mean interaction strength of the connected network, excluding the zero IS values. Connectance is the number of realised interactions from the total possible number (May, 1973), while Link Density is the mean number of edges per node; both computed including the unconnected nodes.

PAP Central grid
The relatively high degree of connectance between morphotypes in the PAP Central grid involves connections between mobile and sessile organisms, and also between morphotypes that may feed on similar detrital fractions. The three most abundant megafaunal organisms on the abyssal plain, Iosactis vagabunda, Elpidiidae spp. and Ophiuroidea (Durden et al., 2015;Durden et al., 2020a), are major components of the network; at least one of these three morphotypes is involved in 12 of the 19 connections in this network. The base morphotype in most chains (with IS value > 0.3) is a sessile organism (Porifera, Cnidaria sp.7, Cnidaria sp.9 or Stalked tunicate), likely attached to hard substrate such as clinker or deposited ice-rafted dropstones (as previously observed in images; Durden et al., 2015), to which one of these three mobile organisms is dependent. This result suggests an aggregation of mobile morphotypes around the hard substrate and/or hard substrateattaching fauna, a behaviour observed for similar benthic megafauna around Antarctic dropstones (Chickering, 1996;Ziegler et al., 2017). An intermediary is Amphianthus bathybium (Fig. 1), which is itself associated with Porifera as it attaches to siliceous sponge spicules (Riemann-Zürneck, 1987). Sponges and hard substrate-attaching fauna provide three-dimensional structures in the deep sea, with which mobile fauna have been observed to be associated (Beaulieu, 2001;Buhl-Mortensen et al., 2010;Lacharité and Metaxas, 2017).
The mobile morphotypes in the network, including the holothurian Oneirophanta mutabilis, are surface deposit feeders, and thus individuals are likely co-located in order to feed on patches of deposited detritus. Some mobile surface deposit feeders alter their feeding behaviour related to the availability of detritus (Durden et al., 2020b), or movements to stay proximate to an area of detritus (Kaufmann and Smith, 1997). Environmental filtering related to detritus availability is likely important; aggregation suggests high competition between organisms for the shared detrital resource with low availability. The holothurians Elpidiidae spp. and Oneirophanta mutabilis differ in terms of trophic level based on δ 15 N values (Iken et al., 2001), the latter grazing approximately 5 times faster (Durden et al., 2019). The trophic levels of specimens of several Ophiuroidea are similar to each of Elpidiidae spp. and O. mutabilis (Iken et al., 2001). The hemisessile Iosactis vagabunda moves at a slower speed, moving burrow locations approximately every 20 days (Durden et al., 2015b) to exploit deposited detritus around the burrow. Its trophic position is considerably higher, likely because of its occasional predation; I. vagabunda occupies a similar trophic position to the Stalked tunicate based on δ 15 N values (Iken et al., 2001), which may explain their association in the network.
The negative interaction between Cnidarian sp.7 and Stalked crinoid indicates spatial separation of these suspension-feeding morphotypes. This negative interaction is likely due to niche filtering since suitable hydrodynamics may be facilitating their suspension feeding. It may be simply related to substrate texture, with Cnidarian sp.7 utilising soft substrate and Stalked crinoid requiring hard substrate.

Hill fine-scale grid
The network of the Hill grid is characterised by the major role of Ophiuroidea, which was the most abundant morphotype at the Hill, and more abundant there than at the sites on the abyssal plain (Table 2). Ophiuroidea at the Hill were more strongly associated with Elpidiidae spp. and Iosactis vagabunda that at PAP Central. These latter two morphotypes were less abundant on the Hill than on the abyssal plain, but were each connected in the network to sessile morphotypes (Porifera, Cnidaria sp.9 and Tunicata), as on the plain. Ophiuroidea at the Hill was also associated directly and individually with 5 other sessile suspension feeding morphotypes and three other mobile ones, some of which it was associated via Elpidiidae spp. at PAP Central (e.g. Stalked crinoid, Aphroditid). One major chain was unbroken from the PAP Central network: from Ophiuroidea to Elpidiidae spp. to Porifera. The 'Ophiuroidea' morphotype consists of all ophiuroids (including at least Ophiocten hastatum and Ophiomusium lymani), and interspecific differences in feeding habits are likelyfour different morphotypes studied by Drazen et al. (2008) and Iken et al. (2001) had different δ 15 N values. Thus, many of the associations with Ophiuroidea may be related to its feeding modes, including surface deposit feeding, or predation/ scavenging, in similarity with the polychaete Aphroditid. Such flexibility of feeding mode may be advantageous in the hill environment, where proportionally less organic matter is deposited on the seafloor, as some suspended material is likely removed by suspension feeders before deposition .

North plain grid
The morphotypes in the North Plain grid were highly connected, as in the other abyssal plain grid, despite the network being less complex. As in the PAP Central grid, strong associations between the mobile deposit feeding Iosactis vagabunda, Elpidiidae spp. and Ophiuroidea are likely related to these organisms aggregating around patches of detritus. Similarly, the network analysis suggests that these mobile morphotypes are found to aggregate with sessile organisms (in this case, Amphianthus bathybium and Stalked tunicate Fig. 2). I. vagabunda was found to have mutual dependencies with A. bathybium and Elpidiidae spp., in contrast to the relationships of these two pairs at PAP Central. At PAP Central, I. vagabunda is dependent on A. bathybium and on Elpidiidae spp., but the latter relationship is via dependencies on Ophiuroidea and Oneirophanta mutabilis (which is not present at the North Plain). Thus, the dependency of Elpidiidae spp. on Ophiuroidea at the North Plain is the reverse of the dependency between these morphotypes at PAP Central. These differences in the network of the North Plain from that of PAP Central suggest that the relative decrease in the density of I. vagabunda and the increases in the densities of Ophiuroidea and Elpidiidae spp. alter the relationships between these morphotypes and potentially the community function there.

Contrasting interspecific distributions
Despite most common morphotypes having non-random distributions across both abyssal hill and plain habitats, some morphotypes had different interspecific distributions between these habitats. The anemone Amphianthus bathybium was non-randomly distributed at PAP Central, but randomly distributed in the other two grids (Table 2, Fig. 1). Porifera, its host, was similarly non-randomly distributed at PAP Central, and also in the Hill grid, but randomly distributed at the North Plain. Again, where the distributions were determined to be random, this randomness may be a true signal or the abundances of these morphotypes and the weakness of the signal may be too low to detect any patterns. Note though that because there is no significant correlation with abundance and non-random spatial behaviour (p = 0.4748) any non-random signal for low abundance morphotype would have to be weak to not be detected. Similarly, the lack of detection of any dependency between these two morphotypes at the Hill and North Plain grids, in contrast to those found at PAP Central, may also be related to relative low densities of one or both morphotypes at these sites. However, both live and dead Porifera provide substrate to which A. bathybium attaches (Riemann-Zürneck, 1987), and the distribution of such stalks was found to be random at another abyssal site (Beaulieu, 2001), complicating the analysis of this morphotype association pair.

Spatial considerations
Organisms may aggregate or associate differently in different areas and at different scales, a situation which occurred at the abyssal plain sites for some organism pairs. At PAP Central, the dependency direction between Cnidarian sp.7 with Peniagone spp., a minor pair in terms of abundance, was scale-dependent. By contrast, at the North Plain site, both dependencies of the dominant Iosactis vagabunda (with Amphianthus bathybium and with Elpidiidae spp.) changed orientation with scale. This suggests that at some scale, Elpidiidae spp. may be highly dependent on others, with the network being highly centred on this organism.
The sizes of the photographic grids and the cells determined the maximum and minimum of spatial scales at which interspecific distributions and intraspecific associations would be detected; the number of grid cells was a balance between having sufficient available data per cell, and sufficient cells for analysis. Thus, our analyses would not detect within-cell (< 100 m) patchiness and aggregations, nor nonrandom behaviour which occurs over the kilometre scales of the finescale transects. Patchiness and aggregations that occur between 100 m and 1 km could be detected. These limitations, in combination with the observed densities and community structure at this scale, determined the included morphotypes and detectable network interactions.
Only a small fraction (13-28%) of the morphotypes at each site were sufficiently abundant to be included in the intraspecific dispersion and interspecific network analyses. This was related to the density dominance of the megafaunal assemblages by a few organisms at this scale, particularly in the two grids on the abyssal plain (Durden et al., 2020a). The included morphotypes are generally small to medium in size, particularly in the two abyssal plain grids; the larger morphotypes, such as Psychropotes longicauda, Paroriza prouhoi, Benthodytes spp. and Benthothuria sp., are present at low densities (Durden et al., 2020a), precluding their inclusion in the analysis. Thus, their distributions remain unknown and the networks represent only a portion of the potential morphotype associations within each assemblage. Two large holothurians, Pseudostichopus aemulatus and Molpadiodemas villosus, were included in the analysis of the Hill grid, but were not found to have non-random distributions or associations with other morphotypes there. As these omitted morphotypes are generally large, mobile fauna, likely with larger ranges given their fast locomotion (Durden et al., 2019), they may have non-random distributions or associations at spatial scales larger than that of this study.
The exclusion of low density morphotypes in combination with the spatial scale of the study resulted in some known associations of megafauna at the PAP not being identified. For example, pairing in the large holothurian Paroriza prouhoi, and non-random co-occurrence between epi-and basibionts, such as P. prouhoi acting as host for the anemone Kadosactis commensalis (Bronsdon et al., 1994;Bronsdon et al., 1993), or the anemone Actinauge abyssorum residing on siliceous sponge fibres and polychaete tubes (Riemann-Zürneck, 1986). These associations involve direct interactions, at scales of less than 1 m, and thus would not have been detected even if the large holothurians (including P. prouhoi) were included in the analysis. The distributions and associations of morphotypes that exist at low densities at this scale, and those that aggregate at small spatial scales, and their impacts on community structure, could be assessed with more high density seabed photography.

Temporal considerations
This analysis provides a snapshot of the community as photographed in summer 2012. However, megabenthic communities in abyssal habitats are known to vary in activity and composition at interannual and seasonal time scales Kuhnz et al., 2014;Durden et al., 2020b), with likely alterations to the patterns of intraspecific and interspecific aggregation. Interannual differences in the structure of the community has been previously found to alter its function in terms of carbon flow at this site  and another abyssal plain site (Dunlop et al., 2016). Densities of Amperima rosea, a small holothurian and component of Elpidiidae spp., and Ophiuroidea have fluctuated significantly interannually at the PAP-SO site (PAP Central here); at the time of this survey, Elpidiidae spp. occurred at densities in between the boom and bust densities found previously Billett et al., 2001). As one of the most common megafaunal morphotypes, significant fluctuations in its density could alter interspecific aggregations in the community. Seasonal reproduction in deep-sea echinoderms (Tyler et al., 1982) may involve periodic intraspecific aggregations of holothurians (e.g. Tyler et al., 1992b). Seasonal detrital inputs may induce aggregations of organisms exploiting the deposited detritus, since many abyssal deposit feeders are selective feeders (FitzGeorge-Balfour et al., 2010;Ginger et al., 2001). Patch selectivity in deposit feeding holothurians (Uthicke and Karez, 1999) is exemplified by the loop or run-and-mill behaviour of Oneirophanta mutabilis, a taxon presumed to feed on patches of detritus (Kaufmann and Smith, 1997). The periodic input of detritus also induces seasonality in the seabed surface activity of some burrowing organisms (Durden et al., 2020b). Thus, the community structure and inter-and intraspecific aggregations likely also vary seasonally.

Conclusions
We found a complex network of interspecific associations in two abyssal plain communities, and a contrasting simple network on a modest abyssal hill (~50 m above the plain), where the community was highly dependent on a single morphotype. The significant differences between these networks suggest that the communities are operating under contrasting regimes, likely with different environmental filtering (related to hard substrate and organic matter availability) or different organizing forces (related to differing faunal function, such as through flexible feeding modes). These results also suggest the vulnerability of the megafaunal community on the abyssal hill to perturbations, as opposed to the relative robustness of the abyssal plain communities. Thus, disturbances to the sedimentary environment, such as from sediment deposition as a result of deep-sea mining (Jones et al., 2017), or alterations to the deposition of detrital material, which is likely to result from climate change (e.g. Jones et al., 2014), may impact the structure and dynamics of communities in these two habitats differently. Furthermore, the spatial scales of heterogeneity and aggregation for organisms of interest will be important to the robust design of any assessment of megafaunal communities in abyssal benthic habitats, including baseline environmental assessments (e.g. for Environmental Impact Assessments; Clark et al., 2017;Durden et al., 2018) and assessments of ecosystem change (e.g. environmental monitoring), or spatial planning for their conservation (e.g. Jones et al., 2018). These results also demonstrate the importance of considering the interactions of communities as a whole when assessing ecosystem vulnerabilities.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.