Spatiotemporal Patterns of Aquatic Product Risks in China

1. Introduction

Aquatic products refer to animals, plants, and microorganisms derived from aquatic sources, as well as various derivative products obtained through industrial processing [1]. As a vital source of protein, aquatic products supply nearly 20% of global animal protein consumption and have become a key component of dietary structures [2,3]. They are important sources of high-quality proteins, essential amino acids, fatty acids, and trace elements [4], characterized by high protein, low fat, easily digestibility, rich minerals, and delicious taste [5]. These products play a critical role in improving cardiovascular and cerebrovascular health, enhancing cognitive function, and promoting neural development, offering particular nutritional benefits for the cardiovascular health in the elderly and the intellectual development in children [6,7,8]. Asia accounts for over 70% of global aquatic product consumption [9]. As a major producer, consumer, and trader of aquatic products worldwide, China achieved a total aquatic product output of 71.1617 million metric tons in 2023, with total imports and exports reaching 10.5605 million metric tons and a total trade value of USD 44.237 billion, while its per capita annual consumption significantly exceeds the global average [10].

Due to the complexity of the aquatic product supply chain, which encompasses aquaculture, fishing, transportation, processing, and dining consumption, it faces multidimensional contamination risks. These risks present a complex feature by the coexistence of traditional and emerging pollutants, along with intertwined biochemical hazards. The quality and safety of aquatic products have attracted widespread attention, with food safety concerns primarily including heavy metal contamination, drug residues, and biological hazards [11,12]. Heavy metals are toxic and potentially carcinogenic. They accumulate in aquatic systems and cause bioaccumulation even at low exposure concentrations, and can damage multiple organs, including the nervous system, liver, lungs, kidneys, stomach, skin, and reproductive system [13]. For example, arsenic can cause reproductive disorders and nervous system damage; lead may impair the nervous system, disrupt brain activity, and impair bone marrow function; cadmium is linked to adverse effects on bones, kidneys, and reproductive systems [14,15]. As aquaculture expands in scale, stocking densities have increased, which in turn leads to frequent outbreaks of various diseases. This has also been coupled with increased overuse of drugs, non-compliant medication, and the use of banned substances, resulting in the residues of these substances in aquatic products and widespread concern [16,17]. Common drug residues in aquatic products primarily include nitrofuran metabolites, chloramphenicol, quinolones, sulfonamides, and tetracyclines [18,19,20]. They pose significant health risks to people, including the development of antibiotic resistance, carcinogenic effects, teratogenic effects, allergic reactions, and disruption of normal gut microbiota [21,22]. In response to the aforementioned risks, scholars have progressively developed a multidimensional assessment framework. Chen et al. [23] comprehensively evaluates heavy metal contamination in aquatic products and associated health risks by integrating single-factor pollution indexes, composite pollution indexes, and safety assessments. Ku et al. [24] proposes a decision tree model for heavy metal hazards to assess the potential health risks of 12 heavy metals in seven categories of aquatic products. The study integrates the statistical data from the fisheries in Taiwan and national food consumption data, and by analyzing heavy metal concentrations in 556 cooking aquatic product samples, it conducts risk assessments using estimated daily intake (EDI) and hazard quotient (HQ). Liu et al. [25]determines indicator weights using the Analytic Hierarchy Process (AHP), and calculates risk indices using comprehensive index and multiplicative synthesis. Taking data from Guangzhou as an example, the annual risk indice and grades are calculated to assess food safety conditions and validate the index methods. The final risk assessment indicator system comprises 5 primary indicators, 13 S-level indicators, and 13 chemical contaminants, covering hazards, vulnerability, monitoring, consumption patterns, public sentiment, and regulatory oversight. However, existing studies are largely limited to static assessments and lack dynamic analysis of spatiotemporal heterogeneity and synergistic effects of pollutants, which constrains the precision of risk warnings.

Therefore, the present study proposed a three-dimensional assessment framework integrating time, space, and risk substances. Based on the supervision data of 1.04 million batches of aquatic products from China's 31 provinces between 2021 and 2023, four categories of indicators were consolidated: non-compliance rate, detection rate, qualification degree, and hazard degree. Methodologically, the entropy weighting method was employed to objectively assign weights, overcoming the limitations of subjective preferences inherent in traditional AHP. By combining the TOPSIS model, the framework quantified the priority of risk substances, and achieved dynamic ranking by calculating the proximity of each sample to the ideal solution. This model aimed to reveal the spatiotemporal distribution patterns of risk substances in aquatic products and the risk transmission pathways within the supply chain. It provides algorithmic support for constructing a comprehensive, intelligent, and dynamic monitoring network spanning the entire chain “from farm to fork”, offering scientific basis for the safety supervision of aquatic products. It drives the transformation of aquatic product safety governance from an experience-driven approach to a data-driven paradigm, facilitating the transition from extensive regulation to precision-based management.

2. Materials and Methods

2.1. Data Sources and Study Methods

The data in this study were sourced from the publicly released supervision and inspection data of aquatic products by the market supervision administrations in various provinces, municipalities, and autonomous regions (excluding Taiwan, Hong Kong, and Macao) between 2021−2023. Based on the research objectives, the supervision and inspection data of aquatic products were organized, and nine items were ultimately retained: sampling location, sampled province, sampling date, sampling stage, regional type of sampled entity, inspection items, sampling results, and determination of sampling results. Four types of indicators were selected: non-compliance rate, detection rate, qualification degree, and hazard degree [26]. The rankings of risk substances were determined using the entropy weight TOPSIS method, and risk classification was performed based on the Pareto principle.

2.2. Analytical Methods

2.2.1. Construction of Risk Classification Model Based on the Entropy Weight TOPSIS Method

