1 Introduction

A central issue of community ecology is to understand why some species, but not others, are observed in a biological community on a certain spatial and temporal scale (Morin, 2011). Although niche differentiation is considered as a key factor in the coexistence of plant species (Tilman, 1984; Grubb, 2010), the important role of niche differentiation in species coexistence is also expressed by some ecologists in varying degrees of skepticism. Functional traits of species are more important than environment or evolutionary relatedness in determining the plant species co-occurrence pattern (Soliveres et al., 2014). Generally, some possible mechanisms for community assembly have been proposed, including the effects of competition, complex relationships between species niches and environmental heterogeneity, and favored and/or unfavored species combination on the co-occurrence patterns of communities (Diamond, 1975). Additionally, different factors may play a major role in structuring fish communities, such as predation and environmental factors (Englund et al., 2009; Giam and Olden, 2016).

A null model analysis is an important statistical test based on the randomization of ecological data or random sampling from a known or imagined species distribution (Gotelli, 2000), and it has been extensively applied to examine the co-occurrence patterns of different animal or plant communities. For example, null model analysis indicates that competitive exclusion is strong among distantly related species in co-occurrence patterns for Bornean vertebrates (Beaudrot et al., 2013); co-occurrence patterns of species composition are random and unstructured for ectoparasites assemblage (Gotelli and Rohde, 2002); and competitive interaction plays an important role in community assembly for the forest trees (Rooney, 2008). Null models are used to study the species co-occurrence of fish communities in lakes in Ontario (Jackson et al., 1992), stream fish communities (Yant and Karr, 1984; Peres-Neto, 2004), demersal fishes in a tropical bay in southeastern Brazil (Azevedo et al., 2006), and fish communities in the litorals of three lagoons of the Orinoco River floodplain (Gabriela and Nirson, 2017).

Temperate coastal waters or bays are generally important spawning, feeding and nursery grounds for many fish and invertebrate species, which are crucial to the recruitment of fish resources and support many local fisheries (Wafar et al., 1983; Chen, 1991). Haizhou Bay, located in the western Yellow Sea, is an important spawning, nursery and feeding ground, and fishing ground for many commercially important fish species in the Yellow Sea (Chen, 1991). Previous studies focused on species diversity, taxonomic diversity and community structure of fish assemblages and their relationships with environmental variables in the Haizhou Bay (Su et al., 2013; Wang, 2013). The fish community in the bay is divided into the northern and southern assemblages characterized by different species composition in spring and autumn (Zhang et al., 2018). With fisheries management moving towards ecosystem-based method, it is more and more important to illuminate the spatial structure of fish community and species coexistence pattern (Roafuentes and Casatti, 2017; Ruppert et al., 2018). Ascertaining coexistence mechanism of fish species is a prerequisite for better understanding of fish community dynamics, and sustainable utilization and conservation of the fish resources in the coastal waters.

Detecting the spatio-temporal variations in co-occurrence patterns of fish communities is helpful to explore the community assembly mechanism, and provides some theoretical supports for fisheries management in the Haizhou Bay. In this study, null model analysis was applied to analyze the presence-absence matrices to examine the fish species co-occurrence patterns in spring and autumn over different years in the temperate bay. The goals of this study are to: 1) test whether fish species are assembled randomly in the Haizhou Bay; 2) explore the spatial and temporal variation of fish species co-occurrence patterns; 3) find significant species pairs structuring the fish assemblages.

2 Materials and Methods 2.1 Study Area

Haizhou Bay is located in the western Yellow Sea (Fig. 1). The depth increases from west to east in the bay. The hydrographic conditions are mainly influenced by the Yellow Sea Coastal Current along the southern Shandong Province and northern Jiangsu Province, the Yellow Sea Warm Current, and the Yellow Sea Cold Water Mass (Lin, 1989). The bay was divided into two parts at a 34°50'N boundary (Fig. 1), which primarily corresponded to waters of the northern and southern fish assemblages defined by Zhang et al. (2018).

2.2 Data Collection

