Introduction

Asia is the main market of dairy products, accounting for about 55% of global imports (Agrimonti et al., 2019). Milk production in China is 36.55 million tons, accounting for about 4.5% of global production after India and the United States, and ranks third in the world (Xu et al., 2022). Dairy products are produced by heating, drying out, and fermentation of milk (Hu et al., 2022). These products are rich in nutrients (proteins, lipids, carbohydrates, essential minerals, and vitamins that are easily absorbed by the body) and are particularly worthy for children, convalescents, and the elderly (Nastasescu et al., 2020). However, the dairy and milk industry is among the most polluting fields of production (Feil et al., 2020); hence, quality and safety of dairy products must be given utmost attention.

The main pollutants of dairy products are microorganisms, pesticides, veterinary drugs, and inorganic pollutants. Microorganisms in dairy products are mainly the colonies of coliforms (Slosarkova et al., 2021), mold, and yeast (Stobnicka-Kupiec et al., 2019). Microorganisms are derived from contaminated colostrums, which enter raw milk from beddings in cowsheds, feeds, animal skin, hair, and milking equipment as well as through workers’ clothes and hands (Stobnicka-Kupiec et al., 2019). Lead (Pb), chromium (Cr), cadmium (Cd), mercury (Hg), and arsenic (As) are inorganic pollutants that contaminate milk and dairy products (Yan et al., 2022). The oxidation state of Cr results in health risks to humans exposed to it, with hexavalent Cr (Cr[VI]) being the most hazardous constituent. Exposure to Cr(VI) causes the internal reduction of chromates to oxidize guanine residues, releasing free radicals and reactive oxygen species (ROS) in the process. This has genotoxic effects, because the reactive forms of Cr are able to damage deoxyribonucleic acid (DNA) bases (Iyer et al., 2023).

The principal sources of these inorganic pollutants in dairy products are connected to environmental conditions, manufacturing methods, and feed and soil quality, which are transferred through nutrient chains (i.e., soil–plant, cattle, milk, dairy–human) (Caggianob et al., 2006; Li et al., 2019; Taghizadeh et al., 2018). In addition to the above-mentioned contaminants, aflatoxin M1 (AFM1) is an important risk factor in the safety of milk and dairy products because of its carcinogenic and teratogenic effects in humans (Chen et al., 2023). When mammals, including cows, consume contaminated feed, their milk is contaminated with AFM1 (Nishimwe et al., 2022).

European monitoring of dairy products focuses on relevant hazards of supply chains with the implementation of good agricultural practices (GAP) and hazard analysis critical control points (HACCP) at all stages of supply chains (van Asselt et al., 2017). Food safety monitoring is performed by sampling feeds from farms, dairy farms and their storage silos, and transportation trucks, considering them as possible control targets (Wang et al., 2020, 2021b).

New Zealand manages dairy production process mainly by including management at dairy breeding stage, raw milk production, and dairy processing, and selects corresponding key targets for supervision and management (Wang et al., 2018). In addition, veterinary drugs are fully monitored in the United States, Israel, and Europe (Bommuraj et al., 2020). Dairy supervision in China involves raw milk production and acquisition (Zuo et al., 2020); however, problems still dwell in the areas of regulatory measures (Chen et al., 2020). Difference between time and space is not considered during the supervision of dairy products. Therefore, this paper uses the entropy weight method (EWM) to calculate the risk of each control point and discusses the time and space distribution characteristics of Chinese dairy products according to the involved risk.

This study aims to: (1) classify potential contaminants in dairy products using EWM and (2) analyze the temporal and spatial patterns of potential contaminants in dairy products. The results would provide a reference for supervising and quality control of dairy products in China.

Research methods and data sources

The following four risk indicators were selected for analyzing dairy products: failure rate, detection rate, qualification degree, and hazard degree. Weight of these four indicators was calculated by the entropy value method. The risk value of each risk indicator was comprehensively calculated, and graded by the Pareto Principle (Benjamin Harvey and Sotardi, 2018; Zou et al., 2018). The temporal and spatial distribution characteristics of each risk indicator thus derived were discussed.

Data sources

