Considering species functional and phylogenetic rarity in the conservation of fish biodiversity

Rare species make substantial contributions to coastal ecosystem functions. Functional rarity (FR) and phylogenetic rarity (PR) are important features for biodiversity conservation. This work aimed to discuss the necessity and reasonableness of conserving fish FR and PR in coastal seas.


| INTRODUC TI ON
Coastal habitats are essential for ecological functions and provide innumerable benefits and services for human beings (Barbier et al., 2011).Fish conservation in megadiverse developing countries is imminent as fish biodiversity faces huge risks from overfishing, species invasion, pollution and climate change (Lövei & Lewinsohn, 2012;Vitule et al., 2009).The China Seas, including the Bohai Sea, the Yellow Sea, the East China Sea and the South China Sea, extend over 38° of latitude and span from tropical to temperate climate zones (Wang et al., 2012).With abundant coastal resources, China makes the largest contribution to global fisheries production (FAO, 2020), which outlines the importance of its fisheries.However, the coastal ecosystem in China Seas has experienced severe deterioration in fish community structure during the past half century (Kang et al., 2021).These pressures compromise coastal fish biodiversity and potentially affect the stability of the ecosystem (Beukhof, 2019).A thorough understanding of biogeographic patterns of fish biodiversity is thus necessary so that appropriate conservation measures can be properly schemed and prioritized.
Functional diversity reflecting how an ecosystem operates or functions (Mouchet et al., 2010) and phylogenetic diversity representing the evolutionary history of the community (Faith, 1992) are both important components of biodiversity (Naeem et al., 2012).
They are related to ecosystem functions, ecological services and evolutionary potentials (Winter et al., 2013).Species contribute disproportionately to functional and phylogenetic diversity (Pool et al., 2014).Rare species possibly support vulnerable functions in ecosystems (Mouillot et al., 2013), represent unique evolutionary branches (Winter et al., 2013) and contribute significantly to functional and phylogenetic diversity (Mi et al., 2012).
However, most conservation efforts usually focus on the threatened species on the Red List of the International Union for Conservation of Nature (IUCN) (Rodrigues et al., 2006), which still ignores the functional and phylogenetic attributes of the species.
The concept of rarity is initially derived from the geographic extent, local population size and habitat specificity of rare species (Gaston, 1994;Violle et al., 2017).It has been a central topic of conservation biology aiming to determine the characteristics that possibly incur risks of species extinction and biodiversity shrinkage (Flather & Sieg, 2007).Functional rarity (FR) can be defined as a pool of species that 'integrates both functional distinctiveness and taxonomic scarcity at the local scale' (Violle et al., 2017).The same schemes can be applied to phylogenetic rarity (PR).The difference is that FR is trait-based, whereas PR is derived from certain forms of evolution information such as a phylogenetic tree (Redding et al., 2010).FR and PR are proved useful as surrogates of biodiversity to guide conservation prioritization (Albuquerque & Astudillo-Scalia, 2020).Previous works have revealed biogeographic patterns of FR and PR of marine fish at global scales (Trindade-Santos et al., 2022).However, rarity is scale-dependent (Hartley & Kunin, 2003), and such a coarse resolution such as a global scale of 2 arc-degree (Trindade-Santos et al., 2022) provides limited information for regional conservation practices.
China has been dedicated to marine biodiversity conservation for nearly six decades (Bohorquez et al., 2021).Until 2019, more than 270 operational marine-protected areas (MPAs) have been constructed, comprising over 5% of its national waters (Li, Ren, et al., 2019;UNEP-WCMC & IUCN, 2022).This coverage is insufficient to effectively protect coastal biodiversity (Edgar et al., 2014), and a reasonable expansion of MPAs for effective systematic conservation is desperately needed (Li, Ren, et al., 2019).This expansion will rely on conservation prioritization, which comprehensively takes biological, ecological, economic, political, historical and geographical factors into an integrative framework (Margules & Pressey, 2000).
FR and PR of coastal fish in China Seas have never been studied at finer scales, impeding the development of effective systematic conservation (Qiu et al., 2009).
By compiling historical fish investigation data in the past decades in coastal China Seas, this work focussed on three objectives with the application of joint pecies distribution modelling (JSDM): (1) explore the pools of functionally and phylogenetically rare species in coastal China Seas, (2) discuss the biogeographic patterns of FR and PR and (3) evaluate the efficiency of current MPA networks in terms of match/mismatch with FR and PR.Few studies have discussed coastal fish rarity in China Seas, so we did not presume any patterns of FR and PR, nor matches or mismatches between MPAs and the hotspots of FR and PR.

