Multivariate analysis of the spatial species diversity of demersal fish assemblages in relation to habitat characteristics in a subtropical national park, Taiwan

To understand the spatial species diversity of demersal fish assemblages in Taijiang National Park (TJNP) of Taiwan, fishes from 44 demersal trawl hauls and environmental data were collected in the nearshore and offshore areas of TJNP from April 2016 to May 2019. In total, fishes of 47 families, 84 genera, and 113 species were recorded. The nearshore and offshore demersal fish assemblages in TJNP exhibited significant variability in species composition assessed via beta diversity. Using distance-based redundancy analysis, we demonstrated that bottom depth and substrate type were significant explanatory variables of spatial species diversity and identified three habitat types (I: shallow soft bottom; II: deeper soft bottom; III: deeper bottom with mixed sand and gravel substrates). The nearshore assemblage was characterized by type I, where Tarphops oligolepis (flounder), Trachinocephalus myops (snakefish), and Liachirus melanospilos (carpet sole) dominated in terms of abundance. The offshore assemblage was characterized by either type II or type III because differences in substrate types among sampling sites were noticeable. At the offshore sites characterized by a deeper soft bottom (type II), Johnius distinctus (croaker), Cynoglossus kopsii (shortheaded tonguesole), and Coelorinchus formosanus (Formosa grenadier) predominated. In contrast, the westernmost sampling site, characterized by type III habitat, exhibited relatively high Shannon indices, and Scorpaena miostoma (scorpionfish), Urolophus aurantiacus (sepia stingray), and Parabothus taiwanensis (lefteye flounder) predominated. Our results provide the first baseline information on the environmental characteristics and spatial species diversity of demersal fish assemblages in TJNP and have implications for biodiversity conservation in existing spatial management areas.