The data for this study were derived from the 2017–2021 dairy product supervision and sampling data of the China Market Regulatory Bureau and the Market Supervision Administration of various provinces and municipalities that function directly under the Central Government and autonomous regions (excluding Hong Kong, Macao, and Taiwan).

Entropy weight method

Entropy weight method, originally a thermodynamic concept, was later introduced into information theory to unify and reflect changes in the information provided by each index, allowing for the avoidance of interference from human factors in the weighting of each evaluation index and producing more objective evaluation results (Zhao et al., 2022). The lower the entropy value, the higher the weight, and the more information the indicator provided (Cunha-Zeri et al., 2022). Matrix X obtained from m risk substances and n indicators is as follows:

The original matrix was normalized, and the range of indicators was based on [0-1] (Equations 2 and 3) (Zhao et al., 2022):

A positive index means the higher the index value, the better the index. A negative index means the lower the index value, the better the index. Each indicator was measured together and the proportion of hazards in the index (P_ij) under an indicator was calculated using Equation 4:

Then, the entropy of the quality index was calculated by Equations 5 and 6:

The weight of each index was calculated according to entropy, and the weight of each quality index was expressed by Equations 7 and 8:

The risk value of each risk substance was calculated using Equation 9:

Spatial autocorrelation analysis

Both correlation of the same variable at various spatial locations and association between a variable's value at one location and values close at hand are examples of spatial autocorrelation. Both geographical clustering and geographical outliers are a part of spatial autocorrelation. A high-concentration variable is clustered spatially if it is surrounded by other high-concentration variables. In contrast, an outlier is clustered spatially if a variable’s high concentration is surrounded by its low concentration. Spatial autocorrelation develops between variables if they display some degree of spatial regularity. Both global and local markers are used in spatial autocorrelation. Global spatial autocorrelation uses a single number to represent the level of autocorrelation in that region and represents the overall distribution of a variable to see whether spatial clusters of that variable exist within a broader region. The local spatial autocorrelation method determines how each spatial cluster is in relation to its surrounding clusters on an attribute (Liu et al., 2013).

Global spatial autocorrelation

Spatial analysis is used to identify clustered areas and observe geographic variations (Zuo et al., 2020). The spatial distribution characteristics of the qualified rate of risk substances with high-risk grade and high nonconforming conditions in dairy products were calculated. Moran's I value was between -1 and 1, where values close to 1 represented positive spatial autocorrelation, values close to -1 represented negative spatial autocorrelation, and 0 is random distribution. The local Moran index was calculated and hot spot analysis was performed to determine location of the cluster. It is usual to associate weight (w_ij) to each pair (x_i, x_j), which quantifies spatial pattern (Asosega et al., 2021) as shown in Equation 10:

where n is the data from Chinese provinces (excluding Hong Kong, Macau, and Taiwan); x_i and x_j represent the pass rate of risk substances in China, and i and j are the average pass rate and the spatial weight matrix, respectively.

Local spatial autocorrelation

Using the local Moran index to determine spatial autocorrelation and detect some spatial clusters with similar adjacent features and outliers, local indicators of spatial association (LISA) was used to study the spatial correlation of attribute values of spatial units in a local range (Wang et al., 2021a), which is shown in Equation 11:

where S² is the variance of x_i or x_j, and the rest is interpreted in the same way as the global Moran's I.

Statistical analysis

Results for each hazardous material were calculated as mean (x̄). Significance analysis was carried out using SPSS (version 23). ArcGIS 10.8 was used to create a map of spatial analysis, and Origin 2021 was used to construct additional figures.

Results and discussion

Precipitation of hazards in dairy products