| Fish occurrence data
Fish occurrence matrix was compiled from 148 fish investigation incidents between 1980 and 2018 (see the fish investigation data reference list in Appendix S4).These investigation incidents occurred independently and were hereinafter referred to as samples.Sampling methods were identified according to FAO definition and classification of fishing gear categories (Nédélec & Prado, 1990).Sampling year was categorized to the period when the fish investigation occurred.There were 14 levels of sampling methods and 62 levels of sampling years (Appendix S2: Figure S1).Original fish taxonomic data were double-checked and validated with FishBase records to avoid invalid species, synonyms and homonyms (Froese & Pauly, 2022).We followed Grenie et al. (2018) and defined the 'threatened' species as those listed as Vulnerable, Endangered and Critically Endangered on the IUCN Red List (IUCN, 2022).There were 384 valid fish species identified altogether, which belong to 108 families (Appendix S5), 12 of which were threatened species.
Historical environmental data for model fitting included multiple sources.Mean water depth (depth) was extracted from the data F I G U R E 1 Research area and sampling sites.Codes, geographic coordinates and other information about the sampling sites are listed in Appendix S1: Table S1.
centre of National Oceanic and Atmospheric Administration (NOAA) (Amante & Eakins, 2009).Types of seafloor surface substrate (sub-  et al., 2018).Depth and substrate remained unchanged among different sampling years because the variations of these two variables were expected to be minor at the same site and no more historical data were available so far.Correlation coefficients between these variables were below .70(Appendix S2: Figure S2a), and values of variance inflation factor were smaller than a threshold of 3 (Appendix S2: Figure S2b), indicating acceptable levels of multicollinearity (Shrestha, 2020).

| Functional traits
We followed Brosse et al. (2021) to perform morphological measurements (Appendix S2: Figure S3).These morphological measurements were recorded by ImageJ Version 1.53k (Schneider et al., 2012).Pictures were taken primarily from regional fish atlases.Only pictures of adults were used and juveniles were not considered because morphological changes can occur during ontogeny.In cases of sexual dimorphism, only male morphology was used, because pictures of females were scarce for most species.
Nine morphological traits (Appendix S1: Table S3) were derived from these morphological measurements, including body elongation, vertical eye position, relative eye size, oral gape position, relative maxillary length, body lateral shape, pectoral fin vertical position, pectoral fin size and caudal peduncle throttling.For those fish with unusual morphologies, we followed Brosse et al. (2021) and Villéger et al. (2010) to apply the following rules: (1) for species with no visible caudal fin (e.g., Sternopygidae, Anguilidae, Plotosidae), caudal peduncle throttling was set to one, assuming that caudal fin depth is equal to caudal peduncle depth; (2) for algae browser species with the mouth positioned under the body (e.g., Loricaridae and some Balitoridae), oral gape position and relative maxillary length were set to zero; (3) for species without pectoral fins (e.g., Synbranchiforms and some Anguiliforms), pectoral fin vertical position was set to zero; and (4) for flatfish, the body depth measure was the body width because the fish lies on one side of its body.
We also followed Trindade-Santos et al. (2022) and extracted other eight biologically and ecologically relevant traits (Appendix S1: Table S3) from the most recent version of the FishBase database (Froese & Pauly, 2022), including maximum lifespan, generation time, food consumption to biomass ratio, maximum body length, position in water column, diet habit, trophic level and body shape.
Considering the inevitable data gaps in a large data set such as FishBase, a random forest algorithm was applied to fill the missing values of the trait matrix (smaller than 10%) (Stekhoven & Buehlmann, 2012).Pearson's correlations between these traits were also tested to ensure complementarity (Appendix S2: Figure S4).These traits were mostly uncorrelated, indicating that they provided complementary information.

| Phylogenetic data
Cytochrome oxidase I (COI) gene sequences were obtained from the GenBank database (Benson et al., 2012).Among the 384 fish species, 355 (92.44%) were found with verified COI gene sequences, but no records were identified for the other 29 species as of 1 July 2021.For a better estimate of the phylogenetic information within the community, surrogate sequences for these 29 species were used by their phylogenetically and morphologically close relatives within the same genera (see the COI accession no. in Appendix S5).
The Akaike information criterion (AIC) was used to select the optimal parameters for the maximum likelihood (ML) phylogeny (Posada & Crandall, 2001).The substitution model 'GTR + I + Γ' was the one that best explains the empirical data in terms of AIC (Appendix S2: Finally, the ML phylogenetic tree with GTR + I + Γ was constructed (Appendix S2: Figure S6) by the executable program 'PhyML 3.1' in R (Guindon et al., 2010).   of China, 2003), the mesh size of most of the sampling incidents was confined to about 20 mm and the sampling periods were restricted to about 1 h for each sampling station at daytime.Therefore, the sampling efforts of these samples were comparable.However, different sampling methods and sampling years might still impose impacts on the patterns of fish communities (Pope & Willis, 1996;Wileman et al., 1996).