The entropy weight TOPSIS model was an integrated approach that combines the entropy weight method with the TOPSIS technique. The entropy weight method determined the weights of various indicators by analyzing the information they provide. Subsequently, the TOPSIS method was used to calculate the distances between each indicator and its ideal solution, thereby assessing the proximity of each evaluation object to the optimal level. This approach mitigated the subjectivity associated with manual weight assignment, reducing the potential for bias from human judgment. Consequently, it yielded evaluation results that were more aligned with objective reality. The model boasted several advantages such as strong interpretability, intuitive geometric interpretation, and flexible computation [27,28].

2.2.1.1. Calculation of Weights by the Entropy Weight Method

For n samples and k indicators, X_ij represented the jth indicator for the ith sample (where i = 1, 2, ..., n; j = 1, 2, ..., k), ${\mathrm{D}}_{\mathrm{ij}}$ denoted the standardized value of each indicator. The matrix X was obtained from i risk factors and j indicators, as shown in Equation 1.

$$\left.\mathrm{X}=\left[\begin{array}{cccc}{\mathrm{x}}_{11}& {\mathrm{x}}_{12}& \dots & {\mathrm{x}}_{1\mathrm{j}}\\ {\mathrm{x}}_{21}& {\mathrm{x}}_{22}& \dots & {\mathrm{x}}_{2\mathrm{j}}\\ \dots & \dots & \dots & \dots \\ {\mathrm{x}}_{\mathrm{i}1}& {\mathrm{x}}_{\mathrm{i}2}& \dots & {\mathrm{x}}_{\mathrm{ij}}\end{array}\right.\right]$$

(1)

Perform dimensionless processing on the original matrix to eliminate the impact of differing units of measurement, yielding the matrix ${\mathrm{D}}_{\mathrm{ij}}$ (Equation 2). Max (X_j) was the maximum value for each indicator, while min (X_j) was the minimum value. Should zero values appear in the standardised matrix after processing, to ensure meaningful data handling, all dimensionless data shall be shifted by a minimum unit value.

$${\mathrm{D}}_{\mathrm{i}\mathrm{j}}=\left\{\begin{array}{c}{\displaystyle \frac{{\mathrm{x}}_{\mathrm{i}\mathrm{j}}-\mathrm{m}\mathrm{i}\mathrm{n}({\mathrm{x}}_{\mathrm{j})}}{\mathrm{m}\mathrm{a}\mathrm{x}\left({\mathrm{x}}_{\mathrm{j}}\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left({\mathrm{x}}_{\mathrm{j}}\right)}} \mathrm{Positive} \mathrm{index}\\ {\displaystyle \frac{\mathrm{m}\mathrm{a}\mathrm{x}\left({\mathrm{x}}_{\mathrm{j}}\right)-{\mathrm{x}}_{\mathrm{j}}}{\mathrm{m}\mathrm{a}\mathrm{x}\left({\mathrm{x}}_{\mathrm{j}}\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left({\mathrm{x}}_{\mathrm{j}}\right)}} \mathrm{Negative} \mathrm{index}\end{array}\right.$$

(2)

Calculate the feature weight or contribution of the i-th item under the j-th indicator (Equation 3):

$${\mathrm{P}}_{\mathrm{ij}}={\displaystyle \frac{{\mathrm{D}}_{\mathrm{ij}}}{{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{D}}_{\mathrm{ij}}}}$$

(3)

The information entropy was calculated (Equations 4 and 5):

$${\mathrm{e}}_{\mathrm{j}}=-\mathrm{k}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{P}}_{\mathrm{ij}}\ln {\mathrm{P}}_{\mathrm{ij}}$$

(4)

$$\mathrm{K}={\displaystyle \frac{1}{\ln (\mathrm{n})}}$$

(5)

The weights were calculated (Equations 6 and 7):

$${\mathrm{g}}_{\mathrm{j}}=1-{\mathrm{e}}_{\mathrm{j}}$$

(6)

$${\mathrm{W}}_{\mathrm{j}}={\displaystyle \frac{{\mathrm{g}}_{\mathrm{j}}}{\sum {\mathrm{g}}_{\mathrm{j}}}}$$

(7)

2.2.1.2. Calculation of Relative Proximity Using the TOPSIS Method

Weighted Standardized matrix, multiply matrix D_ij by the obtained indicator weights W_j to obtain matrix U_ij as shown in Equation 8:

$${\mathrm{U}}_{\mathrm{ij}}={\mathrm{D}}_{\mathrm{ij}}\times {\mathrm{W}}_{\mathrm{j}}$$

(8)

The optimal value vector Ui+ and the suboptimal value vector Ui- for each indicator were determined (Equations 9 and 10):

$${\mathrm{U}}_{\mathrm{i}}^{+}=\max \left\{{\mathrm{U}}_{1\mathrm{j}} ，{\mathrm{U}}_{2\mathrm{j}} ，\dots ，{\mathrm{U}}_{\mathrm{nj}} \right\}$$

(9)

$${\mathrm{U}}_{\mathrm{i}}^{-}=\mathrm{m}\mathrm{in}\left\{{\mathrm{U}}_{1\mathrm{j}} ，{\mathrm{U}}_{2\mathrm{j}} ，\dots ，{\mathrm{U}}_{\mathrm{nj}} \right\}$$

(10)

The distances of each evaluation object from the optimal value vector and the suboptimal value vector were calculated (Equations 11 and 12):

$${\mathrm{D}}_{\mathrm{i}}^{+}=\sqrt{{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{j}}\times {\left({\mathrm{U}}_{\mathrm{ij}}-{\mathrm{U}}_{\mathrm{j}}^{+}\right)}^2}$$

(11)

$${\mathrm{D}}_{\mathrm{i}}^{-}=\sqrt{{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{j}}\times {\left({\mathrm{U}}_{\mathrm{ij}}-{\mathrm{U}}_{\mathrm{j}}^{-}\right)}^2}$$

(12)

The relative proximity was calculated, with values ranging from [0,1], where a larger value indicated greater risk associated with the hazard. The equation is shown in Equation 13:

$${\mathrm{C}}_{\mathrm{i}}={\displaystyle \frac{{\mathrm{D}}_{\mathrm{i}}^{-}}{{\mathrm{D}}_{\mathrm{i}}^{+}+{\mathrm{D}}_{\mathrm{i}}^{-}}}$$

(13)

2.2.2. Spatial Autocorrelation Analysis

2.2.2.1. Global Spatial Autocorrelation Analysis

Spatial analysis was adopted to identify clustered areas and observe geographic spatial variations. By describing the spatial distribution characteristics of non-compliance rates for high-risk and higher-risk substances in aquatic products, we utilized Moran's I statistic, which ranged between −1 and 1. Values close to 1 indicated positive spatial autocorrelation, values close to −1 indicated negative spatial autocorrelation, and 0 indicated random distribution. Moran's I was calculated and hotspot analysis was performed to determine the location of clusters [29]. Typically, a weight (W_ij) was associated with each pair (X_i, X_j) to quantify spatial patterns, as shown in Equation 14:

$$\mathrm{I}={\displaystyle \frac{\mathrm{n}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{ij}}\left({\mathrm{x}}_{\mathrm{i}}-\stackrel{\mathrm{‾}}{\mathrm{x}}\right)\left({\mathrm{x}}_{\mathrm{j}}-\stackrel{\mathrm{‾}}{\mathrm{x}}\right)}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{ij}}\left({\mathrm{x}}_{\mathrm{i}}-\stackrel{\mathrm{‾}}{\mathrm{x}}\right)}}$$