After processing the sampling data of dairy products from 2017 to 2021, 348,069 valid product samples were included for analysis. Hazard detection indicators in the national sampling data include AFM1, dexamethasone, acesulfame potassium, ampicillin, benzylpenicillin, o-clopenicillin, sucralose, sorbic acid, zearagiferol, β-lactamase, aspartame, ethyl carbamate, benzoic acid, defloxacin, coliform bacteria, listeria monocytogene, polychlorinated biphenyls (PCBs), enrofloxacin, chromium, yeast, the total number of colonies, chloramphenicol, mold, natamycin, norfloxacin, lead, melamine, and 34 other indicators: These indicators were divided into four categories: pesticides, veterinary drugs, additives, and microorganisms. Detection and failure rates, as well as the qualification and hazard degrees of each item in different links, were calculated. In addition, the comprehensive risk value of each item was calculated according to the weight, and the hazards were graded according to the ranges summarized in Table 1 (Benjamin Harvey and Sotardi, 2018; Zou et al., 2018). The results showed (Table 2) that high and higher-risk substances in circulation and production were mainly veterinary drugs, additives, microorganisms, and chromium. Veterinary drugs included dexamethasone, ampicillin, benzylpenicillin, enrofloxacin, etc. Additives included saccharin sodium. Microorganisms were mainly coliforms, total bacterial count, molds, and yeasts. High and higher-risk substances in the catering sector were sucralose, coliforms, chromium, total bacterial count, mold, and yeast. In general, high-risk substances in dairy products were mainly microorganisms, veterinary drugs, additives, and heavy metals.

Unqualified items in the national sampling project mainly included five categories: microorganisms, additives, veterinary drugs, plasticizers, and physical and chemical indicators. Among these, microorganisms, additives, veterinary drugs, and plasticizers were quality and safety indicators, which could harm the human body. From 2017 to 2021, the unqualified items of dairy products were microorganisms (47.2%), physical and chemical indicators (28.3%), additives (19.7%), veterinary drugs (3.9%), and plasticizers (0.8%) (Figure 1), indicating that the quality and safety of dairy products are mainly dictated by microorganisms and additives. Microorganisms included coliform bacteria, yeast, total number of colonies, and mold. The main problem was yeast (43.3%), followed by coliform bacteria (21.7%), mold (21.7 %), and, finally, the total number of colonies (13.3 %) (Figure 2b). The additives were mainly benzoic acid and natamycin, which accounted for 63.6% and 42.4% of unqualified products, respectively (Figure 2c). The main unqualified veterinary drug was chloramphenicol (Figure 2d). These results showed that more attention should be paid to microorganisms and food additives in dairy products, and this intended that control of physical and chemical indicators, such as proteins (27.8%), acidity (38.9%), fat (19.4%), and non-fat milk solids (13.9%), in dairy products should be strengthened (Figure 2a).

Preservatives prolong the shelf life of dairy products; commonly used preservatives in dairy products are benzoates, sorbates, and natamycin, which inhibit growth and reproduction of microorganisms (Molognoni et al., 2016). Benzoic acid and sorbic acid have high solubility and high bactericidal effect, while natamycin has low solubility and no bactericidal effect (Tfouni and Toledo, 2002).

Iran prohibits the use of additives in dairy products, and the maximum content of sorbate additives allowed in processed cheese is 1,000 mg/kg (Zamani Mazdeh et al., 2017). However, a small quantity of preservatives and sweeteners is not harmful to the human body, but excessive amounts may harm human organs (Li et al., 2016b). Dairy products are a common source of pesticide contamination and are considered a global health problem. Pesticide residues in dairy products originate from animal feed, environmental contamination (air, soil, and water), and veterinary products formulated with active pesticide ingredients (Schopf et al., 2022). Heavy metal contamination in dairy products is one of the main public health issues. Cadmium is an issue because of its adverse effects on humans; it contaminates dairy and milk products if mammals graze on contaminated land; Cd contamination is also induced from water, equipment, and utensils used in the processing of dairy products (Abd El-Salam et al., 2017).

Yan et al. (2022) studied the content of inorganic pollutants in dairy products in China, and discovered toxic heavy metals, including Hg, Pb, Cd, Cr, and As, in milk and dairy product samples, with Pb levels exceeding European Union (EU) standards (0.02 mg/kg). In China, the maximum residue limits (MRLs) of Pb, Cr, Hg, and As are respectively set as 0.05, 0.3, 0.01, and 0.10 mg/kg. Pb and Hg have detection rates of 29.2% and 28.1%, and concentration range of 2.01–33.37 µg/L and 1.02–7.57 µg/L, respectively. Cr and As have a concentration range of 8.01–84.66 µg/L and 1.03–2.30 µg/L, and detection rate of 12.4 and 9.0%, respectively. Cd was not found in certain samples (Qu et al., 2018). In a study conducted by Su et al. (2021), milk from industrial areas had Cr and As concentrations of 2.41±2.12 µg/kg and 0.44±0.31 µg/kg, respectively, which were considerably higher than the milk obtained from nonindustrial areas (P < 0.01). The content of Pb and Cd in Egyptian dairy products was reported to be higher (Abd El-Salam et al., 2017). A study conducted in Bangladesh found dairy products containing higher levels of Pb but within controllable limits (Aktar et al., 2020). Dairy products are an ideal medium for microbes and are highly susceptible to contamination (Cai and Zhao, 2017).