| Accounting for the non-independence of the data
According to the data structure (Appendix S2: Figure S7), the non-independence of the samples collected in this work can be attributed to the crossed effects of sampling site, sampling method and sampling year.Therefore, three crossed random effects were considered in subsequent analysis, including random intercept of sampling site to account for spatial non-independence of the samples within sites, random intercept of sampling year to account for temporal non-dependence of the samples, and random intercept of sampling method to account for the variance caused by different sampling methods (Faraway, 2016).

| Joint species distribution modelling
Data were fitted with hierarchical modelling of species communities (HMSC) (Ovaskainen & Abrego, 2020), a JSDM that includes a hierarchical layer asking how species respond to environmental covariates (Ovaskainen et al., 2017).This model accounts for species co-occurrence patterns using a latent variable approach following Ovaskainen et al. (2016).By doing so, species interactions help explain the fish community structure and improve the predictive capability of the model (Ovaskainen & Abrego, 2020;Tikhonov et al., 2017).Species traits and phylogenetic information were included in the model to improve its performance (Tikhonov et al., 2020).Bayesian generalized linear mixed model with a probit link and binomial error distribution was used in HMSC models assuming the default prior distributions (Tikhonov et al., 2020).
Random effects of sampling site, sampling year and sampling method were applied in the models.Posterior distribution was sampled with four Markov Chain Monte Carlo (MCMC) chains, each of which was run for 300,000 iterations.The first 50,000 iterations were removed as burn-in.The chains were thinned by 1000 to yield 250 posterior samples per chain and so 1000 posterior samples in total.
Five competing models were constructed for model validation, including: (1) model_null which only included an intercept and three random effects; (2) model_cov which included an intercept and all predictors, but no random effects; (3) model_full which included an intercept, all predictors and three random effects; (4) model_joint which performed a variable selection jointly for all species based on model_full by a spike and slab prior process in the context of a Bayesian framework (Mitchell & Beauchamp, 1988;O'Hara & Sillanpää, 2009); and (5) model_separate which performed a variable selection individually for each species based on model_full also by the method of spike and slab prior in model_joint, but species selected these variables separately.Markov Chain Monte Carlo convergence was checked by examining the potential scale reduction factors of the model parameters (Gelman & Rubin, 1992).Explanatory and predictive powers of these competing models were evaluated in terms of Area Under Curve (AUC; Pearce & Ferrier, 2000) and Tjur's R 2 (TjurR 2 ; Tjur, 2009).To compute explanatory power, model predictions were projected based on the whole data set.To compute predictive power, a fivefold cross-validation was performed, in which the sampling units were assigned randomly to fivefolds, and predictions for each fold were based on a model fitted to the remaining fourfolds (Fushiki, 2011).Therefore, AUC and TjurR 2 were calculated based on all data sets, indicating the explanatory power of the model, while AUC_CV and TjurR 2 _CV were based on a fivefold cross-validation, indicating predictive power.The model, which presented the highest explanatory and predictive powers, was used for prediction.A threshold occurrence probability of .5 was then used to derive the matrix of species occurrences per pixel (Bailey et al., 2002).Procedures for HMSC modelling are detailed in specifics in Appendix S3.

| Species restrictedness
Species restrictedness refers to the geographical extent of a species.

It is computed as:
with R i the restrictedness of species i, K i the number of sites where species i occurs and K tot the total number of sites in the data set.A value close to 1 indicates the species is restricted in an infinitely small range of distribution and a value close to 0 indicates the species is ubiquitously distributed in every pixel.

| Functional distinctiveness
Functional distinctiveness is the average functional distance from a species to all the others in the community space of a community.It is computed as: with D i the functional distinctiveness of species i, N the total number of species, and d ij the functional distance between species i and species j.
The larger the index value, the more distinct the species is in the functional space.Gower's distance of the trait matrix was used because continuous and categorical variables were both included (Pavoine et al., 2009). (1) (2)

| Phylogenetic distinctiveness
Analogously to functional distinctiveness, phylogenetic distinctiveness quantifies how isolated a species is on a phylogenetic tree (Redding & Mooers, 2006).We applied the method of fairproportion (Isaac et al., 2007) to calculate phylogenetic distinctiveness because it was reported to be simpler and more informative (Cadotte & Jonathan, 2010).It is computed as: where ED(T, i) indicates the phylogenetic distinctiveness of species i on a phylogenetic tree T, e is a branch of length λ e , q(T, i and r) indicates the set of branches on the tree T connecting species i with the root of the tree r and S e is the number of species that descend from the branch e.

| Rare species and rarity
We followed the schemes of Trindade-Santos et al. ( 2022) using the quartile criterion proposed by Gaston (1994) to identify rare species.According to the distribution predictions extracted from HMSC based on the data for all 148 fish investigations, species were defined as functionally rare if they were simultaneously valued within the top 25% quantile range of restrictedness and functional distinctiveness.Likewise, species were defined as phylogenetically rare if they were simultaneously valued within the top 25% quantile range of restrictedness and phylogenetic distinctiveness.FR and PR were then defined as the richness of functionally or phylogenetically rare species in each grid cell, respectively.
To determine whether the rarity differed from random expectation, a restrained null model based on the curveball algorithm was constructed (Strona et al., 2014).The species occurrence matrix in each grid cell was randomized 2000 times while maintaining the species richness of each grid cell and the prevalence of each species.Then, the standardized effect size (SES) of the rarity was calculated as follows: where observed indicates the original value of FR or PR, null is the mean value of the null model and SD(null) is the standard deviation of the values from all repetitions (Swenson, 2014).The grid cells with SES ≥2 were defined as hotspots of FR and PR, respectively (Swenson, 2014;Trindade-Santos et al., 2022).

| Matches between rarity and MPAs
Match and mismatch between the MPAs and the hotspots of rarity were calculated as follows: where intersection indicates the shared grid cells between MPAs and the hotspots of rarity, and total indicates the total grid cells occupied by the hotspots of rarity.Equation ( 5) was also applied to the calculation of the shared hotspots between FR and PR (Trindade-Santos et al., 2022).
A randomization procedure was performed (2000 repetitions) while maintaining the total grid cells of MPAs (Jackson et al., 1992), and Wilcoxon's signed-rank test was used to determine whether the match was greater than random expectation (Bauer, 1972).Wilcoxon effect size (Fritz et al., 2012) was also calculated to determine whether the effect was small (0.10-0.30), moderate (0.30-0.50) or large (0.50-1.00).

| Sensitivity analysis and spatial autocorrelation
Sensitivity analysis was performed to check the robustness of the patterns of rarity (Trindade-Santos et al., 2022).Each functional trait was removed once from the trait matrix with replacement.Functional distinctiveness calculated without this trait was plotted against the values calculated with all the traits, and a quantile regression model was fitted to check whether the estimated distinctiveness was sensitive to the selection of functional traits.Analogously, each tip on the phylogenetic tree was removed once with replacement, and phylogenetic distinctiveness was calculated without this tip to check the sensitivity of phylogenetic distinctiveness to the structure of the phylogenetic tree.If the quantile regression was close to the line y = x, functional distinctiveness or phylogenetic distinctiveness was not sensitive to the selection of functional traits or the construction of the phylogenetic tree.
We constructed a distance decay plot to check whether spatial autocorrelation obscured the patterns of FR and PR (Sara et al., 2022).Pairwise Sorenson's dissimilarities between the composition of rare species in each grid cell and other cells were plotted against the pairwise spatial distances, and a quantile regression model (τ = .5)was fitted.If the quantile regression was close to a horizontal line, the effect of spatial autocorrelation can be considered ineffective in explaining the patterns of FR and PR.

| Model validation for the joint species distribution modelling
All indices were the highest for model_full with an average crossvalidated AUC of 0.82 and TjurR 2 of .49(Table 1), indicating that this model possessed the highest explanatory and predictive power compared with other four competing models.Model_joint and model_sep included a spike and slab prior procedure.AUC of these two models was higher than 0.90, but their AUC_CV and TjurR 2 _CV (3) were the lowest, indicating that they were not suitable for prediction.
Model_cov did not include any random effects, and its explanatory power and predictive power were lower than those of model_full with both fixed and random effects, indicating that fish communities in this study have a certain spatial and temporal structure, which can be explained by random effects.Model_null only included an intercept.Compared with model_cov and model_full, its predictive ability is lower, indicating that explaining variables were essential for modelling the community.Therefore, model_full was the optimum model in this work.

| Sensitivity analysis and spatial autocorrelation
The estimated functional distinctiveness after removing the functional traits one by one was strongly correlated with the values calculated with the full trait matrix (Figure 2a).The slope of the quantile regression model was close to one with an intercept of zero, indicating that the estimated functional distinctiveness was not sensitive to the choice of functional traits.The situation was likewise for PR (Figure 2b), which suggests that the values of phylogenetic distinctiveness were not sensitive to the structure of the phylogenetic tree.
The distance decay plot shows a slope of zero between the pairwise dissimilarities of FR and geographical distance (Figure 2c), indicating that spatial autocorrelation did not obscure the biogeographic patterns of FR.Likewise, no obvious pattern was observed between the pairwise dissimilarities of PR and geographical distance (Figure 2d).

| Rare species
Forty-four species were identified as functionally rare (Appendix S1: Table S4), 32 of which are bony fish (Actinopterygii) and 12 are cartilaginous fish (Elasmobranchii).Only one functionally rare species (Japanese eel Anguilla japonica) was categorized as threatened.Twenty-two species were identified as phylogenetically rare (Appendix S1: Table S5), 20 of which are bony fish and 1 is cartilaginous fish.None of these species was categorized as threatened

| Hotspots of rarity
Hotspots of FR (Figure 3a) were mostly observed in three areas, respectively, around 22.94° N, 31.48°N, and 38.49° N (Figure 4a), including the southern and eastern coast of Taiwan, the Yangtze River Estuary and its adjacent waters, and the Yellow River Estuary and its adjacent waters.These hotspots covered 10.27% of the coastal areas.
Compared with FR, hotspots of PR (Figure 3b) were more sporadically distributed in the coastal East China Sea, the Bohai Sea, and the northern Yellow Sea.These hotspots covered 3.06% of the coastal areas, suggesting fewer geographical distributions than FR.
The latitudinal distribution of PR was also multimodal (Figure 4b), but only two peaks were observed at 30.93° N and 38.87° N.
The shared hotspots between FR and PR occupied 1.33% of the coastal areas (Figure 3c).These areas covered 13.29% of the hotspots of FR and 44.58% of the hotspots of PR.Two peaks representing the most frequent congruence between the hotspots of FR and PR were observed at 31.09° N and 38.67° N (Figure 4c).

| Representation by MPAs
The match between the MPAs and hotspots of FR was 16.16% (Figure 5a).These matched areas were significantly higher than random expectation (p < .001)with a large effect (Wilcoxon effect size = 0.88).The most frequently matched areas occurred at 31.40° N and 38.06° N (Figure 6a), which were compliant with those two distribution peaks of FR (Figure 4a).
The match between the MPAs and hotspots of PR was 20.48% (Figure 5b), which was also significantly higher than random expectation (p < .001)with a large effect (Wilcoxon effect size = 0.87).The most frequently matched areas occurred at 30.65° N and 38.71° N (Figure 6b), which was compliant with the two distribution peaks of PR (Figure 4b).

TA B L E 1
Validation of the competing models for hierarchical modelling of species communities.The competing models include (1) model_null which only included an intercept and three random effects; (2) model_cov which include an intercept and all predictors, but no random effects; (3) model_full which included an intercept, all predictors, and three random effects; (4) model_joint which performed a variable selection jointly for all species based on model_full by a spike and slab prior process in the context of a Bayesian framework; and ( 5) model_separate which performed a variable selection individually for each species based on model_full also by the method of spike and slab prior in model_joint, but species selected these variables separately.b Index of Area Under Curve.
c Index of Area Under Curve with a fivefold cross-validation.
d Index of Tjur's R 2 .
e Index of Tjur's R 2 with a fivefold cross-validation.The match between the MPAs and the shared hotspots between FR and PR was 40.54% (Figure 5c).This coverage was significantly higher than random expectation (p < .001)with a large effect (Wilcoxon effect size = 0.87).However, substantial mismatched areas (59.46%) were still observed in the coastal East China Sea, the Bohai Sea, and the northern Yellow Sea.The most frequently matched areas occurred at 31.33° N and 38.18° N (Figure 6c), which was also compliant with the latitudinal distribution of the shared hotspots (Figure 4c).

| Distinctiveness of rare species
Functional and phylogenetic distinctiveness of rare species can support important ecological roles and contribute disproportionately to ecosystem functions (Dee et al., 2019;Violle et al., 2017).Extinction of these species will possibly cause an abrupt loss of functions and services delivered by the ecosystem (Colares et al., 2022).Phylogenetic distinctiveness may attribute rare species with distinctive phenotypes and rare ecological roles (Safi et al., 2011).Low functional or evolutionary redundancy always occurs within rare species, increasing the vulnerability of functional and phylogenetic diversity (Mouillot et al., 2013).Extinction of these species will possibly cause an abrupt loss of functions and services delivered by the ecosystem (Colares et al., 2022).
Many conservation practices mainly focus on rangerestricted species associated with extinction risk (Hoffmann et al., 2008).However, only considering their threat status will possibly overlook the species with high functional and phylogenetic distinctiveness.For example, phylogenetic distinctiveness in amphibians is slightly negatively correlated with their threat status on the IUCN Red List (Morelli & Møller, 2018), and functional distinctiveness of coral reef fishes at global scales is also independent of their threat status (Grenie et al., 2018).The reliance on the current IUCN Red List as the basis for estimating vulnerability of fishes in megadiverse nations such as Brazil and China might cause bias in coastal fish conservation conservation (Lövei & Lewinsohn, 2012).For instance, many Neotropical fishes are not included among IUCN threat categories simply because of major gaps in spatial coverage for the IUCN Red List (Vitule et al., 2017).This situation highlights the importance of the conservation of species with high functional or phylogenetic distinctiveness.Therefore, we advocate integrated fish conservation by incorporating the threat status, functional distinctiveness and phylogenetic distinctiveness at the species level, especially in megadiverse developing countries.

| Biogeographic patterns of functional rarity and phylogenetic rarity
The peak of FR hotspots in the southern and eastern coast of Taiwan (around 22.94° N) was possibly attributed to the Kuroshio current, which is a relatively fast-flowing, narrow, deep warmwater current (Mizuno & White, 1983).It transports high amounts of nutrients and energy from the northward branch of North Equatorial Current off the Philippine coast (Kang et al., 2021) and turns northwards off the eastern coast of Taiwan where it attracts a variety of distinctive species (Denis et al., 2019).The possible reason why PR hotspots were not observed in these areas was that higher phylogenetic redundancy possibly exists in areas with more species richness (Ricotta et al., 2020).For example, higher phylogenetic redundancy was observed in fish communities in the coastal southwestern Atlantic, where species richness is higher than in the north (da Silva et al., 2021).This phylogenetic redundancy reduces the phylogenetic distinctiveness of species within the community, and therefore, no PR hotspots were detected., 1998).
The Yangtze River Estuary supports a wide range of fish species, including economically important species such as Chinese sturgeon and various types of carps, catfish and eels (Kindong et al., 2020).
The Yellow River Estuary supports fish species such as flounder and gobies.This estuary is influenced by the sediment-laden waters of the Yellow River, which can impact water clarity and nutrient levels.
The Pearl River Estuary supports a variety of fish species, including croakers, groupers, snappers and mudskippers (Kuang et al., 2021).
However, this estuary could not support as much FR and PR as the Yangtze River Estuary and the Yellow River Estuary possibly for the following reasons.First, the Yangtze River and the Yellow River carry a significant amount of sediments and nutrients to the estuary, creating a highly productive ecosystem and higher gradients of water temperature and salinity than the Pearl River Estuary (da Silva & Fabre, 2019).Therefore, the Pearl River Estuary may have a more limited range of habitats, which can result in lower FR and PR.
Second, human activities, such as pollution, overfishing and habitat destruction, can affect the fish species composition and richness in estuaries.
Species interactions possibly influenced the patterns of functional and phylogenetic rarity (Vitule et al., 2009).For instance, interspecific competition caused by the invasion of non-native species such as Mozambique tilapia Oreochromis mossambicus possibly excluded some native species with high distinctiveness in the Pearl River Estuary (Gu et al., 2019).Darwin's naturalization hypothesis is more likely to explain the establishment of non-native species, which indicates that non-native species most likely out-compete native species that are functionally similar (Skora et al., 2015).The prediction of the establishment of non-native species, for instance, with the conceptual model proposed by Skora et al. (2015), will possibly contribute to a better understanding of the patterns of functional and phylogenetic rarity.
Phylogeographic assembly rules are applicable at large spatial and temporal scales but environmental filtering for fish traits can function at smaller regional scales (Swenson, 2014).
Phylogeographic assembly processes refer to the restrictions in species composition that result from patterns of speciation and large-scale migration (Marske et al., 2013), whereas environmental filtering shapes the community structure as a result of interactions between species and environmental conditions (Keddy, 1992).For instance, sympatric speciation, that is, speciation without geographic isolation, is not common at small spatial scales for coral reef fish (Rocha & Bowen, 2008), and environmental filtering can dominate community assembly in bays or lagoons in coastal waters (Mouillot et al., 2007).Therefore, PR is dispersed at larger scales than FR.Many works aim to seek congruence between different facets of biodiversity because high management efficiency is possible if targeting the congruent areas with top taxonomic, functional, and phylogenetic diversity with limited conservation resources (Cadotte & Tucker, 2018).