(14)

In the formula, n represented the number of provinces, municipalities, and autonomous regions in China; X_i and X_j denoted the compliance rates of a certain risk factor in regions i and j across the country, respectively; $\stackrel{\mathrm{‾}}{\mathrm{x}}$ represented the average compliance rates; W_ij was the spatial weighting matrix.

2.2.2.2. Local Spatial Autocorrelation Analysis

Local Moran's I was employed to determine local spatial autocorrelation, identifying spatial clusters with similar neighboring characteristics and outliers, as shown in Equation 15:

$$\mathrm{I}={\displaystyle \frac{\left({\mathrm{x}}_{\mathrm{i}}-\stackrel{\mathrm{‾}}{\mathrm{x}}\right)}{{\mathrm{S}}^2}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{ij}}\left({\mathrm{x}}_{\mathrm{j}}-\stackrel{\mathrm{‾}}{\mathrm{x}}\right)$$

(15)

In the formula, S² denoted the variance of X_i or X_j, while the other components were interpreted similarly to the global Moran's I.

2.2.3. Spatiotemporal Analysis

Through statistical analysis, the non-compliance rates of aquatic products across provinces and each month were yielded, and thus the risk adjustment factors for different provinces and months were obtained, as shown in Equations 16−18 [30]. In the formula, P was the risk adjustment factor for a province; Pf was the number of non-compliant cases in that province; Qf represented the total number of non-compliant cases for food products; M denoted the risk adjustment factor for a month; Mf signified the total number of non-compliant cases across all provinces for that month. The maximum relative proximity was taken as the baseline risk for the food, to which a risk adjustment factor was added to calculate the final spatiotemporal risk value for the aquatic product.

$$\mathrm{P}=1+{\displaystyle \frac{{\mathrm{P}}_{\mathrm{f}}}{{\mathrm{Q}}_{\mathrm{f}}}}$$

(16)

$$\mathrm{M}=1+{\displaystyle \frac{{\mathrm{M}}_{\mathrm{f}}}{{\mathrm{Q}}_{\mathrm{f}}}}$$

(17)

$$\mathrm{R}={{\mathrm{C}}_{\mathrm{i}}}_{\max }\times \mathrm{P}\times \mathrm{M}$$

(18)

2.2.4. Statistical Analysis

In this study, the inspection data of aquatic products from 2021 to 2023 were organized by Microsoft Excel 2019. Origin 2021 and ArcGIS 10.8 were used to create maps.

3. Results

3.1. Detection Results of Risk Substances in Aquatic Products

After processing the inspection data of aquatic products from 2021 to 2023, 1,044,312 valid sample data were obtained. The inspection items encompassed six categories: heavy metals, additives, veterinary drugs, prohibited drugs, quality indicators, and organic pollutants, with 34 items in total. These inspections covered both distribution and catering stages. The non-compliance rate, detection rate, qualification degree, and hazard degree for all inspection items across different stages were determined. Using the entropy weight method, the weights for the four indicators—non-compliance rate, detection rate, qualification degree, and hazard degree—were obtained for both the distribution and catering stages. The weight coefficients w for hazard degree, qualification degree, detection rate, and non-compliance rate in the distribution stage were 5.46%, 17.18%, 35.45%, and 41.91%, respectively, while those in the catering stage were 5.72%, 8.14%, 36.89%, and 49.25%, respectively. The TOPSIS model was then employed to comprehensively calculate the proximity of each sample to the ideal solution, achieving the dynamic ranking (Table 1). Among the risk substances in the distribution stage and catering stage, those with higher proximity rankings were cadmium, enrofloxacin, total volatile basic nitrogen.