The data were collected from the bottom trawl surveys of fish species conducted in May and October from 2013 to 2015 in the Haizhou Bay. The surveys followed a stratified random design in which the strata were defined mainly according to the oceanographic, geological, regional, and biological characteristics of the Haizhou Bay (Xu et al., 2015). The survey area was divided into 10'×10' sampling grid cells, and 18 sampling grids were randomly selected in each survey with 2, 4, 2, 7, and 3 sampling stations in stratum A, B, C, D, and E, respectively (Fig. 1). All the bottom trawl surveys were conducted using a 220 kW otter trawler during the daytime with the towing speed at about 2.5 knots and the trawling duration of about 1 hour on average. The open width of the sampling net was 25 m and the mesh size for the cod end was 17 mm. The fish samples were iced, brought back to the laboratory and identified as different species. The catches per haul in weight and number at each sampling station were recorded. The protocol for the stratified random survey design was described in details by Xu et al. (2015). The environmental factors including bottom water temperature, salinity and depth were measured simultaneously by CTD during the fish resources surveys.

2.3 Data Analysis 2.3.1 Co-occurrence patterns of fish assemblage

Fixed-fixed (FF) null model assumed that both the row and column sums of the original matrix were fixed, so that differences in the occurrence frequency of each species (row sums) and differences in the number of fish species per site (column sums) were preserved (Connor and Simberloff, 1979). The FF null model with low Type I error probability had a good power for detecting nonrandomness of species co-occurrence (Gotelli, 2002).

The C-score under FF null model algorithm was used to quantify the co-occurrence patterns of fish assemblage with presence-absence matrix data in the Haizhou Bay in spring and autumn from 2013 to 2015. According to its good power with algorithm described above, the C-score under FF null model had good power in dealing with noisy data (Gotelli, 2002).

The C-score quantified the average number of Checkerboard units that could be found for each species pair, which was calculated as follows (Stone and Roberts, 1990; Gotelli and Entsminge, 2001; Gotelli, 2002):

${C_{ij}} = ({R_i} - S)({R_j} - S), $

(1)

where R_i is the number of occurrences (row total) for species i, R_j is the number of occurrences for species j, and S is the number of sampling sites where both species occur. Let N = n (n-1)/2 be the number of species pairs for n species, then the C-score is:

$C = \sum\nolimits_{i < j} {{C_{ij}}/N} .$

(2)

We used 30000 swaps per randomization and null distributions were generated from 10000 randomization of the focal matrix (Lehsten and Harmand, 2006; Lavender et al., 2016). The null hypothesis will be accepted when the observed value falls within the 95% confidence interval of distribution of 10000 times simulated values; the null hypothesis will be rejected when the observed value falls outside the 95% confidence interval.

2.3.2 Identification of significant species pairs

The whole matrix analysis only gave the pattern in average for co-occurrence of fish communities. Further analysis was performed to find out which particular pairs of species showed significantly positive versus negative associations. The Z-score was used to identify the significance of the C-score of each species pairs. Z-score below -2 or above 2 indicated significant aggregation or segregation (P < 0.05) (Ulrich, 2008; Gotelli and Ulrich, 2009; Echevarría and González, 2017).

The null models analyses were performed using the package EcosimR 0.1.0. (Gotelli et al., 2014). Pairwise species associations were calculated using PAIRS (Ulrich, 2008).

3 Results 3.1 Co-Occurrence Patterns of Fish Communities 3.1.1 Co-occurrence patterns of fish species for the whole fish assemblage

The observed C-score was significantly larger than the expected one by chance and the null hypothesis was rejected (P < 0.05), showing that fish species displayed negative co-occurrence patterns in spring and autumn from 2013 to 2015 in the whole bay (Fig. 2). The observed C-score based on the combined data of three years also showed that fish species displayed negative co-occurrence patterns both in spring and autumn, which confirmed the patterns for the whole fish assemblages in each year (Fig. 3).

3.1.2 Co-occurrence patterns of fish species for the northern and southern assemblages

The observed C-score was significantly larger than the expected one by chance and the null hypothesis was rejected (P < 0.05), which indicated that fish species displayed negative co-occurrence patterns for the northern fish assemblages in spring from 2013 to 2015. The C-score showed that fish species were randomly distributed in 2013 and 2014, and showed negative co-occurrence pattern in 2015 for the southern assemblage in spring (Fig. 4). The observed C-score was not significantly different from the expected by chance, and the null hypothesis was accepted for the northern and southern fish assemblages in autumn over three years (Fig. 5).