Temporal distribution of dairy hazards

Data analysis of sampled dairy products from 2017 to 2021 showed that the nonconformity of dairy products was relatively stable (Table 3). The nonconformity rate fluctuated less, except in 2019 and 2021, when it was slightly above average. The nonconformity rate increased in 2021, indicating that the nonconformity of dairy products went up in recent years; hence, it was necessary to strengthen the management and monitoring of dairy products.

Microorganism

Temporal analysis of microbial occurrence in dairy products (Figure 3a) for 2017–2021 showed microbial nonconformities for each year. The highest occurrence of nonconformities was for yeast (0.973%), followed by the total number of molds and colonies for 2017–2018. All microbial disagreement rates decreased but increased sequentially in 2019 by 0.611%, 0.302%, 0.203%, and 0.169%. Since 2019, the failure rate has dropped significantly.

Presence of coliform bacteria is the initial indicator of fecal contamination of water bodies and is used as an indicator of unsanitary conditions in dairy products and other foods (Martin et al., 2016). These bacteria interfere with immunoglobulin absorption in a variety of ways (Cummins et al., 2017). Many factors affect the total number of colonies in raw milk. Major issues are the milking process, hygiene of milking equipment, disinfection, and milk storage conditions (Elmoslemany et al., 2009).

Studies have shown that the total number of colonies in raw milk is also related to cowshed environment because it demonstrates infection in udder of cows (Van Kessel et al., 2004). Fungi and yeast are the main contaminating bacteria in fermented milk. These can even survive in low pH environment. Fruits and grains added in fermented milk are other major sources of yeast and fungi. The environment and equipment used during processing are also one of the main causes of pollution (Zhang et al., 2019).

Molds produce mycotoxins that are toxic to both humans and animals. Grains, feed, and straws are susceptible to mycotoxin contamination (Akinyemi et al., 2021). Studies conducted in Mexico demonstrated the presence of AFM1 even in milk treated at ultra-high temperatures. All samples of feed, raw milk, and pasteurized milk showed AFM1 contamination (26.0 ± 0.4 µg/kg, 32.0 ± 1.0 ng/L, and 31.3 ± 0.7 ng/L, respectively), and a significant proportion (90.4%, 11.3%, and 10.3%, respectively) exceeded the locally applied maximum permissible limits for feed and milk (20.0 µg/kg and 50 ng/L, respectively; Álvarez-Días et al., 2022).

Veterinary drug residues

Analysis of the failure rate of veterinary drugs is shown in Figure 3b, which indicates that chloramphenicol failure mainly occurred in 2017–2018. The highest failure rate was noted in 2018 at 0.296%, but no failure occurred in the next three years. Antibiotics are used in the dairy industry for prevention and treatment of diseases in animals, including treatment of mastitis in dairy cows (Rajala-Schultz et al., 2021). However, antibiotic resistance in animals has become one of the biggest common health issues. Chloramphenicol is the most commonly used antibiotic in veterinary. Antibiotics are injected in animals and excreted into the environment through feces and urine. Antibiotics that are still in the withdrawal stage may also be excreted through cow milk. Antibiotic residues in the environment or in infected milk lead to the contamination of dairy products (Jindal et al., 2021).

Additives

Additives are mainly benzoic acid, sorbic acid, dehydroacetic acid, and natamycin (Figure 3c). Failure of benzoic acid was observed mainly in 2018 and 2019, with the highest rate being 0.428% in 2018; failure of sorbic acid and dehydroacetic acid occurred in 2019 (0.052%) and 2021 (1.786%), respectively. The maximum failure of natamycin was 0.766% in 2019 and this occurred in other years as well, except in 2018. Regulations for using dehydroacetic acid and natamycin must be strengthened. Additives are added to dairy products for preservation and sterilization effects, but use of too many compounds could be hazardous to the human health.