Relative proximity was classified according to the Pareto principle, with the classification range shown in Table 2, and specific classification results are presented in Table 3. In the distribution stage, cadmium, enrofloxacin, and total volatile basic nitrogen were classified as the highest risk substances; sulfur dioxide, diazepam, methylmercury, malachite green, furazolidone metabolites, chloramphenicol, nitrofurazone metabolites, sodium pentachlorophenolate, metronidazole, nitrofurantoin metabolites, furaltadone metabolites, sarafloxacin, and ofloxacin were categorized as the higher-risk substances; polychlorinated biphenyls, danofloxacin, flumequine, difloxacin, oxolinic acid, inorganic arsenic, chromium and histamine fell into the the medium-risk category; lead, deltamethrin, cypermethrin, trimethoprim, florfenicol, and sulfonamides were categorized as the lower-risk substances; while pefloxacin, oxytetracycline/chlortetracycline/tetracycline (sum), norfloxacin, and lomefloxacin were assigned to the lowest-risk category. In the catering stage, cadmium, enrofloxacin, and total volatile basic nitrogen were classified in the highest-risk category; methylmercury and sodium pentachlorophenolate posed an elevated risk. Inorganic arsenic, chromium, nitrofurazone metabolite, furazolidone metabolite, diazepam, malachite green, chloramphenicol, metronidazole, lead, furaltadone metabolites, nitrofurantoin metabolites, and sarafloxacin were categorized as the medium-risk substances. Ofloxacin, oxolinic acid, difloxacin, pefloxacin, norfloxacin, danofloxacin, flumequine, lomefloxacin, deltamethrin, polychlorinated biphenyls , cypermethrin, florfenicol, oxytetracycline/chlortetracycline/tetracycline (sum), sulfur dioxide, and histamine were considered the lower-risk compounds. Trimethoprim and sulfonamides (total amount) were identified as the lowest-risk substances. Overall, the risk of various risk substances in aquatic products was higher in the distribution stage than that in the catering stage.

An analysis of non-compliant aquatic products (Figure 1) showed that the highest non-compliance rate was observed in seawater crabs (3.19%), followed by other aquatic products (2.25%), and the lowest rate was freshwater crabs (0.04%). An examination was conducted on the non-compliant items detected in different categories of aquatic products. The proportion of each non-compliant item by category is illustrated in Figure 2. The main non-compliant items in freshwater crab were heavy metals and prohibited drugs; in freshwater shrimp, shellfish, and other aquatic products, the main items were heavy metals, veterinary drugs, and prohibited drugs; in seawater fish, heavy metals, veterinary drugs, prohibited drugs, and quality indicators were the primary concerns; seawater shrimp were mainly found to have heavy metals, veterinary drugs, prohibited drugs, and additives; seawater crab was mainly associated with heavy metals; in freshwater fish, veterinary drugs and prohibited drugs were mainly reported. The distribution of non-compliant items detected in different categories of aquatic products is presented in Table 4. For seawater crabs, cadmium was the primary concern, while other aquatic products predominantly showed issues with enrofloxacin, furazolidone metabolites, and nitrofurazone metabolites; for seawater shrimp, it was mainly cadmium and enrofloxacin; freshwater fish was mostly associated with enrofloxacin, diazepam, and malachite green, while shellfish primarily contained chloramphenicol and cadmium; both freshwater shrimp and seawater fish had significant levels of enrofloxacin and furazolidone metabolites, whereas freshwater crabs primarily contained cadmium, nitrofurazone metabolites, and furazolidone metabolites.

An analysis of the non-compliant items in aquatic products revealed that veterinary drugs accounted for the largest proportion at 53.74%, followed by heavy metals at 25.38%. No non-compliant cases were detected for organic pollutants. Based on the risk classification of aquatic products and the non-compliant items in aquatic products (Figure 3), the primary focus should be on veterinary drugs and heavy metals. Among veterinary drugs (Figure 4a), particular attention should be paid to the illegal use of enrofloxacin, which accounted for 92.96% of the non-compliant cases. Regarding heavy metals, only cadmium was found to be non-compliant. In combination with the risk classification, significant attention should be paid to cadmium and methylmercury. The inspection items under prohibited drugs (Figure 4b) all presented relatively high risks and require close attention. Non-compliance in quality indicators was detected for histamine and total volatile basic nitrogen; considering the risk classification, total volatile basic nitrogen requires particular attention. In the category of additives, only sulfur dioxide was found to be non-compliant, which should be closely monitored in light of the risk classification.

3.2. Temporal Distribution of Risky Substances in Aquatic Products

A temporal analysis of the non-compliant cases in the aquatic products from 2021 to 2023 (Figure 4. Non-compliance of inspection tests for veterinary drugs(a) and prohibited drugs (b)) showed a general trend of gradual increase in non-compliance across different months, with a clear seasonal pattern. The non-compliance rate was higher in summer and autumn than in winter and spring. Overall, the non-compliance rate in the second half of the year was higher than that in the first half, exhibiting a gradual upward trend, though the overall fluctuation remained relatively small. In 2021, the lowest non-compliance rate was 0.12% in February, and the highest was 1.73% in December, indicating significant fluctuation. In 2022, the lowest rate was 0.44% in December, and the highest was 1.05% in October. In 2023, the non-compliance rate was lowest in December at 0.23% and highest in November at 0.58%, demonstrating relatively minor fluctuation.

Figure 5. Changes in the non-compliance rate of aquatic products by month.

Figure 5. Changes in the non-compliance rate of aquatic products by month.

3.2.1. Veterinary Drugs

A temporal analysis was conducted on the non-compliance rate of veterinary drugs in aquatic products (Figure a). The main non-compliant substances comprised five categories: enrofloxacin, sulfonamides, trimethoprim, florfenicol, and ofloxacin. Among these, enrofloxacin was non-compliant throughout the entire year, with the highest non-compliance rate observed in December at 4.83%. Trimethoprim also exhibited its peak non-compliance rate of 2.50% in December. Florfenicol reached a non-compliance rate of 1.03% in May, a level that warrants particular attention. Other veterinary drugs demonstrated sporadic non-compliance, with relatively stable non-compliance rates.

3.2.2. Heavy Metals

A temporal analysis of non-compliant cases related to heavy metals in aquatic products (Figure b) indicated that cadmium was the only heavy metal that was non-compliant. The non-compliance rate fluctuated significantly throughout the year, showing an upward trend. It peaked at 5.63% in December and was lowest at 0.92% in February. Continuous attention to cadmium is required.