| Mismatch between high rarity and current MPAs pattern
China's proposal of marine nature reserves began in the 1960s (Hu et al., 2020).Till 2020, the total covered areas of the MPAs in China had exceeded 50,000 km 2 (Pavoine et al., 2009).The matches between current MPAs and the hotspots of FR and PR in this work were significantly higher than random expectation, indicating that the spatial distributions of current MPAs were better designed than random from the perspective of FR and PR.However, these matches were still far beyond the average representations by global MPAs, that is, around 50% for global coastal fish species (Trindade-Santos et al., 2022), necessitating the expansion of current MPAs to better protect marine biodiversity (Li, Ren, et al., 2019).One possible reason was that many MPAs were not originally targeting the conservation of coastal fishes.For example, the vast areas of MPA in the Liaohe Estuary are primarily designed for protecting wetland habitats for rare water birds (Ministry of Ecology and Environment, 2017).
Other reasons might be attributed to the unbalanced spatial distribution of available conservation resources for historical reasons (Hu et al., 2020).For instance, governmental marine conservation institutes and research centers were first developed and more distributed in the north (Li & Liu, 2018).From a quantitative perspective, the number of MPAs in southern China was greater than that in northern China.From the perspective of covered area, the average size of MPAs in the south was much smaller than in the north.
For example, the average MPA area in Taiwan was about 40 km 2 (UNEP-WCMC & IUCN, 2022), but in the Liaoning Province was over 400 km 2 (Hu et al., 2020).
Previous evaluations of MPAs in China Seas exposed that the effectiveness of marine biodiversity conservation is weakened by the lack of systematic planning (Li, Zhang, et al., 2019).Only a small proportion of shared hotspots between FR and PR were observed in this work.Therefore, whether to rank higher priority for FR or PR becomes a pragmatic question for stakeholders when congruence can not be achieved (Devictor et al., 2010).Biogeographical patterns of rarity can be used to enhance the performance of the protected area networks (Albuquerque & Astudillo-Scalia, 2020).Integrative conservation prioritization by incorporating FR, PR and other factors is suggested for fish conservation.