Plasticizers

Dibutyl phthalate (DBP) demonstrated nonconformity in 2019 with a failure rate of 14.286% (Figure 3d). Phthalates are a class of chemicals that are frequently found in polyvinyl chloride plastics, food packaging, and cosmetics. They mostly leave behind DBP and dimethyl phthalate (DMP). More than 100 g/kg of DBP was discovered in milk samples analyzed by Li et al. (2011). Children exposed to diethylhexyl phthalate (DEHP) and benzyl butyl phthalate (BBP) had an increased risk of allergic conditions, including asthma and eczema as well as a possible reduction in anal-genital distance in fetuses and male infants (Ge et al., 2016). DBP can impede with hormones, which might cause cancer risk and developmental flaws (Gurudatt et al., 2022). Plasticizer contamination occurs through packaging materials based in milking machinery. When used as packaging, polymer materials containing these chemicals exude from the materials and released into food (Ge et al., 2016). Therefore, despite the fact that plasticizers deviated from the norm only in 2019, these still require monitoring.

Spatial cluster distribution and geographic features

From 2017 to 2021, the highest failure rate was 1.491% in Ningxia (in 2019). Provinces with the overall higher failure rates in 2017–2019 were located in northwest China (Tibet, Xinjiang, Qinghai) and northern Inner Mongolia and Ningxia. Provinces with high failure rates in 2020 and 2021 gradually moved southward. The highest failure rate in 2020 was noted in Hainan province (0.175%), and in 2021, failure was mainly distributed in Qinghai (0.605%) and Jiangxi (0.568%) provinces.

Popular areas of the Chinese dairy industry are northeast Inner Mongolia and north and south China (Cheng and Pang, 2016). Based on the data for the western provinces of China from 2013 to 2017, Zhou et al. (2021) analyzed the supervision level of dairy products and discovered that supervision level in Chongqing, Sichuan, Guizhou, Yunnan, Shaanxi, and Xinjiang was low, reflecting negativity in the stated region. The unqualification rates of dairy products in this study were mainly concentrated in the western region of China, which is the same as the results of the study conducted by Zhou et al. (2021) (Figure 4). Therefore, it is necessary to focus on the quality and safety of dairy products in Inner Mongolia, Xinjiang, Tibet, Qinghai, Sichuan, Ningxia, Jiangxi, Yunnan, and Hainan.

Spatial autocorrelation analysis

According to the grading results shown in Table 2, high and higher-risk hazards were selected as research objects, including coliform bacteria, yeast, mold, and total number of colonies. The spatial distribution of the qualification rates of these hazards given in Figure 5 shows that although there was no global spatial autocorrelation (P > 0.05) for coliform flora, mold, total number of colonies, and yeast, there could be a local spatial correlation.

The local spatial correlation analysis of coliform flora, mold, total number of colonies, and yeast is shown in Figure 6. Only coliform flora and total number of colonies had high–high aggregation (Figures 6a and 6b), with Shanghai, Anhui, and Jiangsu presenting high qualification rates; the surrounding qualification rates were low and mainly concentrated in the eastern region of China. Only total number of colonies had low–low (L-L) aggregation (Figure 6c), and the total number of colonies adjacent to it was lower than the national qualification rate. Qualification rate presented a low–high correlation for coliform bacteria in Xinjiang, Liaoning, and Anhui provinces; for mold in Zhejiang province; for total number of colonies in Xinjiang province; and for yeast in Xinjiang and Hunan provinces. Their qualification rate was low, although it was relatively high for neighboring provinces. The degree of difference in the pass rate across neighboring provinces was high and belonged to the high–low aggregation correlation. The high pass rate-related areas were mainly distributed in Hainan, Gansu, and Shanxi provinces, but a low pass rate was observed in their neighboring provinces as well.

Since there is no procedure to control the growth of microorganisms, the poor microbiological pass rate of dairy products in northwest China could be related to the preference of nomadic herders for manufacturing spontaneously fermented dairy products (Guo et al., 2021). Knowledge, attitudes, and behaviors of dairy workers in Beijing were investigated by Chen et al. (2018). Only 39% of 194 dairy workers were able to identify specific pathogens (particular bacteria) or diseases that could be spread through milk and dairy products, and only 24% of employees were aware of proper handwashing techniques. The study’s findings showed that dairy workers in northern China had comparatively little awareness about food safety. Hence, this could lead to the microbiological contamination of dairy products.