Introduction
Demersal fish assemblages are formed as a result of the environmental preferences of species; that is, gatherings of species dependably co-occur in constrained subareas and in habitually expansive zoogeographical areas. Nonetheless, coastal marine environments can be highly diverse due to different substrate types and hydrographic settings. For the spatial heterogeneity and habitat use of coastal-offshore demersal fishes, water depth, substrate type, temperature, and salinity have frequently been demonstrated as the most imperative environmental drivers (Araújo et al. 2002;Moore et al. 2009). Understanding the spatial patterns of demersal fish biodiversity across regional gradients, such as depth, is one of the major goals in marine ecosystem studies (Zintzen et al. 2012(Zintzen et al. , 2017. However, only a few studies on beta diversity patterns in marine demersal fish assemblages have been conducted (e.g., Anderson et al. 2013;Zintzen et al. 2017;Vega-Cendejas and Santillana 2019). Moreover, knowledge of the small-scale spatial heterogeneity of demersal fish assemblages and environmental settings in coastal-offshore environments is important for creating management and conservation strategies operating at a range of scales (Chang et al. 2012;Hewitt et al. 2015;Amezcua et al. 2019).

3
Many studies have indicated that the community structure of demersal fish assemblages changes along a depth gradient (e.g., Labropoulou and Papaconstantinou 2004;Fitzpatrick et al. 2012;Stephenson et al. 2020). Labropoulou and Papaconstantinou (2004) reported that depth is a significant factor affecting the diversity patterns of demersal fish assemblages. Using a baited remote underwater stereo-video system, Fitzpatrick et al. (2012) found that the assemblage structure and diversity of demersal fishes are associated with habitat/depth categories, that is, habitats from shallow water to deeper water. Stephenson et al. (2020) reported that the demersal fish assemblages off New Zealand are differentiated mainly by hydrographic conditions such as depth and temperature. Determining major patterns of assemblage variability along various habitat gradients (or types) provides useful information for marine spatial planning as well as for developing monitoring schemes that can adequately assess changes in fish population abundances (García-Rodríguez et al. 2011;Shephard et al. 2011;Amezcua et al. 2019;Moura et al. 2020;Stephenson et al. 2020). Scientific surveys of marine demersal fish assemblages with data on environmental settings serve as essential and necessary foundations for subsequent quantitative research. Comprehensive research on the community structure and spatial species diversity of demersal fish assemblages can provide not only species records but also qualitative classifications of species into ecological (functional) or zoogeographical groups for conservation and management purposes (Amezcua et al. 2019;Moura et al. 2020;Stephenson et al. 2020).
Taijiang National Park (TJNP; https:// www. tjnp. gov. tw/ Eng/), located in southwestern Taiwan, was established in 2009. The legislation associated with TJNP focuses on the protection of wetland biodiversity, the characteristic salt and fishing industries and the history of settlement. In TJNP, the lowest monthly average temperature (17.6 °C in 1981-2010; www. cwb. gov. tw) occurs in January on land, and here, we use the term "subtropical national park" for TJNP. TJNP has important wetland habitats for wintering black-faced spoonbills (Platalea minor). In addition to the terrestrial domains, TJNP also comprises two marine areas, i.e., "marine existing use areas 1 and 2 (Fig. 1)." Marine area 1 (hereafter referred to as the nearshore area) covers the nearshore waters from the coastline of the terrestrial domains to either the 20-m isobaths or a 3-nm offshore distance. The nearshore area has a shallow seabed, with the deepest depth of ca. 85 m at its central-western boundary. Marine area 2 (hereafter referred to as the offshore area) encompasses the region from the central-western boundary of the nearshore area to the southeastern boundary of South Penghu Marine National Park (https:// www. marine. gov. tw/ home-en). The offshore area is located across the Penghu Channel, which has submarine valley morphology and is considered a scour furrow that probably originated from erosion (Huang and Yu 2003). The depths of the offshore-area seabed range from ca. 60 to 200 m, making this area deeper than the nearshore area. The offshore area was established because it was the main water route for early settlers crossing from the island of Dongiyu to Luerhmen, southwestern Taiwan. Water masses in the marine areas of TJNP comprise the two main sources, that is, the South China Sea water in summer and the Kuroshio Branch water in the other seasons (Jan et al. 2002). Most previous studies investigating the fish biodiversity of TJNP primarily focused on the species catalogues and distributions in the nearshore water and terrestrial domains of TJNP (e.g., Kuo and Shao 1999;Kuo et al. 2001;Chen et al. 2013). Currently, however, studies on the environmental settings and spatial species diversity of demersal fish assemblages in the marine areas of TJNP, as well as the links between demersal fish assemblages and environmental settings, are limited (Chen et al. 2020). Moreover, an understanding of the spatial species diversity of demersal fish assemblages and environmental settings based on recent survey data is crucially important as baseline information against which future developments in TJNP can be monitored and assessed.
Therefore, the main goals of this study were to obtain the first baseline information on the spatial species diversity of demersal fish assemblages and environmental settings in the nearshore and offshore areas of TJNP and to link environmental (abiotic/biotic) variables of interest to species diversity. To achieve these goals, we tested whether demersal fish species exhibit different patterns of habitat use that are reflected in community structure in terms of spatial species diversity. Here, we addressed three questions: (1) What are the habitat characteristics of the nearshore (shallower) and offshore (deeper) areas in this subtropical park? (2) Do the nearshore and offshore demersal fish assemblages demonstrate significant variability in species compositions in terms of beta diversity? (3) What are the different habitat types and their dominant demersal species and environmental predictors representing the associated fish assemblages? To answer these questions, we used practical permutation methods and multivariate ordination approaches that have been frequently used to elucidate species compositions among different areas in marine ecosystems (Piet and Rijnsdorp 1998;Williams et al. 2010;Amezcua et al. 2019;Yildiz et al. 2019) and to reveal the relationships between demersal fish assemblages and environmental variables (De Leo et al. 2012;Schultz et al. 2014;Galaiduk et al. 2017;Wellington et al. 2018;Amezcua et al. 2019).

Study area
Seven sites located in the nearshore and offshore areas of TJNP (Taiwan) were used for sampling demersal fishes and for collecting environmental data. Specifically, three sites (C1-C3) were nearshore (shallow) sites, and four sites (T1-T4) were offshore (deeper) sites (Fig. 1). Nearshore sites C1 and C2 were located off the coast of Chiku (or Chigu/Qigu) Lagoon, and site C3 was located north of the Tsengwen River estuary (Fig. 1). Offshore sites T1-T4 were situated at the cross-section of the Penghu Channel, the gate of the Taiwan Strait (Fig. 1). All the sampling sites were located more than ca. 3 km from the coastline. Note that none of the sampling sites was chosen in the southern portion of the nearshore area because artificial reef barriers have been placed on the seafloor to prevent bottom trawling (Fig. 1).

Collection of environmental data
Data on five hydrographic variables (i.e., oceanographic conditions including depth, temperature, salinity, density, and dissolved oxygen) of the water column were collected at sampling sites using a shipboard rosette system with conductivity-temperature-depth (CTD) and auxiliary sensors. Near-bottom water samples were collected at the trawling sites with the rosette water sampler system. Concentrations of nutrients (including nitrate, nitrite, 1 3 phosphate, silicate, and ammonia), suspended particles, and chlorophyll-a in the near-bottom water samples were analyzed according to the methods of Meng et al. (2008). Additionally, sampling of surface sediments (substrates) at each site was conducted using a Smith-McIntyre grab sampler (Cat. No.: 5144-BH; Rigosha & Co., Ltd.) before beam trawling. The compositional distributions of substrate grain size were obtained using a Beckman-Counter particle counter (model: LS-100 or LS-13 320). For each trawl haul, the median of the substrate grain size distribution was calculated based on the average of 3 replicates. The classification criteria for substrate grain sizes smaller than 1 mm are referred to in Wentworth (1922). Furthermore, the percentage of organic matter in substrate materials for each haul was conventionally estimated by burning out the organic matter in 3 dry substrate subsamples. Two trawl hauls without data on substrate grain size and organic content were compared with the most recent observations at the same sites in this study. Note that 2 sites (T2 and T4) were surveyed with gravel materials by beam trawling. Therefore, we treated substrate type as a categorical environmental variable for each trawl haul: 1-mixed sand and gravel, 2-sandy (median grain size: 62.5-1000 μm), and 3-muddy (median grain size: < 62.5 μm). Overall, the environmental (including abiotic and biotic) characterization for each trawl haul included the following variables: trawling depth; temperature; salinity; density; dissolved oxygen; Shannon index (fish assemblage variable); substrate type (3 binary categories), substrate organic content and nitrate, nitrite, phosphate, silicate, ammonia, suspended particles, and chlorophyll-a concentrations.

Demersal fish sampling and species identification
Fish sampling (44 trawl hauls) was conducted using a 6-m experimental beam trawl (body: 20-mm mesh; codend: 15-mm mesh) onboard RV Ocean Researcher III (Taiwan) at the seven sampling sites from April 2016 to May 2019 ( Fig. 1). For each trawl haul, the beam trawl was deployed to touch the seafloor at the sampling site, with a towing speed of ca. 2 knots for a period of 30 min. The fish specimens caught by each trawl haul were immediately extracted from the catches, roughly classified by species and frozen at − 20 °C on board. Subsequently, all fish specimens were transferred to the laboratory on land for taxonomic identification and morphological measurement (e.g., fish length and body mass).
The identification of fish specimens was conducted to the lowest level of taxonomic hierarchy possible based on the following references: Shen (1993), Nakabo (2013), Froese and Pauly (2021), and Shao (2021). Fish species nomenclature was mainly accomplished with reference to Froese and Pauly (2021) and Shao (2021), unless recently published articles on the collected fish species could be obtained. Moreover, the total length (T L , mm) and body mass (M B , g) of each fish specimen were measured to the nearest 1 mm and 0.1 g, respectively. Two types of species composition data (presence/absence data and abundance) were used for haul-specific species data in statistical analysis. The species abundance data of each haul were standardized as "abundance/swept area" according to the method of Pacunski et al. (2016) to represent the number of individuals in square kilometers for each species (i.e., ind. km −2 ).

Measurement of alpha diversity of demersal fish assemblages
Using Primer, alpha-diversity indices of the nearshore and offshore assemblages were computed based on haul-specific species abundance data. More specifically, we used species richness [d = (S − 1)/ln(N), where S = total species in each trawl haul and N = total individuals in each haul; Margalef 1951], the Shannon index [H′ = -SUM(P i × ln(P i )), where P i = the proportion of individuals for species i in a haul sample; Lloyd et al. 1968], evenness ([J′ = H ′ /ln(S); Pielou 1966] and the Simpson index [1 − λ′ = 1 − SUM(N i × (N i − 1)/ (N × (N − 1)), where N i = the number of individuals for species i in a haul sample; Simpson 1949] for alpha-diversity indices. We used box plots and the Mann-Whitney U test to compare these indices between the nearshore and offshore assemblages.

Multivariate analyses
Principal coordinate analysis (PCO or PCoA), an unconstrained ordination method, was used to visualize the patterns of environmental data between areas (as well as among the seven sites). We conducted PCO based on haul-specific data on environmental variables using Euclidean distance. Note that the concentrations of nitrite, phosphate, ammonia, and chlorophyll-a, which were below the minimum detection level (MDL), were denoted as the MDL and used in the multivariate analyses. Environmental (abiotic and biotic) data on trawling depth (or bottom depth/water depth) and the concentrations of nitrate, nitrite, phosphate, silicate, suspended particles, and chlorophyll-a were log-transformed (i.e., log(X)) before normalization to construct the resemblance matrix because these data were highly right-skewed. Additionally, we conducted cube transformation for salinity data because the data were left-skewed.
Another unconstrained ordination method, nonmetric multidimensional scaling (NMDS), was used to visualize the patterns of species composition between areas (as well as among the seven sites) based on the resemblance matrix constructed by the Jaccard similarity measure of the haul-specific presence/absence species data. Subsequently, a permutation test was conducted using the constructed Jaccard similarity matrix to test for homogeneity of the compositional species variability between the nearshore assemblage and the offshore assemblage in terms of beta diversity (Anderson et al. 2006).
To link species abundance data to environmental variables, we used the 40 most important species from the 44 trawl hauls as inputs of the species abundance data in multivariate analyses to reduce the influence of rare species on statistical results. These 40 most important species were selected from the haul-specific species abundance data using PRIMER v7 (PRIMER-E: Plymouth, UK; Anderson et al. 2008). Furthermore, a distance-based linear model (DistLM) was used to examine the significance of each environmental variable contributing to species abundances (Legendre and Anderson 1999;McArdle and Anderson 2001). We further used model selection procedures based on both the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) to identify a parsimonious model representing the relationships between species abundances and environmental variables. Consequently, according to the parsimonious model selected, distance-based redundancy analysis (dbRDA), a constrained ordination method, was conducted to gain deep insight into the relationships between species abundances and the contributed environmental variables visually. For plots of dbRDA (or PCO) results, the overlay vectors of environmental variables and fish species were produced based on Pearson correlations greater than 0.4 to represent simple linear correlations of individual variables with the dbRDA (or PCO) axes. We performed all the diversity index calculations, multivariate analyses, and permutation tests (9999 permutations) presented above using PRIMER v7. Guidelines for using this statistical software can be found in Anderson et al. (2008) and Clarke and Gorley (2015).
For further insight, the PERMDISP routine was conducted with a permutation test to examine the homogeneity of species compositions between areas (i.e., beta diversity) based on the resemblance matrix constructed by the Jaccard similarity measure using "area" as a grouping factor. The Bray-Curtis similarity of the log-transformed haulspecific species abundances (i.e., log(X + 1)) was used for the dbRDA. To link species abundances to environmental variables, we first fit a DistLM based on the R 2 selection criterion using the environmental and abundance data of the 40 most important species. Subsequently, we consulted the marginal test for each environmental variable shown in the DistLM results. We then discovered significant environmental variables (P < 0.05). To save time when running model selection procedures, we used the significant environmental variables to run the "best" selection procedure based on either the AIC or the BIC to determine the best combinations of environmental variables. We extracted the 20 best results obtained from the AIC and BIC procedures and created a biplot of these results to identify a parsimonious model. After choosing a parsimonious model from the biplot, we fit a DistLM again using the environmental variables selected by the parsimonious model and obtained dbRDA plots with overlay vectors of both environmental variables and fish abundances. Ultimately, fish habitat types in the study areas were identified according to the dbRDA results, and permutational multivariate analysis of variance (PERMANOVA) tests were conducted to compare the fish abundance data between habitat types. Similarity percentage analysis (SIMPER analysis) was then used to identify the percentage contribution of each species to the Bray-Curtis dissimilarity metric between pairs of habitat types. The relative abundances (RAs) of each fish species by percentage were calculated for all the habitat types identified, and the top 10 most dominant species were determined according to these RAs.

Habitat characteristics
Nearshore sites C1 and C2 were the most intensively trawled sites (i.e., 13 hauls at each site) because the seabed topography around these sites is shallow and flat with soft-bottom substrates (Table S1; Fig. 1), making them more suitable for bottom trawling than the other sites. Four trawl hauls at site C3 were conducted on the inclined seabed topography off the Tsengwen River estuary ( Fig. 1). At the deeper offshore sites (T1-T4), a total of 14 trawl hauls were conducted (Table S1; Fig. 1). Overall, most of the site trawls revealed soft-bottom habitats with sandy and muddy substrates, except for 2 hauls at T2 and 4 hauls at T4, which revealed a substrate type of mixed sand and gravel (Table S1).
As shown in Table 1, nearshore sites C1-C3 had relatively shallow trawling depths (i.e., 12-76 m) compared to those (93-167 m) of offshore sites T1-T4. On average, site T4 had the largest median substrate grain size, while site T1 had the smallest median substrate grain size. Organic content was found to be higher at sites C1, C2, T1, and T4 than at sites C3, T2, and T3. According to hydrographic data from the near-bottom water of the sampling sites, it was clear that the deep-water sites (i.e., T1-T4) had relatively low water temperature, dissolved oxygen levels, and chlorophyll-a concentrations but relatively high salinity levels, densities, and nitrate and silicate concentrations. The concentrations of nitrite, phosphate, and ammonia tended to be low at all sampling sites. The concentrations of suspended particles exhibited no clear spatial trends. The chlorophyll-a concentrations were higher at sites with shallower water than at those with deeper water.
The unconstrained ordination by PCO of environmental variables (Fig. 2a) demonstrated some probable relationships among them. For instance, the hydrographic settings of nearbottom water at the nearshore sites were characterized by shallower water depths, higher temperatures and dissolved oxygen and lower salinity levels, densities, and nitrate, phosphate, and silicate concentrations than those at the offshore sites. The concentrations of nitrate, phosphate, and silicate in near-bottom water tended to be positively correlated with water depth. Conspicuously, the temperature and dissolved oxygen of near-bottom water were negatively correlated with water depth. Additionally, the substrate type of the nearshore sites was dominated by sandy substrates. In contrast, the substrate types differed more conspicuously among offshore sites; that is, there were 2 hauls (at T1) with muddy substrates, 6 hauls (2 at T2 and 4 at T4) with mixed sand and gravel substrates, and 6 hauls (T2, T3, and T4) with sandy substrates (Table S1).

Species composition and spatial diversity of demersal fish assemblages
As listed in Table S2, 30 trawl hauls were conducted at the nearshore sites (C1-C3), and 726 fish individuals belonging to 28 families, 45 genera, and 61 species were collected. At the offshore sites (T1-T4), 427 fish individuals belonging to 32 families, 54 genera, and 67 species were obtained from 14 hauls. In total, fish specimens belonging to 47 families, 84 genera, and 113 species were identified from the nearshore and offshore assemblages. Moreover, an additional 15 fish individuals that could not be identified to the species level in the 44 haul samples was excluded from the statistical analysis because these fishes contributed little to influencing the statistical results of this study. For the nearshore assemblage, the dominant families (see Table S2) in terms of species number (in parentheses) included Bothidae (8), Cynoglossidae (8), Playcephalidae (4), Sciaenidae (4), Apogonidae (3), Leiognathidae (3), Paralichthyidae (3), and Soleidae (3). For the offshore assemblage, the dominant families comprised Bothidae (8), Scorpaenidae (5), Playcephalidae (4), Sciaenidae (4), Paralichthyidae (4), Soleidae (4), and Serranidae (3). The 40 most important species (Table 2) used in the multivariate analyses accounted for 81.5% of the total number of fish individuals. NMDS ordination (Fig. 2b) revealed the patterns of dissimilarity in the assemblage compositions between the nearshore area and the offshore area (as well as among the seven sites). The demersal fishes at nearshore sites C1 and C2 formed a group separate from those at the other sites (Fig. 2b). Additionally, the demersal fishes from site C3 served as an intermediate group between the group of nearshore sites C1 and C2 and the group of offshore sites T1-T4 (Fig. 2b).
The four classic alpha-diversity indices revealed no apparent differences between the nearshore assemblage and the offshore assemblage (Mann-Whitney U test: P > 0.05; Fig. 3). Beta diversity, as demonstrated by a permutation test, however, showed significant differences in species composition between the nearshore and offshore assemblages (P < 0.05).

Environmental influences on species distribution and community structure
The DistLM marginal test (Table 3) indicated that the environmental (abiotic and biotic) variables, including depth, substrate type, organic content, temperature, salinity, density, dissolved oxygen, nitrate concentration, silicate concentration, and the Shannon index, were significant predictors of fish abundances. As shown in Fig. 4, a parsimonious model with only four predictor variables, i.e., depth, substrate type (1), nitrate concentration, and the Shannon index, was identified, and this information was used to create the dbRDA plots in Fig. 5. Additionally, these four variables were all significant predictors, as shown in the sequential test (Table 4). The first two axes of the dbRDA explained ca. 30% of the variability in the fitted model (Fig. 5). A depth gradient was revealed across the nearshore and offshore sites (Fig. 5c), and a relatively high Shannon index was noticeably associated with site T4 (Fig. 5d), which was also associated with substrate (1), i.e., mixed sandy and gravel substrate. In Fig. 5a, the overlay vector (using Pearson correlation) of nitrate concentration was shorter than Additional details of environmental variables on a are given in Tables 1 and S1. Note that one outlier of trawl data from site T2 on the NMDS plot was removed to identify the area-based pattern the other vectors. According to the dbRDA results (Fig. 5), three habitat types (I: shallow soft bottom; II: deeper soft bottom; III: deeper bottom with mixed sand and gravel substrates) were identified in the study areas (permutational t test: I vs II, P < 0.001; I vs III, P < 0.05; II vs III, P < 0.05). The fish species contributing prominently to the nearshore assemblage (habitat type I; see Fig. 5 and Table S3) were the large-tooth flounder Tarphops oligolepis (#27), the snakefish Trachinocephalus myops (#3), and the carpet sole Liachirus melanospilos (#34). For the offshore assemblage, the sharptooth seabass Synagrops philippinensis (#14), the ovate sole Solea ovata (#35), and the shortheaded tonguesole  Cynoglossus kopsii (#39) contributed conspicuously to the group of offshore sites characterized by a soft bottom (habitat type II; Fig. 5 and Table S3). In contrast, the sepia stingray Urolophus aurantiacus (#1), the scorpionfish Scorpaena miostoma (#7), the stonefish Scorpaenodes crossotus (#8), Whitehead's basslet Plectranthias whiteheadi (#15), the grub fish Parapercis sexfasciata (#25), and the lefteye flounder Parabothus taiwanensis (#32) predominated in the group of offshore sites characterized by a substrate type of mixed sand and gravel (habitat type III; Fig. 5 and Table S3). The top 10 dominant species in terms of RA for the three habitat types identified in this study are presented in Table 5. In habitat type I, the large-tooth flounder T. oligolepis was the most dominant species, followed by the snakefish T. myops and the carpet sole L. melanospilos. In habitat type II, the croaker Johnius distinctus predominated in terms of abundance, followed by the shortheaded tonguesole C. kopsii and the Formosa grenadier Coelorinchus formosanus (Table 5). In habitat type III, the scorpionfish S. miostoma predominated in terms of abundance, followed by the sepia stingray U. aurantiacus and the grub fish P. sexfasciata (Table 5). Additionally, several deep-sea benthopelagic and benthic fishes were collected in habitat types II and III (Table S2), e.g., Uroconger lepturus (slender conger; T1), Hime japonica (Japanese thread-sail fish; T4), C. formosanus (Formosa grenadier; T2), Coelorinchus multispinulosus (spearnose grenadier; T2), Neomerinthe procurva (curvedspine scorpionfish; T4), Neomerinthe rotunda (round scorpionfish; T4), Acropoma japonicum (glowbelly; T4), S. philippinensis (sharptooth seabass; T1 and T3), Champsodon snyderi (Snyder's gaper; T2-T4), Acanthaphritis barbata (barbel duckbill; T4), Osopsaron formosensis (Formosa duckbill; T4), and Japonolaeops dentatus (lefteye flounder; T2).

Discussion
The present study reports the first baseline information on the habitat characteristics of demersal fish assemblages in the nearshore and offshore areas of TJNP. The species compositions of demersal fishes differed between the nearshore assemblage (habitat type I) and the offshore assemblage (habitat types II and II), as demonstrated in a multivariate manner ( Fig. 5; Table S3). In total, fishes belonging to 47 families, 84 genera, and 113 species were identified in the nearshore and offshore assemblages. Among the fishes collected, no invasive species are present. As a result of habitat heterogeneity caused by different bottom depths and substrate types, the offshore assemblages (habitat types II and III) possessed more reef-associated and deep-sea fish species than did the nearshore assemblage (habitat type I). In the nearshore assemblage, the first three dominant species in terms of RA were T. oligolepis, T. myops, and L. melanospilos, while S. miostoma, U. aurantiacus, C. kopsii, and J. distinctus predominated in the offshore assemblages, with conspicuous substrate-dependent disparity in species distribution patterns. Among the environmental variables used as predictors of species abundance in multivariate analyses, bottom depth, substrate type, nitrate concentration, and the Shannon index were the most significant predictors. As shown by the DistLM marginal test (Table 3), various hydrographic variables of interest (e.g., temperature, salinity, density, dissolved oxygen, and nutrient concentrations) were also significant environmental predictors of species abundances. Most of these significant hydrographic variables, however, demonstrated strong depth-dependent properties (multicollinearity-strong intercorrelations and an amount of redundancy among these variables) and contributed little to the explained total variations in the dbRDA, as "depth" was incorporated as a predictor. The westernmost site, T4, which is located close to Dongiyu Island (Fig. 1), possessed a mixed sand and

Fig. 4 A scatterplot of the Akaike information criterion (AIC) and
Bayesian information criterion (BIC) values based on the model selection procedure of either AIC or BIC to obtain the top 20 best results from the distance-based linear models (DistLM) of modelling demersal fish abundances in relation to environmental variables in this study. Nine of these best models overlap and are indicated in curly brackets with numbers. Numbers in curly brackets by sequence indicate the environmental variables selected in each parsimonious model (1-bottom depth; 2-mixed sand and gravel substrate; 3-sandy substrate; 4-muddy substrate; 5-nitrate concentration in near-bottom water; 6-Shannon index). The parsimonious model, identified by the combination of four variables (1, 2, 5, 6), was used to demonstrate the dbRDA results in Fig. 5 1 3 gravel substrate and had higher Shannon indices than the other offshore sites (Fig. 5). Note that the sea area around Dongiyu Island has one of the healthiest coral reef ecosystems of Taiwan (Chen et al. 2019; https:// www. marine. gov. tw/). At site T4, four scorpionfish species and three grouper species were collected, revealing that more reef-associated fish species inhabited this site than the other offshore sites.

Constraints and limitations
The beam trawling method used in the present study provides a series of snapshots of life sampled over a stretch of approximately one nautical mile of the bottom. As a result, we investigated the diversity of demersal fish assemblages among sites rather than among time periods. Beam trawling can be used to target a wide range of bottom-living species,  Table 2; c Bubble plots of haul-specific depths of trawling; d Bubble plots of haul-specific H′ values. Note that a zero value of H′ indicates only one fish species caught from a trawl haul including fishes and various invertebrate species (Philippart 1998), e.g., shrimps, crabs, and cephalopods. Nonetheless, small eel-like fishes and large benthopelagic fishes that can swim at high speeds, such as Epinephelus awoara (yellow grouper) and Parupeneus spilurus (blackspot goatfish), were caught infrequently in the present study. Although bottom trawling targets demersal fishes, several pelagic fishes can be caught in nets during lowering and retrieving; for example, Commerson's anchovy, Stolephorus commersonnii, was caught at site C3. Moreover, S. commersonnii can inhabit water depths up to 50 m (Froese and Pauly 2021). We did not remove S. commersonnii from the database (#2 in Table 2) because omission of this species had no appreciable effects on the results of multivariate analyses; as a result, this species was not shown in Figs. 5b and Table 5. Because nearshore sites shallower than 10 m in depth and the southern portion of the nearshore area were not surveyed in the present study, the inferences here are necessarily limited to the fishes obtained from the surveyed habitats across the nearshore area. Latitudinal and longitudinal changes in the community structure of demersal fish assemblages were not analyzed because these geographic factors can confound with the factors of depth gradient and substrate type. Our study sites are not located in brackish estuaries, which can be largely influenced by tides, waves, freshwater flow, and anthropogenic forcing and classified into different estuary types (e.g., Pease 1999; Roy et al. 2001;Saintilan 2004;Amezcua et al. 2019). As a result, anthropogenic-, seasonal-and year-dependent variabilities in demersal fish assemblages were not analyzed in our study because of limitations on data collection.

Multivariate analyses of environmental influences on demersal fish assemblages
We conducted PCO of haul-specific environmental data based on Euclidean distance, and the distribution pattern of these data was equivalent to the pattern shown by principal component analysis (Anderson et al. 2008). The first three axes of the PCO explained an acceptable portion of the total variance (> 59%) for environmental data. Additionally, the dbRDA plot (Fig. 5b) demonstrated a pattern similar to that in the PCO plot of the environmental variables (Fig. 2a). This result suggests that the four-variable dbRDA model indeed represents the overall patterns of the major prominent variability for the environmental variables (Anderson et al. 2008). Table 4 Sequential test results from the distance-based linear modelling (DistLM) analysis of demersal fish abundances in Taijiang National Park, Taiwan, using trawling depth (D T ), substrate type (1), nitrate concentration in near-bottom water and Shannon index (H′) as predictor variables (see Tables 1 and S1

3
In the present study, the main faunal change in demersal fish assemblages was attributed to the effects of bottom (trawling) depth and substrate type ( Fig. 5; Table 1). Changes in the community structure of demersal fish assemblages caused by different bottom depths and/or substrate types in marine ecosystems have been well documented in the literature (e.g., Rainer and Munro 1982;Bianchi 1991;Fujita et al. 1995;Connell and Lincoln-Smith 1999;Hyndes et al. 1999;Demestre et al. 2000;Araújo et al. 2002;Beentjes et al. 2002;Prista et al. 2003;Damalas et al. 2010;García-Rodríguez et al. 2011;Reum and Essington 2011;Schultz et al. 2014;Zintzen et al. 2017;Stephenson et al. 2020). Nevertheless, the present study is the first to report the species compositions of demersal fish assemblages in relation to various hydrographic variables and substrate types in TJNP, a subtropical national park in Taiwan. Here, we found that bottom depth explained 19.6% of the total variance in the DistLM (Table 4) of environmental influences on demersal fish assemblages in this subtropical national park. Substrate type 1 (mixed sand and gravel substrates) and the fish assemblage variable "Shannon index" were the other important predictors, explaining ca. 9% of the total variance, indicating that patterns of substrate preference differed among fish species and thus affected spatial species diversity among sites. Although the Shannon index is sensitive to sample size (Magurran 2004), this index can be easily calculated and has frequently been used to measure the structure complexity of fish assemblages (e.g., Sánchez-Hernández et al. 2018). Here, we simply use this index to look at how this fish assemblage variable co-vary with environmental variables as well as trawling sites. Higher Shannon indices tended to be related to the trawl hauls at site T4 and several hauls at site C2 ( Fig. 5; Table 1), indicating that demersal fishes at these two sites seem to be more diverse in terms of structure complexity than the other sites.
Depth-dependent hydrographic variables, including temperature, salinity, density, dissolved oxygen, and nitrate and silicate concentrations, contributed only a minor portion of the explained variation when bottom (trawling) depth was incorporated as a predictor. The relatively low concentrations of dissolved oxygen in near-bottom water at the deeper offshore sites (Table 1) are considered to pose little risk to demersal fish survival. Although the parsimonious model incorporated the depth-dependent variable of nitrate concentration (Figs. 4 and 5a), this variable explained only ca. 3% of the total variance. Why was nitrate concentration not excluded from the parsimonious model when bottom depth was incorporated as a predictor? One possible explanation is that some sporadic high-nitrate-concentration events occurred in the near-bottom water at shallow site C2 (Table 1), resulting in additional variability that could not be explained by bottom depth in the multivariate analyses. The seabed topology of site C2 is shallow and flat, but this site is close to an area with an inclined seabed (Fig. 1). Sporadic high-nitrate-concentration events in near-bottom water at site C2 could be attributed to various sources caused by biogeochemical cycling, e.g., nutrient entrainment from deeper water (coastal upwelling) and terrestrial contaminants. Nonetheless, verifying the sources of these events is beyond the scope of our study.
Trachinocephalus myops was used as the nomenclature for the snakefish in the present study, while Polanco et al. (2016) suggested that T. trachinus be used for the populations of the T. myops species complex from the Indo-West Pacific Ocean. The demersal fishes in our dataset included two recently described new species, i.e., the sandperch Parapercis moki (Ho and Johnson 2013) and the lizardfish Synodus pacificus (Ho et al. 2016). Both species were collected at site T4 (sampling depths of 105-116 m; see Table S2), suggesting that site T4 and its surrounding water possess comparatively great potential for the discovery of new fish species as a result of higher demersal fish biodiversity than observed at the other sites.

Life-history traits and habitat use of dominant demersal fish species
The fish specimens collected in the nearshore area were mostly small-and medium-sized fish (Table S2). The largetooth flounder T. oligolep was the most dominant species in terms of abundance in the nearshore assemblage (habitat type I; Table 5), but it was also caught with lower abundances at sites C3, T3, and T4, with water depths reaching up to 158 m (Table S2). Several imminent-spawning (hydrated oocyte-possessing) adults of T. oligolepis were present at site C2 in April 2018, suggesting that April is one of the months in which this species spawns in the study area. The snakefish T. myops was the second most dominant species in the nearshore assemblage. Juveniles and subadults of this species were frequently caught at sites C1 and C2 in summer and autumn, but they were not caught at sites C3 and T1-T4 (Table S2). Additionally, only a few adults of these species were collected during the study period (Table S2). These results suggest that T. myops exhibits seasonal distribution patterns in habitat use in the coastal waters of western Taiwan. Moreover, both T. oligolep and T. myops have rarely been found to use Chiku Lagoon and the Tsengwen River estuary (Fig. 1) as habitats (Kuo and Shao 1999;Kuo et al. 2001), indicating that they are exclusively marine species. Additionally, these two species might be useful as indicators for monitoring changes in the community structure of demersal fish assemblages and the associated environmental conditions in TJNP due to their dominance in the nearshore assemblage.
In the offshore demersal fish assemblage (habitat types II and III), more reef-associated and deep-sea fishes were caught because of the more diverse habitats (i.e., depth gradient and substrate types) compared to those in the nearshore assemblage (habitat type I). The sharptooth seabass S. philippinensis, ovate sole S. ovata, and shortheaded tonguesole C. kopsii were the major species in the trawl haul group characterized by a habitat with a deeper soft bottom, i.e., habitat type I (Fig. 5b). Among these three species, S. philippinensis was collected with low abundances at sites T1 and T3, while S. ovata and C. kopsii were collected with relatively high abundances (Table S2). Moreover, C. kopsii was the most widely distributed species among the study sites, but this species was not collected at site T4 (Table S2). Juveniles, subadults, and adults of C. kopsii were all found at the study sites, indicating that it is a resident species. The reef-associated scorpionfishes S. miostoma and S. crossotus predominated at sites T2 and T4 (Table S2), which were characterized by mixed sand and gravel substrates rather than sandy and muddy substrates (Table S1). Both S. miostoma and S. crossotus were rarely found at the other sites, indicating the strong substrate dependence of habitat use by these species (Table S2). Scorpaena miostoma was caught at various life-history stages at sites T2 and T4 (Table S2), suggesting that it is also a resident species. The serranids (groupers), namely, E. awoara, P. whiteheadi, and S. analis, were the other reef-associated species collected, but only P. whiteheadi dominated at both T2 and T4 (Table S2). Similar to S. miostoma, P. whiteheadi is suggested to be a resident species because its individuals have been recorded at various life history stages (Table S2). Similar to S. miostoma and P. whiteheadi, the sepia stingray (oriental stingaree) U. aurantiacus was exclusively found at sites T2 and T4, indicating comparable substrate preferences among these species. Urolophus aurantiacus is distributed in the adjacent waters of the East China Sea and prefers inhabiting rocky habitats (Last and Marshall 2006). This species has been evaluated to be "Near Threatened" on the IUCN Red List of Threatened Species (Last and Marshall 2006). In our study, this species was a dominant species in the western part of the offshore area (Table S2), and one pregnant individual was collected at site T4 in April 2018. In summary, the mixed sand and gravel substrates at sites T2 and T4 created a rough seabed with gravel materials providing suitable microhabitats (e.g., diverse patches of shelter habitat) that facilitated colonization by reef-associated fishes and other marine organisms.

Implications for biodiversity conservation in the existing spatial management areas
Understanding the spatial distribution patterns of biodiversity is fundamental to the development of management and conservation strategies at a regional scale, but this information is lacking for the marine areas of TJNP. In previous studies, data on fish occurrence were routinely collected to develop species catalogues, and these data might be useful for other purposes, e.g., assessing spatial community structure and species-habitat linkage (Kuo and Shao 1999;Kuo et al. 2001;Chen et al. 2020). To manage large areas of the seafloor, species distribution and habitat suitability data obtained by mapping are likely to be collected in the future (Hewitt et al. 2015;Costello and Chaudhary 2017). We suggest that creating indices of coastal habitat health in the marine ecosystems of TJNP is required for monitoring and assessing the environmental conditions suitable for demersal fish survival.
In recent decades, fishing and other anthropogenic impacts have led to declines in fish stocks and habitat destruction in the coastal environment of Taiwan (Shao 2009;Chen 2010;Kuo and Booth 2011). Therefore, efforts to restore marine fish resources have focused on spatially explicit management methods to protect fishes and the habitats they require for survival (Shao 2009), e.g., the design and implementation of marine protected areas (MPAs). The inability to discern the relationships between fish and their environment can clearly obstruct conservation and management measures for fish resources (Elliott et al. 2016). To meet conservation and management goals, there have been increasing investigations aimed at mapping fish habitats vulnerable to anthropogenic forcing and regional climate change and at identifying fish resource needs (Hsieh et al. 2007;Shao 2009). The nearshore area (habitat type I) should be reserved exclusively for artisanal and small-scale fisheries using nontrawling types of gear. In the offshore 1 3 area (including habitat types II and III), westernmost site T4 demonstrated a bottom type of mixed sand and gravel (creating microhabitats for marine organisms) and relatively high Shannon indices, suggesting that the water adjacent to this site represents an area of relatively high biodiversity in TJNP. Currently, the offshore area is not an MPA or a sanctuary; rather, it is termed an "existing use" area, which has received little management or control related to marine resources. To enhance MPA networks and to conserve biodiversity of demersal resources, the western part of the offshore area might be expanded in scope to cover more microhabitats, serving as an MPA or a sanctuary to diminish fishing pressure from commercial bottom trawling fisheries. Moreover, efforts to evaluate possible approaches and research requirements across biodiversity organizations can facilitate the planning of future studies and promote more effective communication among stakeholders, scientists and policy makers. Establishing a general approach that captures how various field studies have focused on different biodiversity components can also contribute to meta-analyses of worldwide experience in scientific research to support ecosystem-based management (Ellis et al. 2011). Ultimately, the establishment of a marine environmental monitoring system should be emphasized as one of the major goals of current conservation and management plans. In terms of implementing robust monitoring programs for TJNP, spatial stratification of survey data might thus be necessary to represent possible confounding effects of fishing-induced and environmentally induced heterogeneity in the indicator. Such stratification should be applied to evaluate environmental factors that could drive changes in the spatial species diversity and community structure of demersal fish assemblages. For example, climate change will likely shift the defined fish assemblages northward, resulting in the immigration of new species and thus altering species distribution patterns and community structure (Shephard et al. 2011). Apart from implications for monitoring methodology, research efforts and the productivity and spatial pattern of demersal fish assemblages under a changing climate should be redefined (Free et al. 2019). Environmental education regarding the sustainable development of marine ecosystems should be routinely enhanced by TJNP headquarters to increase public awareness of biodiversity conservation issues.
In conclusion, the present study provides the first baseline information of the environmental conditions in the marine ecosystems of TJNP, a subtropical national park in Taiwan. The spatial species diversity of the nearshore and offshore demersal fish assemblages in TJNP is demonstrated across multiple depths and environmental covariates using practical multivariate analyses. The results offer insight into the general ecological features of demersal fish assemblages in this subtropical national park. Specifically, three habitat types (I: shallow soft bottom; II: deeper soft bottom; III: deeper bottom with mixed sand and gravel substrates) in the study areas are identified primarily according to bottom depth (trawling depth) and substrate type. Such habitat heterogeneity increases the spatial species diversity of demersal fish assemblages. Moreover, hydrographic setting differences between the nearshore and offshore assemblages are noticeable in the multivariate environment-and taxonbased analyses. In addition to water depth, several hydrographic variables of interest, including temperature, density, dissolved oxygen, and nitrate and silicate concentrations in near-bottom water, are significant predictors of fish diversity. Nonetheless, most of these significant variables are highly correlated with each other and exhibit strong depth dependence. Ultimately, the identification of species distribution patterns and habitat characteristics here pinpoints priorities for future investigators that will further enhance our understanding of the demersal fish assemblages in this subtropical national park and the drivers that act to shape the assemblages.
Data availability All data generated or analyzed during this study are included in this published article and its supplementary information files.
Author contribution M-H Chen and C-Y Chen conceived and designed research. H-S Chen, K-S Chen and Y-L Su conducted field sampling, species identification, and statistical analysis. P-J Meng dealt with hydrographic data. K-S Chen and M-H Chen wrote the manuscript. All authors read and approved the manuscript.
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/.