| Limitations
Anthropogenic influences, especially marine pollution, overfishing and management strategies such as summer fishing moratoriums in coastal China Seas also impose certain impacts on coastal fish communities (Kang et al., 2021).These factors are difficult to be properly accounted for in current studies mainly because of scarce availability of useful historical data tracing back to the 1980s.The influence of these non-included factors remained in the residuals or was partly explained by the random effect of spatial autocorrelation.Although the explanatory and predictive powers of HMSC used were acceptable, inclusion of these variables in future studies will undoubtedly improve the performance of community modelling.

| CON CLUS IONS
Only considering threatened species in conservation practices will possibly omit functionally and phylogenetically rare species because they are not necessarily estimated as threatened on the IUCN Red List.Most FR hotspots do not overlap with PR hotspots.
FR hotspots geographically converge in the southern and eastern coastal seas of the Taiwan Province, the Yangtze River Estuary and the Yellow River Estuary.PR hotspots sporadically distribute in the coastal East China Sea, the Bohai Sea and the northern Yellow Sea.
Although the matches between MPAs and the hotspots of FR and PR were significantly higher than random expected, the low coverage indicates substantial insufficiency for the conservation of FR and PR.The expansion of MPAs is imminent for better protection of FR and PR.This work provides complementary information for fish biodiversity conservation, which will be preliminary for the development of a scientific and systematic conservation plan in coastal China Seas.