3.2.3. Prohibited Drugs

A temporal analysis was conducted on the non-compliance rate of prohibited drugs in aquatic products (Figure c). The main non-compliant substances were furazolidone metabolites, malachite green, chloramphenicol, nitrofurazone metabolites, sodium pentachlorophenolate, and diazepam. Among them, diazepam had the highest non-compliance rate in November at 1.14%; chloramphenicol peaked in May at 0.72%; malachite green reached its highest level of 0.71% in April. Other prohibited drugs showed sporadic non-compliance with relatively stable rates. Overall, the non-compliance of prohibited drugs was minimally influenced by seasonal variations.

3.2.4. Quality Indicators

A temporal analysis was conducted on the non-compliance of quality indicators in aquatic products (Figure d), and the only non-compliant item among the quality indicators was total volatile basic nitrogen. Overall, the risk of non-compliance for total volatile basic nitrogen was relatively low. In May, the highest non-compliance rate was 0.43%, which may be attributed to temperature elevation causing excessively high values in quality indicators, requiring targeted monitoring.

3.2.5. Additives

A temporal analysis of non-compliance related to additives in aquatic products (Figure e) showed that the only non-compliant additive was sulfur dioxide. The non-compliance rate fluctuated significantly, reaching a peak of 3.78% in May. The non-compliance rate remained relatively high from September to November, at 1.44%, 1.89%, and 1.51%, respectively.

Figure 6. Trend of non-compliance rate of veterinary drugs (a), heavy metals (b), prohibited drugs (c), quality indicators (d), additives (e) in aquatic products.

Figure 6. Trend of non-compliance rate of veterinary drugs (a), heavy metals (b), prohibited drugs (c), quality indicators (d), additives (e) in aquatic products.

3.3. Spatial Distribution of Risk Substances in Aquatic Products

An analysis was conducted on the non-compliant sample data of aquatic products across provinces during 2021−2023. A classification map of the non-compliance rate was generated using the equal interval method in ArcGIS (Figure 7), illustrating the overall national non-compliance rate of aquatic products (a), and non-compliance rate of heavy metals (b), quality indicators (c), additives (d), veterinary drugs (e), and prohibited drugs (f). In the analysis of the overall non-compliance rate across cities, a general pattern was observed, with higher rates in the east and lower rates in the west, and coastal areas showing higher rates than inland areas. High-risk areas were concentrated in coastal cities of East and South China (such as the Yangtze River Delta and Pearl River Delta), while provinces in the central and western regions (such as Qinghai and Tibet) were classified as low-risk areas. The non-compliance rates in Zhejiang, Chongqing, and Jiangxi were higher than those in other provinces, at 1.66%, 1.52%, and 1.48%, respectively, with Zhejiang being the province with the highest number of non-compliant aquatic products nationwide. An analysis of non-compliance rate of heavy metals across cities indicated that high-risk areas were concentrated in northern cities, medium risk in coastal cities, and low risk in the northwest region. Regarding the non-compliance rate of quality indicators across cities, there was a scattered distribution without significant large-scale clustering. High-risk areas were sporadically distributed in the northeastern industrial zone, the North China Plain, and the Sichuan Basin, while most regions across the country exhibited low risk. For the non-compliance rate of additives across cities, it revealed isolated high-value points without marked large-scale clustering and no evident spatial patterns. By analysis of the non-compliance rate of veterinary drugs across cities, it was found that non-compliance mainly occurred in aquaculture-intensive areas in South and East China (such as Guangdong, Fujian, and Jiangsu), while the risk was relatively low in northern pastoral areas and southwestern mountainous regions. An analysis of the non-compliance rate of prohibited drugs across cities displayed that high-risk areas were identified in provinces (autonomous regions and municipalities) such as Beijing, Sichuan, Guangxi, and Zhejiang, while most cities across the country exhibited relatively low risks.

3.3.1. Spatial Autocorrelation Analysis

According to Table 3 that classifies the risk substances in aquatic products, the hazards with the highest and higher risks were selected as the focus. Among them, there were no non-compliant samples for methylmercury and sarafloxacin; furaltadone metabolites was found to have non-compliance issues only in Jiangxi; nitrofurantoin metabolites showed non-compliance only in Inner Mongolia and Jiangxi. Therefore, we analyzed the spatial distribution patterns of the following risk substances: cadmium, enrofloxacin, total volatile basic nitrogen, sulfur dioxide, diazepam, malachite green, furazolidone metabolites, chloramphenicol, nitrofurazone metabolites, sodium pentachlorophenate, metronidazole, and ofloxacin.

3.3.1.1. Global Spatial Autocorrelation Analysis

A spatial autocorrelation analysis was performed on the non-compliance rate of 12 risk substances in aquatic products across provinces in China. The global Moran's I values for each risk substance in aquatic products were obtained (Table 5), among which the p-values for total volatile basic nitrogen, sulfur dioxide, diazepam, malachite green, nitrofurazone metabolites, sodium pentachlorophenate, metronidazole, and ofloxacin were all greater than 0.05, suggesting no global spatial autocorrelation in the non-compliance rate of these eight risk substances in aquatic products across provinces. However, local spatial correlation may exist. The p-values for cadmium, enrofloxacin, furazolidone metabolites, and chloramphenicol in aquatic products across provinces were less than 0.05, indicating that these four risk substances exhibited global spatial correlation in non-compliance rate across provinces, demonstrating a clustered distribution pattern.

3.3.1.2. Local Spatial Autocorrelation Analysis