In spring, the observed C-score based on the combined data of three years was significantly larger than the expected one by chance and fish species exhibited negative co-occurrence patterns for the southern and northern fish assemblages, which confirmed the patterns in each year. In autumn, the C-score based on the combined data of all years showed that fish species in the southern and northern fish assemblages were randomly distributed (Fig. 6).

3.2 Species Pairs with Significant Co-Occurrence Pattern

The numbers of fish species varied from 43 to 50 in spring and from 48 to 51 in autumn in the whole bay. In spring, 35-41 fish species in northern assemblages and 24-32 species in southern assemblages occurred from 2013 to 2015. In autumn, the numbers of fish species were 38-45 for northern assemblage and 30-40 for southern assemblage over three years. There were 30 warm temperate fish species, 12 cold temperate species, and 22 warm water species in 2013. In 2014, 27 warm temperate fish species, 9 cold temperate species, and 23 warm water species were caught. There were 22 warm temperate fish species, 10 cold temperate species, and 26 warm water species in 2015. There were more cold temperate species in the northern assemblage (Table 1).

More significantly positive species pairs than negative ones were observed for the fish communities in the whole bay in both spring and autumn in different years (Table 2). Most of the segregated species pairs were composed of cold temperate species, such as Fang's blenny (Pholis fangi), Tanaka's snailfish (Liparis tanakai), fat greenling (Hexagrammos otakii) and Pacific sandlance (Ammodytes personatus), and other fish species (most species belong to warm water species or warm temperate species) for the fish communities in the whole bay. All these cold temperate species showed segregated patterns with red tongue sole (Cynoglossus joyneri) in spring over different years (Table 3). The numbers of segregated pairs of cold temperate species with other species were less in autumn than those in spring for fish communities in the whole bay.

Many significantly positive species pairs were observed in the northern assemblage, which indicated a segregated pattern. However, few segregated species pairs were found in the southern assemblage (Table 4). Most segregated species pairs in the northern assemblage were formed by cold temperate species (including Fang's blenny, Tanaka's snailfish, fat greenling, Pacific sandlance) with other species. Six segregated species pairs were selected in the northern assemblage in spring of 2013, and all the pairs were formed by Pacific sandlance with others (Table 5).

4 Discussion

This study showed that species co-occurrence patterns of fish assemblages were non-random and highly structured, and the fish species co-occurrence patterns kept relatively consistent and stable in different seasons and years in the whole Haizhou Bay. However, fish species co-occurrence patterns were different for the northern and southern assemblages in different seasons and years.

This study provided new evidence supporting the preponderance of non-random structure for fish communities that have been observed in different waters (Jackson, 1992; Echevarría and González, 2017). Species co-occurrence of fish assemblages were non-random and structured in the whole bay, which was similar with those of birds and mammals (Gotelli and Rohde, 2002; Pitta et al., 2012). However, it is hard to determine the extent to which the negative association is due to competitive species interactions or isolation of species (Gotelli and Rohde, 2002). Species co-occurrence patterns might vary with spatial scales (Mcnickle et al., 2017). Hoeinghaus et al. (2006) showed that fish community in temperate systems tends to show a nonrandom pattern at a regional scale, whereas it shifts to randomness at fine scale. The spatial and temporal variations of fish community structure should be incorporated into the analysis of fish species co-occurrence patterns.