ACK N O WLE D G E M ENTS
This study was supported by the National Natural Science Foundation of China (no.U20A2087 and no.41976091).We ap-

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors have no conflicts of interest to declare.

PE E R R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/ api/ gatew ay/ wos/ peer-review/ 10. 1111/ ddi.13804 .

DATA AVA I L A B I L I T Y S TAT E M E N T
The data and R codes underlying this study are openly available in figshare at https:// doi.org/ 10. 6084/ m9.figsh are.22193 386.v1.

R E FE R E N C E S
Albuquerque, F., & Astudillo-Scalia, Y. (2020).The role of rarity as a surrogate of marine fish species representation.PeerJ, 8, e8373.https:// doi.org/ 10. 7717/ peerj.8373 strate) were provided by the State Oceanic Administration of China (2017).The north component of sea surface current (current), sea surface temperature (temperature) and sea surface salinity (salinity) were provided by National Science and Technology Infrastructure-National Marine Data Center (2022).Sea surface NPP (NPP) and sea surface pH (pH) were obtained from Copernicus Marine Environment Monitoring Service (2019).Monthly data of temperature, salinity, NPP and pH were averaged over the sampling period and sampling area.Present environmental variable layers for model prediction were obtained from Bio-ORACLE (Assis

Figure
Figure S5), where 'GTR' designates the general time-reversible model, meanwhile '+I' and '+ Γ' indicate that proportion of variable size and that gamma rate parameter get optimized in the GTR model.