When there is a global Moran's I index with statistical significance (p<0.05), further analysis of the local Moran's I index can be conducted to investigate whether there is regional spatial clustering in the study subjects. A local spatial autocorrelation analysis was performed on cadmium, enrofloxacin, furazolidone metabolites, and chloramphenicol in aquatic products across provinces in China. First of all, the p-values of the local Moran's I index were analyzed. If the p-value is less than 0.05, it indicates the presence of a regional spatial correlation; conversely, if the p-value is greater than 0.05, it indicates the absence of spatial correlation (Table 6). A high-high clustering signifies that both the province and its neighboring provinces were at a “high level”; a low−high clustering meant that the province was at a “low level”, but its neighboring provinces were at a “high level”; a low−low clustering indicated that both the province and its neighboring provinces were at a “low level”; a high−low clustering denoted that the province was at a “high level”, but its neighboring provinces were at a “low level”. For the Moran scatter plot, the X-axis (z-score of outliers) represented the distance between data values and the mean; therefore, data points located further to the right indicated higher non-compliance rate. The Y-axis (Spatial Lag) represented the spatial lag value, with higher values indicating that the neighbors (surrounding areas) of the study subject had higher non-compliance rate.

From Figure 8, the compliance rate for cadmiumnon exhibited significant low−low clustering in Chongqing, Sichuan, Guizhou, Yunnan, Tibet, Shaanxi, Gansu, Qinghai, and Xinjiang. In contrast, it showed significant high−high clustering in Beijing, Tianjin, Hebei, Inner Mongolia, Liaoning, Jilin, Heilongjiang, and Shandong. Obvious low−high clustering was observed in Shanxi. For enrofloxacinnon, the compliance rate displayed significant low−low clustering in Hebei, Shanxi, Inner Mongolia, and Liaoning, and significant high−high clustering in Shanghai, Jiangsu, Zhejiang, Anhui, Fujian, Jiangxi, Hubei, and Hunan. Additionally, marked low−high clustering was observed in Guizhou. Regarding the non-compliance rate of furazolidone metabolites, there was remarkable low−low clustering in Hebei, Inner Mongolia, and Jilin, and significant high-high clustering in Fujian, Shanghai, Guangdong, Hainan, and Zhejiang, while remarkable low−high clustering was observed in Jiangxi and Jiangsu. For chloramphenicolnon, significant low−high clustering was detected in Yunnan and Guizhou, significant high−low clustering in Sichuan, and significant high−high clustering in Guangdong and Hainan.

3.3.2. Sampling Location Analysis

The sampling locations for aquatic products included urban areas, township areas, tourist attractions (urban), tourist attractions (townships), vicinity of schools (townships), and vicinity of schools (urban). Non-compliant aquatic product samples were primarily found in urban areas (Figure 9), accounting for 81.34%, while township areas accounted for 16.04%.

3.3.3. Sampling Venue Analysis

Two sampling stages were included in the inspection of aquatic products: distribution and catering. The highest incidence of non-compliance occurred in the distribution stage, accounting for 85.71%, followed by the catering stage, which accounted for 14.29% (Figure 10). An analysis of sampling venue (Figure 11) showed that in the distribution stage, non-compliant cases mainly occurred in farmers' markets and supermarkets, accounting for 38.75% and 26.48%, respectively. Fewer non-compliant cases were found in wholesale markets, other locations, vegetable markets, online purchase, small grocery, and shopping malls, accounting for 10.39%, 11.16%, 5.55%, 5.41%, 1.99%, and 0.27% respectively. Therefore, special attention should be paid to the quality safety of aquatic products in farmers' markets and supermarkets. In the catering stage, nearly half of the non-compliant cases occurred in medium-sized restaurants, accounting for 45.69%. The remaining non-compliant cases occurred in small restaurants, large restaurants, other locations, extra-large restaurants, snack bars, and catering service provider, which accounted for 33.68%, 15.81%, 1.95%, 1.33%, 0.51%, and 1.03%, respectively. Special attention should be given to the quality safety of aquatic products in small restaurants and medium-sized restaurants.

3.4. Spatiotemporal Analysis of Risk Substances in Aquatic Products

According to the evaluation results by the entropy weight TOPSIS method, the maximum relative proximity value ${\mathrm{C}}_{\mathrm{i}}$ of cadmium was 0.707, which served as the baseline risk for aquatic products. After incorporating the risk adjustment factor, the final spatiotemporal risk values of aquatic products were calculated, as shown in Table 7. Based on the monthly final spatiotemporal risk values for each province, a heatmap(Figure ) was generated to illustrate the risk variations of aquatic products across different provinces and time periods by changes in color intensity, with red indicating high risk and blue indicating low risk. The figure showed that the range of risk values was concentrated between 0.71 (low risk) and 0.92 (high risk), with a relatively narrow span, indicating that while there were regional differences in the risk of aquatic products nationwide, the overall risk was controllable. From August to November, the risks were higher across provinces (municipalities, autonomous regions), while from December to April, the overall risk was lower. Among them, Chongqing, Jiangxi, and Zhejiang exhibited higher risks throughout the year, whereas Tibet, Ningxia, and Gansu demonstrated lower annual risks.

Figure 12. Heatmap of changes in aquatic products risk in different provinces and months. Note: In the figure, values ranging from 0.71 to 0.92 indicate low−high risk.

Figure 12. Heatmap of changes in aquatic products risk in different provinces and months. Note: In the figure, values ranging from 0.71 to 0.92 indicate low−high risk.

4. Discussion

4.1. Core Risk Substances and Their Health Risks