Despite the tendency of negative co-occurrence patterns for the fish species, several random and aggregated species pairs were found for the fish assemblage in the Haizhou Bay. Therefore, non-random structure might result from non-random associations of only a subset of species (Krasnov et al., 2011). Cold temperate species (Fang's blenny, Tanaka's snailfish, fat greenling, Pacific sandlance) were most frequently segregated with several other species such as red tongue sole, filamentous shrimpgoby (Cryptocentrus filifer) in spring (Table 3). However, cold temperate species played less roles in significant species pairs, and most segregated species pairs were formed between cardinal fish (Apogonichthys lineatus), jewfish (Johnius belangeri) and other fish species in autumn in the bay. In the northern assemblage in different years, some cold temperate species were segregated with other species, whereas low number of pairs with significant co-occurrence patterns were found in the southern assemblage. These cold temperate species were strongly related to water temperature, which were identified as significant species pairs in the northern assemblage. It was proved that the environmental variability between the northern and southern parts of the bay was a factor resulting in the variation in the spatial distribution of cold temperate species like Fang's blenny and warm water species like red tongue sole in the bay (Wang, 2013). It showed that the water temperature and substrate were different in northern and southern areas of Haizhou Bay, and the northern assemblage was characterized by most cold temperate fish species. The distribution of Tanaka's snailfish, Fang's blenny, fat greenling, and anglerfish (Lophius litulon) shows negative correlations with temperature and positive correlations with depth in different seasons (Wang, 2013).

Environmental variables such as bottom water temperature and salinity show distinct spatial heterogeneity in spring and autumn in the Haizhou Bay (Wang, 2013; Su et al., 2013). The north-eastern waters is characterized by relatively low temperature and high salinity, and the southwestern waters has characteristics of high temperature and low salinity in the Haizhou Bay in spring. The salinity in autumn shows a similar pattern with that in spring in the bay. Bottom water temperature has no obvious pattern in autumn although it is still low in northeastern waters of the bay (Wang, 2013). The environmental heterogeneity on the bay scale has an effect on the spatial distribution of fish species with different ecological traits. For example, cold temperate species such as Pacific sandlance, Tanaka's snailfish, Fang's blenny and fat greenling widely distributed in the northern part of Haizhou Bay, whereas temperate or warm water species including red tongue sole and filamentous shrimpgoby, widely distributed in the southern part of the bay (Wang, 2013). In addition to the bay-scale heterogeneity, high and low spatial variations in environmental variables exist in the northern and southern parts of bay, respectively, which might result in negative co-occurrence of fish species in the northern waters and randomly distribution in the southern waters. The distributions of significant species pairs in the northern and southern assemblages were mainly related with the water temperature differences in the northern and southern waters of the bay. In terms of the distribution of species pairs and their habitat characteristics, environmental heterogeneity could play an important role in structuring the co-occurrence patterns for the northern and southern assemblages in the Haizhou Bay.

Gotelli and Ellison (2002) suggested that harsh environments (habitats) are the primary filtering factors for assembly rules. Azevedo et al. (2006) used unconstrained and constrained null models to examine the co-occurrence of fish communities in different seasons and three zones during summer in a tropical bay, and stated that habitat segregation could explain the pattern of reduced co-occurrence. Environmental filtering is by far the most important process structuring the fish communities (Azevedo et al., 2006). The importance of environmental filtering in driving the fish community structure is also proved by co-occurrence of fish species in tropical streams (Peres-Neto, 2004). In the Haizhou Bay, environmental heterogeneity might be a factor affecting the co-occurrence patterns of northern and southern fish assemblages. The differentiation of co-occurrence patterns between the northern and southern assemblages can lead to different implementation schemes for the conservation and management of fish resources in the northern and southern parts of the bay.

In summary, fish communities in the whole bay showed non-random species segregation patterns in the spring and autumn over different years. However, the northern and southern assemblages showed different patterns compared with the fish community in the whole bay. Positive species pairs were not always more than negative ones in significantly co-occurring species pairs. Environmental conditions had an effects on the co-occurrence patterns of fish communities in the bay. The deep understanding of co-occurrence patterns of fish communities would be helpful in the ecosystem-based management in the Haizhou Bay. Further study will be needed to explore to what extent the environmental variables or biotic interaction plays a role in affecting the co-occurrence patterns of fish species in the Haizhou Bay. These results expose the need to further explore the co-occurrence patterns in temperate coastal waters or bays.

Acknowledgements

We are grateful to all scientific staff and crew for their assistance in data collection during the surveys. This work was funded by the National Natural Science Foundation of China (No. 31772852), the Fundamental Research Funds for the Central Universities (Nos. 2015 62030, 201612004), and the Public Science and Technology Research Funds Projects of Ocean (No. 201305030).