Coliform bacteria, mold, total number of colonies, and yeast showed random distribution in China. It was recommended that the intensity of sampling for coliform bacteria must be strengthened in Xinjiang, Liaoning, and Anhui provinces. It was further urged to strengthen sampling for mold in Zhejiang province, for total number of colonies in the surrounding provinces of Xinjiang, Gansu, and Ningxia, and for yeast in Xinjiang and Hunan provinces.

The results of risk classification showed that the main high-risk hazards in dairy products are microorganisms. It was recommended to strengthen the sampling and control of microorganisms. Microbial contamination in dairy products mainly occurs through raw materials, production links, and workers. Therefore, it is necessary to strengthen microbial control during the production processing of raw milk and dairy products, which could be possible by thorough cleaning of equipment during production and by strengthening the health management of staff and sites. In addition, it is necessary to prevent collision and extrusion and reduce mechanical damage to packaging material during storage and transportation (Li et al., 2016a). In addition to the managing of equipment and personnel, it is recommended to use heat and high-pressure treatments and rational use of chemical additives to destroy bacteria or inhibit microbial growth. The packaging must be carried out in hygienic and sterile conditions to reduce the risk of biological pollution (Garnier et al., 2017).

Conclusion

Based on the sampling data of dairy products in China from 2017 to 2021, the present study presented risk classification and spatiotemporal distribution characteristics of hazardous substances in dairy products. The results showed that the main high and higher-risk substances in circulation and production were veterinary drugs, additives, microorganisms, and heavy metal chromium, while the main high and higher-risk substances in catering were additives, microorganisms, and chromium. Main unqualification items in dairy products were microorganisms, additives, veterinary drugs, plasticizers, and physical and chemical indicators of five categories.

From the perspective of time distribution characteristics for 2017–2021, it was observed that unqualified microorganisms occurred every year. The highest unqualification rate was presented by yeast (0.973%), followed by mold and total bacterial count. The unqualification rate of chloramphenicol occurred in 2017–2018, with the highest unqualification rate being 0.296% in 2018. The unqualification rate of benzoic acid was highest in 2018 and 2019, with the highest being 0.428% in 2018. The respective highest unqualification rate of sorbic acid and dehydroacetic acid was in 2019 (0.052%) and 2021 (1.786%). DBP showed nonconformity in 2019, with an unqualification rate of 14.286%.

According to the spatial distribution characteristics for 2017–2021, the highest unqualification rate was 1.491% in Ningxia (in 2019). Provinces having higher unqualification rates of dairy products during 2017–2019 were located in the northwest of China (Tibet, Xinjiang, Qinghai) and northern Inner Mongolia and Ningxia. From 2020 to 2021, provinces with high unqualification rates gradually moved southward. The highest unqualification rate was recorded for Hainan province (0.175%) in 2020. In 2021, it was mainly distributed in Qinghai (0.605%) and Jiangxi provinces (0.568%).

Microorganisms presented the highest of all risks. Therefore, appropriate measures must be taken to control them via sanitary management of equipment and personnel during production as well as thermal and high-pressure treatment of dairy products or through the rational use of chemical additives. In addition, it is necessary to strengthen the sampling inspection of dairy products in the western region. Through the classification of risk substances and the characteristics of related time and space distribution, the dynamic adjustment of targeted regulatory opinions and regulatory plans was proposed to provide scientific basis for the safe supervision of dairy products.

This study was based on the data for risk classification collected through the monitoring and sampling program of Chinese dairy products from 2017 to 2021. The study’s findings excluded information on raw material links, making it impossible to estimate risk level of each raw material and classify dangerous compounds. In order to undertake full analysis, subsequent study should incorporate involvement between the supervision and sampling data and agricultural productivity.

Author Contributions

Guangcan Tao: Funding acquisition, visualization, and supervision. Shengmei Liao and Kang Hu: Writing, including original draft preparation, software, formal analysis, and investigation. Qianghai Zhang: Project administration. Su Xu: Resources. All authors read and agreed to the final version of the manuscript.