Until 1
January 2022, the World Database on Protected Areas (WDPA) had only archived 15 Chinese MPAs (UNEP-WCMC & IUCN, 2022).In case of underestimation, we used the MPA data set recently compiled by Bohorquez et al. (2021), which recorded 273 operational MPAs in China Seas as of 1 October 2021.In case of overestimation, the integrity and authenticity of these WPAs were verified after cross-checking with the List of Nature Reserves in China (Ministry of Ecology and Environment, 2017).Polygons of these MPAs were then rasterized and cropped to align with the resolution of the research area.
Regulated by the national legitimated marine resources survey standard of China, the Technical Specification for Marine Fishery Resources Survey (Ministry of Agriculture of the People's Republic both functionally and phylogenetically rare, including milkfish Chanos chanos, six-lined trumpeter Helotes sexlineatus, Indo-Pacific tarpon Megalops cyprinoides, narrow-barred Spanish mackerel Scomberomorus commerson, Korean rockfish Sebastes schlegelii, and Bluespotted stargazer Xenocephalus elongatus.None of these species were threatened species.

F
Sensitivity analysis for functional (a) and phylogenetic distinctiveness (b), and distance decay plots for functional (c) and phylogenetic rarity (b).The X-axis in (a) and (b) indicates the distinctiveness calculated with all traits, and Y-axis indicates distinctiveness calculated after the removal of each trait one by one with replacement.Y-axis in (c) and (d) refers to pairwise Sorenson's dissimilarity of compositions of rare species.The red dotted lines in each plot are quantile regression lines y = x in (a) and (b) and y = 1.00 in (c) and (d).F I G U R E 3 Biogeographic patterns of functional rarity (a), phylogenetic rarity (b), and the shared hotspots of functional and phylogenetic rarity (c).The black circle labelled with YTR indicates the Yangtze River Estuary, and YR indicates the Yellow River Estuary.Legend colours in (a) and (b) indicate the standardized effect size (SES) of rarity.Hotspots in red colour in (c) indicate areas with SES > 2.0 for both functional and phylogenetic rarity.

F
Latitudinal density plots of the hotspots of functional rarity (a), phylogenetic rarity (b), and the shared hotspots of functional and phylogenetic rarity (c).F I G U R E 5 Spatial match and mismatch between the marine protected areas (MPAs) and the hotspots of functional rarity (a), hotspots of phylogenetic rarity (b), and the shared hotspots of functional and phylogenetic rarity (c).F I G U R E 6 Latitudinal density plots of the matched areas between the marine protected areas (MPAs) and the hotspots of functional rarity (a), hotspots of phylogenetic rarity (b), and the shared hotspots of functional and phylogenetic rarity (c).
preciate the data support by China National Science & Technology Infrastructure-National Marine Data Center and E.U.Copernicus Marine Service Information.

Models a AUC b AUC_CV c TjurR 2d TjurR 2 _CV e
The most frequent FR and PR occurred at approximately the same latitudes, that is, 30.93-31.48°N and 38.49-38.87°N.These two areas correspond to the two peaks of FR and PR hotspots in the Yangtze River Estuary, the Yellow River Estuary and their adjacent waters.The Yangtze River Estuary and the Yellow River Estuary are the first two biggest estuaries in China (Editorial Board of Bays of China