This study determined that high-risk substances in aquatic products include cadmium, enrofloxacin, and total volatile basic nitrogen by the entropy weight TOPSIS method. While substances posing higher risks include sulfur dioxide, diazepam, methylmercury, malachite green, nitrofuran compounds, chloramphenicol, sodium pentachlorophenate, metronidazole, sarafloxacin, and ofloxacin. A study analyzing the monitoring data of heavy metals in Chinese aquatic products has found that among different aquatic products-heavy metals combinations, the sea crab-cadmium combination posed the highest risk[31], which aligns with the findings of the present study that seawater crabs had the highest proportion of non-compliance rates, with cadmium posing the highest risk. Ding et al.[32] conducted testing on 300 batches of commercially available aquatic product samples for veterinary drug residues and heavy metals, and revealed that key risk parameters included cadmium, malachite green, enrofloxacin, nitrofurazone metabolite, and ofloxacin, which are consistent with the results of the present study. This confirms the reliability of the aquatic product classification model developed herein. Based on the risk classification, it is essential to enhance inspections for cadmium, enrofloxacin, total volatile basic nitrogen, sulfur dioxide, diazepam, and methylmercury in the distribution stage of aquatic products, and to increase inspections for cadmium, enrofloxacin, volatile basic nitrogen, methylmercury, and sodium pentachlorophenate in the catering stage.

Among heavy metals, cadmium warrants particular attention. Cadmium in the environment exhibits a bioamplification effect, which can accumulate through food sources and subsequently transfer to humans via the food chain. Fish, crustaceans, and bivalves are the most common aquatic products, and are also among the main sources of cadmium exposure[33,34]. These aquatic products can accumulate cadmium in their tissues through body surface permeation and ingestion, and absorb cadmium from water[35]. A study from Laizhou Bay, China, indicates that cadmium is the second most abundant metal and the most easily accumulated element in seafood, with crabs exhibiting a much higher cadmium enrichment capacity than other aquatic products[36]. Heavy metals can move upward along the food chain, posing potential health risks to humans[1]. Since food is an important carrier for the transfer of heavy metals from the environment to the human body, it is crucial to control the content of heavy metals in aquatic products and thus reduce harm to human health. Long-term consumption of high-risk species should be avoided, and it is advisable to choose freshwater fish or small-sized shellfish with low bioaccumulation capacity for heavy metals. Soil and water quality investigations for heavy metal pollution should be conducted in aquaculture areas to trace and remediate agricultural land contamination. Strict limitations should be imposed on the discharge of industrial wastewater and agricultural non-point source pollution to reduce heavy metals in water bodies, thereby lowering the bioaccumulation risk of aquatic products from the source[37].

Among veterinary drugs, enrofloxacin should be a primary focus, while particular attention should be given to high-risk compounds such as diazepam, furazolidone metabolites, malachite green, chloramphenicol, and sodium pentachlorophenate among prohibited substances. Enrofloxacin, due to its broad-spectrum antibacterial properties, excellent pharmacokinetic characteristics, low risk of cross drug resistance, and low cost, is widely used in the prevention and control of diseases in livestock, poultry farming, and aquaculture. It is one of the most frequently detected veterinary drugs in aquatic products[38,39]. The Ministry of Agriculture and Rural Affairs of China has set the maximum residue limit of enrofloxacin in fish meat at 100 μg/kg[40]. Diazepam[41], chloramphenicol[42], malachite green[20], sodium pentachlorophenate[43], and nitrofuran compounds[44] are drugs that have been banned from detection in aquatic products, but illegal use still exists, resulting in the residue or even excessive levels of harmful substances in aquatic products, posing health risks. To reduce the occurrence of veterinary drug residues and non-compliance situations for prohibited drugs in aquatic products, enterprises need to optimize farming management by maintaining water quality, ensuring reasonable feeding practices, and implementing disease prevention strategies, thus reducing drug use, and strictly adhering to withdrawal periods[45].Governments can establish cross-provincial and even cross-border traceability mechanisms to enhance information sharing between production and consumption areas, thereby preventing regulatory gaps. Promoting the aquatic products traceability platforms that record the entire process of medication use in aquaculture, transportation, and sales, will facilitate the rapid recall of problematic products.

Total volatile basic nitrogen should be a key focus in quality indicators. It is the primary indicator reflecting the freshness of aquatic products[46], enabling timely determination of whether aquatic products have been contaminated during transportation, food service, and other stages. It is commonly used as an important monitoring indicator for the freshness and safety of aquatic products, and changes in its concentration can objectively reflect the spoilage and deterioration process of aquatic products[47,48]. To ensure the freshness and quality safety of aquatic products, quality control should be implemented across the entire industry chain. At the raw material stage, a stable and reliable supplier selection mechanism should be established, along with rapid detection of total volatile basic nitrogen during fishing and aquaculture stages to ensure the content does not exceed the standard. At the processing stage, optimizing processing techniques (such as pre-cooling treatment[49,50] and mild processing methods[51]) should be implemented to shorten processing time. Technologies such as modified atmosphere packaging[52] and biological preservation [53] can be adopted to delay spoilage. The logistics end can enhance the cold chain delivery system and equip it with temperature and humidity monitoring devices to ensure temperature control throughout the entire process from catch to retail.

Sulfur dioxide warrants particular attention among food additives. It effectively inhibits microbial growth and delays food spoilage. Some enterprises, to improve the color of aquatic products and extend the shelf life, use excessive amount of additives for bleaching or preservation, which may lead to residues that pose potential health risks. Long-term excessive intake can cause damage to the digestive system and liver function[54]. Therefore, it is necessary to strictly control the amount and residue standards of such additives in aquatic products. Through the synergistic effect of technological alternatives, regulatory reinforcement, and consumer education, the risks of sulfur dioxide and other additives can be systematically reduced, thus ensuring the safety of aquatic products and the sustainable development of the industry.

4.2. Risk Variations Across Different Stages and Their Spatiotemporal Distribution Patterns

Risks associated with aquatic products differ between the distribution stage and catering stage, with distribution posing higher risks than catering. In the distribution stage, risks are primarily concentrated in farmers' markets and supermarkets, where common issues include “failure of cold chain management” and “contamination from multiple-party handling”. Aquatic products requires low temperatures throughout the entire process. However, in farmers’ market, temperature fluctuations often occur in the seafood section due to untimely replenishment of ice packs and frequent opening and closing of refrigerated display cases. Moreover, supermarkets face the risk of damaged packaging during multi-tiered transportation from “production site→distribution center→storefront”. Both scenarios accelerate microbial proliferation and protein degradation, leading to exceedance indicators such as volatile basic nitrogen and total bacterial count. Risks in the catering stage are primarily concentrated in medium-sized and small restaurants. The core reason for the overall low risk in this stage is the immediate sensory feedback from consumers. The freshness of aquatic products can be quickly assessed through appearance, smell, and texture. Consumers reject any abnormal products they detect, which pressures food service businesses to promptly remove non-compliant items. However, risks in the catering stage remain significant and should not be overlooked. Regulatory authorities should focus on key sampling locations such as farmers' markets, supermarkets, and small-to-medium-sized restaurants to reduce the risks associated with aquatic products.

In terms of time, aquatic product risks exhibit significant seasonal variations, with the risk values from August to November being significantly higher than those from December to April. The high temperatures and humidity in summer provide optimal conditions for the proliferation of microorganisms, accelerating protein degradation. Conversely, the low temperatures in winter inhibit microbial activity and slow down the migration of contaminants, thereby reducing risk. Regulatory authorities should optimize sampling strategies for the summer and autumn, enhance the coverage of cold chain transportation during summer, and promote technologies such as chilled packaging to extend the shelf life of aquatic products. In spatial distribution, the risks exhibit a pattern of being higher in eastern regions than western regions and higher in coastal areas than inland areas. Zhejiang, Chongqing, and Jiangxi are classified as high-risk provinces. These areas are either major aquaculture production zones (characterized by high farming density and frequent disease outbreaks) or distribution hubs (involving multiple distribution stages and significant cold chain pressure). Strengthened sampling and inspection efforts should be implemented in high-risk provinces. The distribution of non-compliant heavy metals indicates that northern regions and coastal provinces should be given particular attention. The distribution of non-compliant quality indicators and additives shows scattered occurrences without discernible spatial patterns. Regarding non-compliant veterinary drugs and prohibited substances, coastal and southwestern regions require heightened focus and strengthened supervision and sampling to reduce risks. The global spatial autocorrelation analysis reveals that only cadmium, enrofloxacin, furazolidone metabolites, and chloramphenicol have clustered distribution patterns, while other risk substances show random distribution. The local spatial autocorrelation analysis suggests that inspection efforts for cadmium should be strengthened in Beijing, Tianjin, Hebei, Inner Mongolia, Liaoning, Jilin, Heilongjiang, Shandong, and their surrounding provinces. It is also recommended to strengthen the sampling intensity of enrofloxacin in Shanghai, Jiangsu, Zhejiang, Anhui, Fujian, Jiangxi, Hubei, and Hunan, as well as their neighboring provinces. Sampling intensity for furazolidone metabolites should be diminished in Hebei, Inner Mongolia, Jilin, and their neighboring provinces. While we should increase the sampling intensity for chloramphenicol in Guangdong, Hainan, and their neighboring provinces. These characteristics provide precise targets for “regional, pollutant-specific regulation”.

5. Conclusions

Various hazards are present throughout the entire supply chain of aquatic products, spanning aquaculture, capture, transportation, processing, and food service. However, the risk levels associated with different risk substances vary significantly across stages and spatiotemporal dimensions. This study has clarified the rankings of risk substances in Chinese aquatic products along with their spatiotemporal evolution patterns. While regional variations exist in the risks, the overall situation remains manageable. The top three risk substances in aquatic products are cadmium, enrofloxacin, and total volatile basic nitrogen. Among the high-risk and higher-risk substances in aquatic products across Chinese provinces, cadmium, enrofloxacin, furazolidone metabolites, and chloramphenicol have significant global spatial autocorrelation in non-compliance. Non-compliant aquatic products primarily occur in urban areas in China, with the distribution stage identified as the main process and the primary venues being farmers' markets and supermarkets. Among the sampled provinces, Zhejiang exhibits the highest non-compliance rate. Spatiotemporal analysis found that aquatic products pose higher risks during summer and autumn, while risks are lower in winter and spring. Based on the risk classification and spatiotemporal analysis, future monitoring strategies should be optimized to implement targeted and differentiated inspection of aquatic products, and thus enhance the efficiency of regulatory oversight and mitigate quality and safety risks associated with aquatic products. Additionally, monitoring frequency should be intensified in key high-risk provinces such as Zhejiang, particularly during summer and autumn and prioritizing critical venues including supermarkets and farmers' markets. For major risk substances such as cadmium and enrofloxacin, a comprehensive traceability system covering the entire supply chain—from aquaculture sources to point of sale—should be established. Strengthening cross-departmental collaborative oversight will facilitate the efficiency of the risk warning and response mechanisms. To strengthen comprehensive control across the entire global supply chain, customs authorities should implement differentiated sampling inspections for high-risk product categories and origins, leveraging digital technologies to enhance the precision of port inspections. Furthermore, there is a need to promote the establishment of unified cross-border traceability standards to ensure information traceability throughout the product circulation process from exporting to importing countries. We should, by deepening international regulatory collaboration, construct a transnational risk information-sharing mechanism to promptly exchange global updates on detected risk substances in aquatic products and disease outbreak information. By jointly conducting source risk assessments, it will reduce safety risks associated with aquatic product .

This study conducted risk classification based on the supervision data of aquatic products in China from 2021 to 2023. However, there lacks data from production stage in the findings, which prevents the assessment of risk classification for substances in that stage. Furthermore, the study failed to differentiate risk variations among the subcategories of aquatic products and examine the synergistic effects between different contaminants. Future studies should integrate data across the entire aquatic product supply chain, from aquaculture to consumption.Extending the data collection period, incorporating risk analysis for various aquatic product types, and exploring synergistic mechanisms among different risk substances will better guide the prevention and control of seafood risks.