*Article* **Spatial-Temporal Characteristics in Grain Production and Its Influencing Factors in the Huang-Huai-Hai Plain from 1995 to 2018**

#### **Chunshan Zhou <sup>1</sup> , Rongrong Zhang <sup>1</sup> , Xiaoju Ning 2,\* and Zhicheng Zheng <sup>3</sup>**


Received: 30 October 2020; Accepted: 7 December 2020; Published: 9 December 2020

**Abstract:** The Huang-Huai-Hai Plain is the major crop-producing region in China. Based on the climate and socio-economic data from 1995 to 2018, we analyzed the spatial–temporal characteristics in grain production and its influencing factors by using exploratory spatial data analysis, a gravity center model, a spatial panel data model, and a geographically weighted regression model. The results indicated the following: (1) The grain production of eastern and southern areas was higher, while that of western and northern areas was lower; (2) The grain production center in the Huang-Huai-Hai Plain shifted from the southeast to northwest in Tai'an, and was distributed stably at the border between Jining and Tai'an; (3) The global spatial autocorrelation experienced a changing process of "decline–growth–decline", and the area of hot and cold spots was gradually reduced and stabilized, which indicated that the polarization of grain production in local areas gradually weakened and the spatial difference gradually decreased in the Huang-Huai-Hai Plain; (4) The impact of socio-economic factors has been continuously enhanced while the role of climate factors in grain production has been gradually weakened. The ratio of the effective irrigated area, the amount of fertilizer applied per unit sown area, and the average per capita annual income of rural residents were conducive to the increase in grain production in the Huang-Huai-Hai Plain; however, the effect of the annual precipitation on grain production has become weaker. More importantly, the association between the three factors and grain production was found to be spatially heterogeneous at the local geographic level.

**Keywords:** grain production; spatial–temporal characteristic; influencing factors; the Huang-Huai-Hai Plain

### **1. Introduction**

China is an important food-producing country in the world, as well as a large food consumer. China's food self-sufficiency rate has reached more than 95% [1]. Although the current grain supply and demand in China maintains a balance of total quantity and a surplus in harvest, the small per capita arable land area, low mechanization level, and family-based business units determine that the current arable land has an extremely limited potential for increasing grain production [2,3]. From 1850–1900 to 2006–2015, the mean land surface air temperature increased by 1.53 ◦C [4]. Climate change has already affected food security due to warming, changing precipitation patterns, and the greater frequency of some extreme events [4]. Under the dual constraints of climate change and the process of urbanization, unpredictable meteorological disasters, limited arable land resources, huge population pressure, and the diversified consumer demand of residents directly generated strong demands for stable grain production [5,6]. Additionally, the outbreak of COVID-19 in 2019 led to labor shortages

and supply chain disruptions, which affected the food security of some countries and regions [7]. Certain countries have even banned the export of food, leading to large fluctuations in global food prices and casting a shadow over the world's food crisis. Although China's major food crops, such as rice and wheat, are less dependent on the international market and national food security will not be affected by the agricultural trade restrictions brought about by the COVID-19, ensuring China's food security in the later period of the epidemic is of great strategic significance for domestic economic recovery and social stability.

The Huang-Huai-Hai Plain (HHH), located in the north of the country, is a high-yield agricultural region [8], accounting for 70% of wheat and 30% of maize production in China [9]. Because of its importance in grain production and China's self-sufficiency in food, it is known as the 'Breadbasket of China' [10,11]. Changes in grain production in the HHH can have direct impacts on both the national economy and food security of China [12]. Therefore, it is of great significance to understand its temporal and spatial characteristics and the influencing factors of grain production in the Huang-Huai-Hai plain in order to ensure national food security.

Previous studies have pointed out that while grain production in the Huang-Huai-Hai Plain is increasing, the difference in grain yields in various regions is gradually shrinking, but the spatial distribution pattern of grain production in the southern region of the Huang-Huai-Hai Plain is higher than that in the northern region, which remains unchanged [13–15]. Later studies attribute the distribution characteristics of grain production in the HHH to socio-economic factors [13,14,16]. The increase in grain production in the HHH is mainly due to the improvements of crop varieties, fertilizers and effective irrigation area. For example, the amount of fertilizer used has increased by about 400% and the effective irrigation area has increased by about 20% [10]. Liu, Tang [15] used the spatial lag model to reveal the factors affecting the differentiation in grain yield in the HHH from 1995 to 2010, and found that farmers' per capita net income, effective irrigation area ratio, and industrial structure had significant positive effects on grain yield. At the same time, studies have also found that although grain production and fertilizer input in the HHH are still significantly positively correlated, the current application of fertilizers in actual production is extremely unbalanced and unreasonable [17]. It is essentially important to improve fertilizer use efficiency in order to save resources as well as increase yield [10].

A number of studies have also emphasized the role of climate factors, especially precipitation, in grain production in the HHH [18–22]. Qu, Li [20] indicated that the increase in precipitation in the HHH can significantly increase the food production in the HHH, but the increase in thermal resources will increase the shortage of water resources and offset the impact of the increase in temperature. Xiao, Qi [23] reveled that climate change reduced the potential winter wheat yield of 80% of the stations by 2.3–58.8 kg·yr−<sup>1</sup> , while at the same time it is pointed out that increasing the heat time of the wheat growth period is essential to alleviate the impact of the shortening of the growth period caused by warming climatic conditions. However, if the advancement of agricultural technology and other non-climatic factors are taken into account, for every 1 ◦C increase in the average temperature of the Huang-Huai-Hai Plain, the winter wheat yield in the north will increase by 2.1%, and the yield in the south will decrease by 4.0% [22].

Therefore, most of the previous research on grain production in the Huang-Huai-Hai Plain mainly focuses on the factors affecting grain production, that is, mainly from two aspects: socio-economic factors and climatic factors. Further, these studies are only carried out on one level, which separates the comprehensive impacts of climate and socio-economic factors on grain production. This will inevitably affect the final assessment results. To complete these data, this paper aims to explore the influencing factors of grain production in HHH by incorporating the climate factors and socio-economic factors into the models so as to provide a reference for ensuring food security and relevant departments to make decisions.

#### **2. Study Area and Data Sources 2. Study Area and Data Sources**

#### *2.1. Study Area 2.1. Study Area*

The Huang-Huai-Hai Plain (HHH), the second largest plain in China, is located at 32◦–40◦ N and <sup>114</sup>◦–121◦ E, with a land area of about 4 <sup>×</sup> <sup>10</sup><sup>5</sup> km<sup>2</sup> , spanning seven provinces and cities of Beijing, Tianjin, Hebei, Shandong, Henan, Anhui, and Jiangsu (Figure 1). The HHH belongs to the warm temperate semi-humid monsoon climate zone and is one of the most sensitive areas to climate change in China. The winter climate in this area is typically dry and cold, spring is dry, with less rain and much evaporation, and summer is characterized by high temperatures and heavy rainfall, including high intense rainfall that often leads to summer floods [24]. As an important agricultural region, this area has a long history of farming and is an important grain production base for food security in China, with its sown area of 20.4% of the nation's farmland and 23.6% of the whole nation's grain yield [25]. The Huang-Huai-Hai Plain (HHH), the second largest plain in China, is located at 32°–40° N and 114°–121° E, with a land area of about 4 × 105 km2, spanning seven provinces and cities of Beijing, Tianjin, Hebei, Shandong, Henan, Anhui, and Jiangsu (Figure 1). The HHH belongs to the warm temperate semi-humid monsoon climate zone and is one of the most sensitive areas to climate change in China. The winter climate in this area is typically dry and cold, spring is dry, with less rain and much evaporation, and summer is characterized by high temperatures and heavy rainfall, including high intense rainfall that often leads to summer floods [24]. As an important agricultural region, this area has a long history of farming and is an important grain production base for food security in China, with its sown area of 20.4% of the nation's farmland and 23.6% of the whole nation's grain yield [25].

Here, wheat and maize occupy a larger proportion in the structure of grain production. The annual double-cropping system of winter wheat and summer maize is the most popular planting pattern. In addition, the yield of wheat and maize accounts for approximately 61% and 31% of the total national output, respectively [26]. Therefore, it is necessary to determine the pattern change rules and influencing factors of grain production in upgrading HHH grain production and ensuring national food security. Here, wheat and maize occupy a larger proportion in the structure of grain production. The annual double-cropping system of winter wheat and summer maize is the most popular planting pattern. In addition, the yield of wheat and maize accounts for approximately 61% and 31% of the total national output, respectively [26]. Therefore, it is necessary to determine the pattern change rules and influencing factors of grain production in upgrading HHH grain production and ensuring national food security.

**Figure 1.** Location map of the Huang-Huai-Hai Plain (HHH). **Figure 1.** Location map of the Huang-Huai-Hai Plain (HHH).

#### *2.2. Data Sources and Processing 2.2. Data Sources and Processing*

The historical socio-economic data used in this paper were collected from the Statistical Yearbook published by the National Bureau of Statistics of China, which includes annual data on the grain yield per unit area, the sown area of grain crops, the amount of agriculture fertilizer application, the amount of effective irrigation area, the amount of pesticide application, the amount of mechanical power, and the per capita annual income of rural residents in the HHH from 1995 to 2018. According to China's statistics, grain production includes corn, rice, soybeans, wheat, potatoes, and sweet potatoes. The historical socio-economic data used in this paper were collected from the Statistical Yearbook published by the National Bureau of Statistics of China, which includes annual data on the grain yield per unit area, the sown area of grain crops, the amount of agriculture fertilizer application, the amount of effective irrigation area, the amount of pesticide application, the amount of mechanical power, and the per capita annual income of rural residents in the HHH from 1995 to 2018. According to China's statistics, grain production includes corn, rice, soybeans, wheat, potatoes, and sweet potatoes.

The historical climate data were collected from the Chinese meteorological data hub (https://data.cma.cn), including the annual precipitation, annual temperature, and annual sunshine duration of 123 meteorological stations in the HHH from 1995 to 2018. The meteorological data need The historical climate data were collected from the Chinese meteorological data hub (https: //data.cma.cn), including the annual precipitation, annual temperature, and annual sunshine duration of 123 meteorological stations in the HHH from 1995 to 2018. The meteorological data need to be processed by SQL-Server (Microsoft Corporation, Redmond, WA, USA) and ArcGIS10.2 software (Environmental

Systems Research Institute Inc, Redlands, CA, USA), which can analyze the average annual temperature, annual precipitation, and average annual sunshine hours of each city in different years.

Grain production is influenced by a variety of natural and socio-economic factors. In this paper, based on the previous literature [27–29], the grain yield per unit sown area (kg/km<sup>2</sup> ) was taken as the dependent variable, and climate and socio-economic factors were taken as the independent variables. Among them, the climate indexes included the annual average temperature (◦C), the annual average precipitation (mm), and the annual average sunshine duration (h), and the social-economic factors included the proportion of effective irrigation area (%), amount of fertilizer application per unit sown area (t/km<sup>2</sup> ), amount of pesticide application per unit sown area (t/km<sup>2</sup> ), amount of mechanical power per sown area (kw/km<sup>2</sup> ), and per capita annual income of rural residents (RMB). Table 1 summarizes the data used in this study.


#### **3. Methods**

#### *3.1. Exploratory Spatial Data Analysis (ESDA)*

Exploratory spatial data analysis is a collection of techniques for describing and visualizing spatial distributions, determining atypical locations or spatial outliers, discovering spatial associations, clusters, or hot spots, and to infer spatial characteristics or other forms of space heterogeneity [30]. In general, global and local spatial autocorrelation (or hot spots analysis) is often used to explore the spatial characteristics of observations [31].

#### 3.1.1. Global Spatial Autocorrelation

Global spatial autocorrelation is used to test the spatial correlation of the observations of spatial units within the study area [32], and is mainly measured by the Global Moran's I, which was first proposed by Moran [33]. The Moran's I can be calculated using Equation (1):

$$I = \frac{n\sum\_{i=1}^{n}\sum\_{j=1}^{n}w\_{lj}(\mathbf{x}\_{i}-\overline{\mathbf{x}})\left(\mathbf{x}\_{j}-\overline{\mathbf{x}}\right)}{\sum\_{i=1}^{n}\sum\_{j=1}^{n}w\_{lj}\sum\_{i=1}^{n}\left(\mathbf{x}\_{i}-\overline{\mathbf{x}}\right)^{2}}\tag{1}$$

where *I* represents Moran's *I*, *n* is the number of spatial units (in this study, *n* = 59), *x<sup>i</sup>* and *x<sup>j</sup>* are the observations of spatial units *I* and *j*, respectively, *x* is the average value of observations of spatial units, and wij is the spatial weight matrix, where *wij* = 1 if spatial units *I* and *j* share a common border and *wij* = 0 otherwise. The values of Global Moran's *I* range from −1 to 1. If *I* < 0, it means there is a negative spatial correlation in the space, while if *I* > 0, it means there is a positive spatial correlation, an d if *I* = 0, it means there is no spatial correlation.

The significance of Moran's *I* is usually measured by *Z* statistics using Equation (2):

$$Z(I) = \frac{I - E(I)}{\sqrt{Var(I)}}\tag{2}$$

where *E*(*I*) and *Var*(*I*) are the expected value and variance of Moran's *I*, respectively.

#### 3.1.2. Hot Spot (Getis-Ord *G<sup>i</sup>* \* ) Analysis

The Getis-Ord *G<sup>i</sup>* \* is commonly used for hot spot analysis, which can identify clustering relationships at different spatial locations. Compared with the local spatial autocorrelation, the Getis-Ord *G<sup>i</sup>* \* is more sensitive to the identification of cold and hot spots, and can fully reflect the high or low value distribution relationship between a certain geographic element and other surrounding elements [34]. The formula is [34–37]:

$$G\_{l}^{\*} = \frac{\sum\_{j=1}^{n} w\_{lj} (\mathbf{x}\_{j} - \overline{\mathbf{x}})}{\sqrt{\frac{1}{n} \sum\_{j=1}^{n} \mathbf{x}\_{j}^{2} - \overline{\mathbf{x}}^{2}} \cdot \sqrt{\frac{n}{n-1} \sum\_{j=1}^{n} w\_{lj}^{2} - \frac{n}{n-1} \left(\sum\_{j=1}^{n} w\_{lj}\right)^{2}}} \tag{3}$$

where *x* = <sup>1</sup> *n* P*n j*=1 *xj* , *n* is the number of spatial units (in this study, *n* = 59), and *wij* is the spatial weight matrix, where *wij* = 1 if spatial units *I* and *j* share a common border and *wij* = 0 otherwise.

The significance of *G<sup>i</sup>* \* is usually measured by *Z* statistics using Equation (4):

$$Z(G\_i^\*) = \frac{[G\_i^\* - E(G\_l^\*)]}{\sqrt{Var(G\_l^\*)}} \tag{4}$$

where *E*(*G<sup>i</sup>* \* ) and Var(*G<sup>i</sup>* \* ) are the expected value and variance of *G<sup>i</sup>* \* , respectively. If *Z*(*Gi*\* ) is significantly positive, it indicates that the observations around the spatial unit *i* are relatively high (higher than the average), and are high-value clusters in the space, belonging to hot spots; on the contrary, if *Z*(*G<sup>i</sup>* \* ) is significantly negative, it indicates that the observations around the spatial unit *i* are relatively low (lower than the mean), and are low-value clusters in the space, belonging to cold spots. The larger (or smaller) the *Z*(*G<sup>i</sup>* \* ) is, the more intense the clustering of high (or low) values. A *Z*(*G<sup>i</sup>* \* ) near zero indicates no apparent spatial clustering.

#### *3.2. Gravity Center Model*

The gravity center model is used to measure the overall distribution of a certain attribute in a region. It can provide a concise and accurate feature of the distribution of the attribute in the space, and can indicate the general trend and central location of its distribution. We assumed that a large region (such as an administrative region) consists of several subregions, and so the gravity center of grain production in the region could be calculated by the grain production and geographic coordinates of each sub-region. The formula is [38]:

$$X = \frac{\sum\_{i=1}^{n} M\_i X\_i}{\sum\_{i=1}^{n} M\_i}; Y = \frac{\sum\_{i=1}^{n} M\_i Y\_i}{\sum\_{i=1}^{n} M\_i} \tag{5}$$

In Equation (2), *X<sup>i</sup>* and *Y<sup>i</sup>* represent the geographic coordinates of the ith subregion. *M<sup>i</sup>* represents the grain yield per unit sown area of the subregion. *X* and *Y* represent the gravity center of grain production in a large region. Using formula (3), the moving distance of the gravity center in grain production can be obtained, which can reflect the evolution of the gravity center of a property in a region;

$$D\_{ij} = R \times \sqrt{\left(X\_i - X\_j\right)^2 + \left(Y\_i - Y\_j\right)^2} \tag{6}$$

In Equation (3), *Dij* is the gravity center movement distance (km) of grain production from *j* to *i* years. (*X<sup>i</sup>* , *Y<sup>i</sup>* ) and (*X<sup>j</sup>* , *Y<sup>j</sup>* ) are the gravity center coordinates of grain production in the *i* and *j* years. R is typically 111.111, which represents the coefficient of the spherical longitude and latitude coordinates converted to plane distance.

#### *3.3. Spatial Panel Data Model*

When the data have spatial autocorrelation effects, the residuals are no longer independent of each other; thus, it is not appropriate to use the ordinary least square regression (OLS) model. Instead, the spatial lag model (SLM) or the spatial error model (SEM) should be used for analysis. The formula is [39–41]:

$$\text{SLM}: \quad \text{y} = \rho w\_{\text{lj}} \text{y} + \text{x}\emptyset + \mu \tag{7}$$

$$\begin{array}{c} \text{SEM}: \quad \left\{ \begin{array}{l} y = \mathsf{x}\pounds + \varepsilon\\ \varepsilon = \lambda w\_{ij} + \mu \end{array} \right. \\\\ \left. \end{array} \tag{8} \\\\ \tag{9}$$

In Equations (4) and (5), *y* is the dependent variable, *x* is the explanatory variables, *Wij* is the space weight matrix, ρ is the spatial hysteresis parameter, β is the parameter vector, µ is the random interference term, ε is the regression residual vector, and λ is the autoregression parameter.

#### *3.4. Geographically Weighted Regression (GWR) Model*

Geographically weighted regression models are superior to traditional regression models such as ordinary least squares (OLS). The geographical weighted regression (GWR) model can fully consider the spatial characteristics of each influencing factor, and more accurately show the spatial relationship between independent and dependent variables [42]. The form of a GWR model is as follows:

$$\mathcal{Y}\_{i} = \beta\_{0}(u\_{i\prime}v\_{i}) + \sum\_{\lambda} {}^{n} \beta\_{\lambda}(u\_{i\prime}v\_{i}) {}^{\lambda}\_{i\lambda} + \varepsilon \tag{9}$$

In Equation (6), Y*<sup>i</sup>* represents the grain production in region *i*, β0*(ui,vi)* represents a constant, βλ*(ui,vi)* represents the regression coefficient, *(u<sup>i</sup>* ,*vi )* represents the geographic location of the cities *i*, *Xi*<sup>λ</sup> represents the parameter value of the λ independent variable of city *i*, and ε represents the random error. The optimal bandwidth distance can be obtained automatically in GWR4.0 corrected by finite correction of the Akaike Information Criterion (AICc). The smaller the AICc value, the higher the goodness of fit of the model will be [43].

#### **4. Results**

#### *4.1. Temporal Changes of Grain Production in the HHH*

From 1995 to 2018, the grain yield per unit of sown area in the HHH has shown a steady increase (Figure 2), which can be divided into three stages, as follows: the fluctuating growth stage (1995–2005), the steady growth stage (2005–2015), and the slow descent stage (2015–2018). Firstly, in the fluctuating growth stage (1995–2005), although the grain yield per unit sown area has increased in this stage, the fluctuation range is relatively large and the grain yield per unit of sown area gradually stabilized in 2005. Secondly, during the steady growth stage (2005–2015), the grain yield per unit of sown area showed a characteristic of small fluctuation growth. Finally, a slow decline began to appear in the grain yield per unit of sown area in the HHH during the slow descent stage (2015–2018).

**Figure 2.** The temporal changes of grain production in the HHH from 1995 to 2018. **Figure 2.** The temporal changes of grain production in the HHH from 1995 to 2018. **Figure 2.** The temporal changes of grain production in the HHH from 1995 to 2018.

#### *4.2. Spatial Characteristics of Grain Production in the HHH 4.2. Spatial Characteristics of Grain Production in the HHH*

*4.2. Spatial Characteristics of Grain Production in the HHH*  We classified the grain yield data of each urban unit of sown area into four types according to 2000–3500, 3500–5000, 5000–6500, and 6500–8000 kg/km2 using ArcGIS10.2 software. For the convenience of analysis, four time sections of 1995, 2005, 2015, and 2018 were selected for study, for We classified the grain yield data of each urban unit of sown area into four types according to 2000–3500, 3500–5000, 5000–6500, and 6500–8000 kg/km<sup>2</sup> using ArcGIS10.2 software. For the convenience of analysis, four time sections of 1995, 2005, 2015, and 2018 were selected for study, for which the spatial distribution characteristics of the grain production pattern were discussed (Figure 3). We classified the grain yield data of each urban unit of sown area into four types according to 2000–3500, 3500–5000, 5000–6500, and 6500–8000 kg/km<sup>2</sup> using ArcGIS10.2 software. For the convenience of analysis, four time sections of 1995, 2005, 2015, and 2018 were selected for study, for which the spatial distribution characteristics of the grain production pattern were discussed (Figure 3).

which the spatial distribution characteristics of the grain production pattern were discussed (Figure 3). It can be seen from Figure 3 that the spatial variation in grain yield per unit of sown area in each city is very significant, and the overall grain yield shows an increasing trend. Specifically, the number of cities in the HHH where the grain yield per unit of sown area remained within the range of 2000– 3500 kg/km2 continued to decrease, and was mainly distributed in Zhangjiakou, Cangzhou, Zhengzhou, Luoyang, Nanyang, Sanmenxia, and Bozhou in 1995, Zhangjiakou and Sanmenxia in It can be seen from Figure 3 that the spatial variation in grain yield per unit of sown area in each city is very significant, and the overall grain yield shows an increasing trend. Specifically, the number of cities in the HHH where the grain yield per unit of sown area remained within the range of 2000–3500 kg/km<sup>2</sup> continued to decrease, and was mainly distributed in Zhangjiakou, Cangzhou, Zhengzhou, Luoyang, Nanyang, Sanmenxia, and Bozhou in 1995, Zhangjiakou and Sanmenxia in 2005, It can be seen from Figure 3 that the spatial variation in grain yield per unit of sown area in each city is very significant, and the overall grain yield shows an increasing trend. Specifically, the number of cities in the HHH where the grain yield per unit of sown area remained within the range of 2000–3500 kg/km<sup>2</sup> continued to decrease, and was mainly distributed in Zhangjiakou, Cangzhou, Zhengzhou, Luoyang, Nanyang, Sanmenxia, and Bozhou in 1995, Zhangjiakou and Sanmenxia in 2005, and only in Zhangjiakou in 2015. In 2018, the number of cities of this type decreased to zero.

2005, and only in Zhangjiakou in 2015. In 2018, the number of cities of this type decreased to zero. The number of cities within the range of grain yield per unit of sown area of 3500 to 5000 kg/km2 also showed a downward trend, and the type area shrank from 30 cities in 1995 to 5 cities in 2018, and only in Zhangjiakou in 2015. In 2018, the number of cities of this type decreased to zero. The number of cities within the range of grain yield per unit of sown area of 3500 to 5000 kg/km<sup>2</sup> also showed a downward trend, and the type area shrank from 30 cities in 1995 to 5 cities in 2018, and presented a layout trend of agglomeration to dispersion in spatial distribution. The number of cities within the range of grain yield per unit of sown area of 3500 to 5000 kg/km<sup>2</sup> also showed a downward trend, and the type area shrank from 30 cities in 1995 to 5 cities in 2018, and presented a layout trend of agglomeration to dispersion in spatial distribution.

and presented a layout trend of agglomeration to dispersion in spatial distribution. The number of cities in the range of 5000–6500 kg/km2 showed a rising–falling–rising trend. Among them, the number of cities in this type of area increased from 19 to 35 in 1995–2005, reduced from 35 to 28 in 2005–2015, and increased from 28 to 40 in 2015–2018. The spatial distribution of this type overall showed a tendency varying between scattered and clustered development, and appeared The number of cities in the range of 5000–6500 kg/km<sup>2</sup> showed a rising–falling–rising trend. Among them, the number of cities in this type of area increased from 19 to 35 in 1995–2005, reduced from 35 to 28 in 2005–2015, and increased from 28 to 40 in 2015–2018. The spatial distribution of this type overall showed a tendency varying between scattered and clustered development, and appeared The number of cities in the range of 5000–6500 kg/km<sup>2</sup> showed a rising–falling–rising trend. Among them, the number of cities in this type of area increased from 19 to 35 in 1995–2005, reduced from 35 to 28 in 2005–2015, and increased from 28 to 40 in 2015–2018. The spatial distribution of this type overall showed a tendency varying between scattered and clustered development, and appeared roughly from east to west, and from south to north.

roughly from east to west, and from south to north. The trend of unit of area sown to grain production maintained in the 6500–8000 kg/km2 interval of urban change was different from the above three kinds. In 1995, the number of cities within this range was 3, which increased to 5 in 2005 and 24 in 2015, and reduced to 14 in 2018, which reflects a characteristic of a sharp increase followed by a slow decline. However, they were still mainly roughly from east to west, and from south to north. The trend of unit of area sown to grain production maintained in the 6500–8000 kg/km<sup>2</sup> interval of urban change was different from the above three kinds. In 1995, the number of cities within this range was 3, which increased to 5 in 2005 and 24 in 2015, and reduced to 14 in 2018, which reflects a characteristic of a sharp increase followed by a slow decline. However, they were still mainly The trend of unit of area sown to grain production maintained in the 6500–8000 kg/km<sup>2</sup> interval of urban change was different from the above three kinds. In 1995, the number of cities within this range was 3, which increased to 5 in 2005 and 24 in 2015, and reduced to 14 in 2018, which reflects a characteristic of a sharp increase followed by a slow decline. However, they were still mainly distributed in the central and southern parts of the HHH.

distributed in the central and southern parts of the HHH. In terms of spatial distribution, the areas with high grain yield per unit of sown area in the Huang-Huia-Hai Plain were mainly distributed in the east and south, while the areas with low grain yield per unit of sown area were mainly distributed in the west and north of the HHH, which indicated that the grain production capacity of the eastern and southern regions of the HHH was distributed in the central and southern parts of the HHH. In terms of spatial distribution, the areas with high grain yield per unit of sown area in the Huang-Huia-Hai Plain were mainly distributed in the east and south, while the areas with low grain yield per unit of sown area were mainly distributed in the west and north of the HHH, which indicated that the grain production capacity of the eastern and southern regions of the HHH was higher, while that of the northern and western regions was lower. In terms of spatial distribution, the areas with high grain yield per unit of sown area in the Huang-Huia-Hai Plain were mainly distributed in the east and south, while the areas with low grain yield per unit of sown area were mainly distributed in the west and north of the HHH, which indicated that the grain production capacity of the eastern and southern regions of the HHH was higher, while that of the northern and western regions was lower.

higher, while that of the northern and western regions was lower.

**Figure 3.** The spatial distribution of grain production in the HHH in (**a**) 1995, (**b**) 2005, (**c**) 2015, and (**d**) 2018. **Figure 3.** The spatial distribution of grain production in the HHH in (**a**) 1995, (**b**) 2005, (**c**) 2015, and (**d**) 2018.

#### *4.3. Dynamic Change of Barycenter of Grain Production in the HHH*

*4.3. Dynamic Change of Barycenter of Grain Production in the HHH*  Figure 3 reflects the static distribution pattern of grain production in the HHH in 1995, 2005, 2015, and 2018, but does not reflect the dynamic change trend. Therefore, we used Equation (2) to calculate the barycenter of grain production in the HHH from 1995 to 2018 (Figure 4), and then used Equation (3) to calculate the barycenter movement distance (Table 2) to analyze the dynamic change Figure 3 reflects the static distribution pattern of grain production in the HHH in 1995, 2005, 2015, and 2018, but does not reflect the dynamic change trend. Therefore, we used Equation (2) to calculate the barycenter of grain production in the HHH from 1995 to 2018 (Figure 4), and then used Equation (3) to calculate the barycenter movement distance (Table 2) to analyze the dynamic change characteristics of the grain production pattern.

characteristics of the grain production pattern. According to Figure 4, the grain production center in the HHH generally shifted from southeast to northwest in Tai'an, and gradually stabilized in the central area of the HHH with the passage of time. Specifically, from 1995 to 1997, the grain production center of HHH was distributed in Jining city, which moved first to the northwest and then to the southwest. From 1998 to 2000, the movement direction of the barycenter in grain production remained highly stable; that is, it continued to move in the northwest direction. The barycenter of grain production began to enter Tai'an City. From 2000 to 2001, the movement direction of the barycenter in grain production was reversed and shifted to According to Figure 4, the grain production center in the HHH generally shifted from southeast to northwest in Tai'an, and gradually stabilized in the central area of the HHH with the passage of time. Specifically, from 1995 to 1997, the grain production center of HHH was distributed in Jining city, which moved first to the northwest and then to the southwest. From 1998 to 2000, the movement direction of the barycenter in grain production remained highly stable; that is, it continued to move in the northwest direction. The barycenter of grain production began to enter Tai'an City. From 2000 to 2001, the movement direction of the barycenter in grain production was reversed and shifted to the southeast.

the southeast. In 2001–2003, the barycenter of grain production shifted first to the southwest and then to the northeast, and in 2003–2005, it shifted first to the southwest and then to the southeast. From 2005 to 2011, the barycenter of grain production fluctuated in all directions. From 2011 to 2013, the grain production center assumed a similar change trend to that from 2001 to 2003. From 2013 to 2015, the In 2001–2003, the barycenter of grain production shifted first to the southwest and then to the northeast, and in 2003–2005, it shifted first to the southwest and then to the southeast. From 2005 to 2011, the barycenter of grain production fluctuated in all directions. From 2011 to 2013, the grain production center assumed a similar change trend to that from 2001 to 2003. From 2013 to 2015, the grain production center first moved to the northeast and then to the northwest.

grain production center first moved to the northeast and then to the northwest. From 2015 to 2018, the grain production center of gravity showed characteristics of moving from From 2015 to 2018, the grain production center of gravity showed characteristics of moving from southeast to northwest, but it was still stable in the territory of Tai'an City, that is, the southeast of the

HHH. The dynamic shift in the grain production center in the HHH indicates that the regional grain production capacity has the characteristics of non-stationarity in time and non-equilibrium in space at the same time, and the shift in the grain production center from southeast to northwest indicates that the grain production capacities in the western and northern parts of the HHH were continuously enhanced. In addition, the grain production center in the HHH tended to be stable over time, and it was concentrated in the border area between Jining and Tai'an. the HHH. The dynamic shift in the grain production center in the HHH indicates that the regional grain production capacity has the characteristics of non-stationarity in time and non-equilibrium in space at the same time, and the shift in the grain production center from southeast to northwest indicates that the grain production capacities in the western and northern parts of the HHH were continuously enhanced. In addition, the grain production center in the HHH tended to be stable over time, and it was concentrated in the border area between Jining and Tai'an.

*Int. J. Environ. Res. Public Health* **2020**, *17*, x 9 of 19

**Figure 4.** The track of the gravity center of grain production in the HHH from 1995 to 2018. **Figure 4.** The track of the gravity center of grain production in the HHH from 1995 to 2018.

Table 2 shows that the barycenter of grain production in the HHH as a whole moved from the southeast to northwest from 1995 to 2005, with a distance of 17.7 km. From 2005 to 2018, the barycenter of grain production moved to the northwest with a distance of 9.4 km, which was significantly smaller than that from 1995 to 2005, which confirmed that the barycenter of grain production in the HHH showed good time-stability characteristics over time; however, this could not cover up the disequilibrium in the spatial characteristics of grain production. On the whole, from 1995 to 2018, the center of gravity of grain production moved 26.1 km to the northwest. The center was still stable in the territory of Tai'an city, that is, to the southeast of the HHH. Table 2 shows that the barycenter of grain production in the HHH as a whole moved from the southeast to northwest from 1995 to 2005, with a distance of 17.7 km. From 2005 to 2018, the barycenter of grain production moved to the northwest with a distance of 9.4 km, which was significantly smaller than that from 1995 to 2005, which confirmed that the barycenter of grain production in the HHH showed good time-stability characteristics over time; however, this could not cover up the disequilibrium in the spatial characteristics of grain production. On the whole, from 1995 to 2018, the center of gravity of grain production moved 26.1 km to the northwest. The center was still stable in the territory of Tai'an city, that is, to the southeast of the HHH.


**Table 2.** The gravity center changes in grain production in 1995–2018.

"-" means the item does not exist.

#### *4.4. Spatial Correlation Characteristics of Grain Production Pattern in the HHH*

#### 4.4.1. Global Spatial Correlation Characteristics

Based on the grain yield data per unit sown area, the Moran's *I* value, *Z* statistic, and *P*-value were calculated using Geoda software, and the spatial correlation characteristics of grain production are shown in Table 3.

It can be seen from Table 3 that the Moran's *I* values from 1995 to 2018 were all greater than 0 and significant at the threshold level of 5%, indicating that the grain yield per unit sown area in the HHH was not randomly distributed but positively correlated. This indicates that the grain production layout showed strong spatial clustering characteristics. From 1995 to 1997, the Moran's *I* value decreased from 0.4114 to 0.1718. This indicates that the agglomeration and development trend in grain production in the HHH weakened during this period. From 1997 to 2002, the Moran's *I* value showed a fluctuating trend, ranging from 0.1718 to 0.3356. In this period, the grain production experienced a process of both agglomeration and dispersion development, but the change trend was small.


**Table 3.** Moran's *I* value of grain production in HHH.

\*\*\* means significant at the 1% threshold level; \*\* means significant at the 5% threshold level.

From 2002 to 2006, the Moran's *I* value showed "up–down" fluctuation characteristics twice, and the change trend was more intense. This indicates that the grain production pattern of the HHH was greatly changed during this period. From 2007 to 2011, the increase in Moran's *I* value indicated that the grain production distribution in the HHH was in an increasingly concentrated state in this period. From 2011 to 2012, the Moran's *I* value changed from 0.2897 to 0.2401, a decrease by 0.0496 units, indicating that the clustering characteristics of grain production layout weakened during this period.

The Moran's I value showed a continuous rising trend, indicating that the clustering characteristics of grain production distribution were enhanced from 2012 to 2014. The value of Moran's *I* changed from 0.3727 to 0.2315 in 2014–2018, indicating that the agglomeration characteristics of grain production distribution during this period weakened. On the whole, the Moran's *I* value experienced a change process of "down–up–down" in this period; however, the agglomeration and distribution characteristics of grain production did not change.

#### 4.4.2. Local Spatial Correlation Characteristics

and (**d**) 2018.

The evaluation of the global spatial correlation feature has the defect of ignoring the instability of local spatial processes. Therefore, the local spatial correlation characteristics of grain production in the HHH can be analyzed by observing the *G<sup>i</sup>* \* index in 1995, 2005, 2015, and 2018 (Figure 5). *Int. J. Environ. Res. Public Health* **2020**, *17*, x 12 of 19

**Figure 5.** Spatial correlation characteristics of grain production in the HHH. (**a**) 1995, (**b**) 2005, (**c**) 2015, **Figure 5.** Spatial correlation characteristics of grain production in the HHH. (**a**) 1995, (**b**) 2005, (**c**) 2015, and (**d**) 2018.

*4.5. Influencing Factors for the Changes of Grain Production in the HHH*  4.5.1. Analysis of the Spatial Spillover Effect of the Influencing Factors First, SPSS software was used to eliminate the collinearity of variables, and finally four indexes were extracted, namely, the effective irrigation area (EIA), amount of fertilizer application per unit of sown area (AFA), per capita annual income of rural residents (PCI), and annual average precipitation (PRE). In addition, 3.3 proves that the grain production pattern in the HHH has dependent According to Figure 5, the numbers of hot spots, sub-hot spots, sub-cold spots and cold spots were 13, 21, 21, and 4, respectively, in 1995, and the cold spots were mainly distributed in Sanmenxia, Luoyang, and Nanyang. The sub-cold spots were distributed in Beijing, Zhangjiakou, Chengde, Baoding, Jiaozuo, Xinxiang, Zhengzhou, Kaifeng, Zhoukou, Zhumadian, Xinyang, Fuyang, Bengbu, Huainan, and Bozhou. Sub-hot spots were concentrated in Tianjin, Qinhuangdao, Langfang, Cangzhou, Hengshui, Xingtai, Shijiazhuang, Dezhou, Binzhou, Dongying, Heze, Anyang, Hebi, and Puyang. Hot spots were distributed in the eastern part of the HHH. In particular, Yantai, Qingdao, Weifang, Jinan, Rizhao, Linyi, Zaozhuang, Jining, Lianyungang and Xuzhou were the most typical ones.

characteristics; therefore, the influence of the spatial spillover effect cannot be ignored. Secondly, the GeoDa software was used to obtain the parameter estimation results of the spatial lag model, spatial In 2005, the numbers of cities in the four categories were 16, 14, 16, and 13. Compared with 1995, the number of cities in hot spots and cold spots increased, indicating that the agglomeration of grain

irrigation and water conservancy facilities in the HHH is the reason for this phenomenon.

According to Table 4, the coefficients of EIA were 0.505, 0.415, 0.532, and 0.588 in 1995, 2005, 2015, and 2018, respectively, and were effective at the 1% threshold level. This indicates that there was a positive correlation between the EIA and grain yield; that is, an increase in the EIA will bring about an increase in the grain yield. At the same time, the coefficient of EIA on the whole was on the rise, indicating that it has a stronger positive effect on grain yield. The continuous improvement of

The coefficients of fertilizer application per unit of sown area in 1995 and 2005 were 0.008 and 0.208, respectively, which were both significant at the 1% threshold level, and the coefficients in 2015 and 2018 were 0.124 and 0.023, respectively, and were significant at the 5% threshold level. This shows that the positive effect of chemical fertilizer application per unit of sown area on grain production increased from 1995 to 2005. The cold spots were concentrated in the north and southwest of the HHH. The sub-cold spots were distributed in Qinhuangdao, Tangshan, Tianjin, Langfang, Baoding, Xinyang, Fuyang, Zhoukou, Xuchang, Zhengzhou, and Jiaozuo. The sub-hot spots continued to shrink, and were mainly distributed in Shijiazhuang, Cangzhou, and Hengshui. The spatial distribution of hot spots is clear, and the main distribution was in the middle and eastern part of the HHH. In 2015, the numbers of hot spots, sub-hot spots, sub-cold spots, and cold spots were 13, 20, 14, and 12 cities, respectively.

Compared with 2005, the number of cities in hot spots and cold spots decreased by 3 and 1, respectively, indicating that the local spatial clustering characteristics of grain production in the HHH were weakened from 2005 to 2015. In 2018, the hot spots, sub-hot spots, sub-cold spots, and cold spots included 10, 17, 22, and 10 cities, respectively. Compared with 2015, the number of cities in hot spots and cold spots decreased by three and two, and the number of cities in sub-hot spots increased by eight. This indicates that the local spatial agglomeration characteristics of grain production in the Huang-Huai-Hai Plain gradually weakened from 2005. In general, hot spots spread from the east of the HHH to the southeast and the central region, and sub-hot spots were mainly distributed in the central region. The cold spots and sub-cold spots were mainly distributed in the north and south regions.

#### *4.5. Influencing Factors for the Changes of Grain Production in the HHH*

#### 4.5.1. Analysis of the Spatial Spillover Effect of the Influencing Factors

First, SPSS software was used to eliminate the collinearity of variables, and finally four indexes were extracted, namely, the effective irrigation area (EIA), amount of fertilizer application per unit of sown area (AFA), per capita annual income of rural residents (PCI), and annual average precipitation (PRE). In addition, 3.3 proves that the grain production pattern in the HHH has dependent characteristics; therefore, the influence of the spatial spillover effect cannot be ignored. Secondly, the GeoDa software was used to obtain the parameter estimation results of the spatial lag model, spatial error model, and OLS model. Finally, the statistical values of LMLAG and R-LMLAG in the OLS results were significantly higher than those of LMERR and R-LMERR at the 10% level; as such, it was appropriate to use the spatial lag model to explore the key factors affecting grain production (Table 4).

According to Table 4, the coefficients of EIA were 0.505, 0.415, 0.532, and 0.588 in 1995, 2005, 2015, and 2018, respectively, and were effective at the 1% threshold level. This indicates that there was a positive correlation between the EIA and grain yield; that is, an increase in the EIA will bring about an increase in the grain yield. At the same time, the coefficient of EIA on the whole was on the rise, indicating that it has a stronger positive effect on grain yield. The continuous improvement of irrigation and water conservancy facilities in the HHH is the reason for this phenomenon.

The coefficients of fertilizer application per unit of sown area in 1995 and 2005 were 0.008 and 0.208, respectively, which were both significant at the 1% threshold level, and the coefficients in 2015 and 2018 were 0.124 and 0.023, respectively, and were significant at the 5% threshold level. This shows that the positive effect of chemical fertilizer application per unit of sown area on grain production experienced a changing process of decline after rising first, reflecting that the increase in chemical fertilizer application per sown area from 1995 to 2005 significantly increased the grain production in the HHH.

From 2005 to 2018, the boosting effect on grain production was alleviated. The reason for this phenomenon may be that the increase in the use of chemical fertilizers in the low-level agricultural development stage has a significant effect on increasing grain production. As the amount of fertilizer input remains at a high level, the effect of increasing the amount of fertilizer input in the future is limited as regards improving the grain yield. The coefficients of rural residents' per capita annual income in the four time sections were 0.405, 0.170, 0.124, and 0.050, which passed the significance tests at the critical value levels of 1%, 5%, 10%, and 10%, respectively.

The increase in per capita annual income continued to weaken the boost to grain production. The reason is that agriculture has been the main source of income for Chinese farmers to maintain

their livelihoods for a long time. As farmers' income levels increase, they have more funds to purchase agricultural machinery, fertilizers, pesticides, seeds, and other production materials, thereby achieving the goal of increasing grain yield. With the rapid advancement of urbanization, farmers' income channels have become increasingly diversified, which has greatly reduced their dependence on agriculture, and the tendency of farmers to invest in non-agriculture has become more obvious. Therefore, to some extent, the increase in PCI has an inhibitory effect on grain output.

The coefficient of PRE was 0.508 in 1995 and passed the significance test of the 1% critical value level, indicating that the increase in PRE in that year played a promoting role in grain production. The coefficients of PRE were 0.016 and 0.148 in 2005 and 2015, both of which failed to pass the 10% significance test, indicating that the PRE increase did not have an obvious positive promotion effect on grain production. The underlying reason may be that with the continuous improvement of irrigation facilities, the impact of changes in PRE on grain production became weaker.


**Table 4.** Regression analysis results of the spatial lag model.

\*\*\* means significant at the 1% threshold level; \*\* means significant at the 5% threshold level; \* means significant at the 10% threshold level.

#### 4.5.2. Analysis of Spatial Heterogeneity of Influencing Factors

Based on the above research findings, the EIA, AFA, and PCI had a significant impact on the grain production in the HHH. However, SLM cannot explain the specific degree and scope of the impact of these three factors in space, and so a GWR model was explored to investigate the spatial difference in the influence of these three factors. The three factors—EIA, AFA, and PCI—were put into the model according to the criteria of AICc minimization.

The spatial distribution of Local R-squared values derived from the GWR model is displayed in Figure 6. The GWR results show that three factors explain 81.1% of the variance in grain production. Geographic variations in these factors describe a difference in the combined statistical influence of the three variables on grain production across cities in the HHH, from a very weak relationship (near 0.30) to a strong relationship (>0.80). We found that 50.94% of cities maintained local R-squared values of more than 50%. The predictive power of the model shows characteristics of increasing from east to west. The local R-squared map suggests that the predictive power of the analysis was greatest in relation to the northern part (Chengde, Baoding, and Zhangjiakou) and the southwestern part (Jiaozuo, Jiyuan, and Sanmenxia). The lower R-squared values demonstrate a poorer regression fit in the eastern parts of the HHH, such as Weihai, Yantai, and Qingdao.

To explore the strength of the influence of each of the three factors on grain production, we created maps for each factor, which represent the geographic distribution of their regression coefficient values across HHH, according to the results of the GWR modeling. The mapped regression coefficients are divided into five classifications through Natural Breaks. In Figure 7, white is used to indicate cities without data, and gradient shading is used to show cities with a significant relationship between variables and grain production.

**Figure 6.** R-squared values derived from the geographical weighted regression (GWR) model. **Figure 6.** R-squared values derived from the geographical weighted regression (GWR) model.

To explore the strength of the influence of each of the three factors on grain production, we created maps for each factor, which represent the geographic distribution of their regression coefficient values across HHH, according to the results of the GWR modeling. The mapped regression coefficients are divided into five classifications through Natural Breaks. In Figure 7, white is used to indicate cities without data, and gradient shading is used to show cities with a significant relationship between variables and grain production. The proportion of the effective irrigated area had a positive impact on grain production in the HHH, and its impact intensity presents a spatial distribution characteristic of "low west and high east". A total of 75.47% of the cities showed a significant positive relationship between the effective irrigate area and production, mainly in Hebei and Henan provinces. The reason may be that these two provinces are mostly inland, with relatively drier climate and less rainfall; therefore, irrigation is mainly used to supply the water requirement of crops.

The proportion of the effective irrigated area had a positive impact on grain production in the HHH, and its impact intensity presents a spatial distribution characteristic of "low west and high east". A total of 75.47% of the cities showed a significant positive relationship between the effective irrigate area and production, mainly in Hebei and Henan provinces. The reason may be that these two provinces are mostly inland, with relatively drier climate and less rainfall; therefore, irrigation is mainly used to supply the water requirement of crops. AFA had a positive effect on 94.34% of the cities in the HHH, which can significantly increase the food production of 39.62% of the cities, mainly in the central, northern, and southern regions of the HHH, such as Zhangjiakou, Puyang, and Fuyang. AFA had a negative effect on 5.66% of the cities, AFA had a positive effect on 94.34% of the cities in the HHH, which can significantly increase the food production of 39.62% of the cities, mainly in the central, northern, and southern regions of the HHH, such as Zhangjiakou, Puyang, and Fuyang. AFA had a negative effect on 5.66% of the cities, mainly in the eastern part of the HHH, such as Qingdao and Weifang. The impact of AFA on grain production is limited. In particular, with the long-term investment of chemical fertilizers by Chinese farmers on cultivated land, the impact of various chemical fertilizers applied by farmers on soil fertility has become nearly saturated. The majority of the chemical fertilizers play a role in maintaining soil fertility after application, so even if the input of chemical fertilizers is increased, the positive effect on crop yield is not clear enough.

mainly in the eastern part of the HHH, such as Qingdao and Weifang. The impact of AFA on grain production is limited. In particular, with the long-term investment of chemical fertilizers by Chinese farmers on cultivated land, the impact of various chemical fertilizers applied by farmers on soil fertility has become nearly saturated. The majority of the chemical fertilizers play a role in maintaining soil fertility after application, so even if the input of chemical fertilizers is increased, the positive effect on crop yield is not clear enough. The per capita income of rural residents had a positive impact on 86.79% of the cities in the HHH, of which only 9.43% passed the significance test, mainly in the northeastern part of the HHH, such as Chengde, Qinhuangdao, Tangshan, Langfang, and Cangzhou. The per capita income of rural The per capita income of rural residents had a positive impact on 86.79% of the cities in the HHH, of which only 9.43% passed the significance test, mainly in the northeastern part of the HHH, such as Chengde, Qinhuangdao, Tangshan, Langfang, and Cangzhou. The per capita income of rural residents had a negative impact on 13.21% of cities, and was concentrated in the southwestern region of HHH, that is, the western and southern regions of Henan Province, such as Sanmenxia, Nanyang, and Xinyang. As rural residents flow into developed urban areas, the non-agricultural income they earn from moving into cities has gradually become an important source of livelihood, and their dependence and emphasis on agriculture has gradually decreased; abandonment can even occur, which causes a huge negative impact on grain production.

residents had a negative impact on 13.21% of cities, and was concentrated in the southwestern region of HHH, that is, the western and southern regions of Henan Province, such as Sanmenxia, Nanyang, and Xinyang. As rural residents flow into developed urban areas, the non-agricultural income they earn from moving into cities has gradually become an important source of livelihood, and their dependence and emphasis on agriculture has gradually decreased; abandonment can even occur,

**Figure 7.** The spatial distribution of regression coefficients for the three factors based on a GWR model for 2018 in the HHH. **Figure 7.** The spatial distribution of regression coefficients for the three factors based on a GWR model for 2018 in the HHH.

#### **5. Discussions**

**5. Discussions**  Even our results suggested that the effect of precipitation on grain production became weaker, we should also be aware of the relationship between rainfall and groundwater; that is, rainfall becomes groundwater through infiltration, providing sufficient water for irrigating crops. With the improvement in irrigation facilities, irrigation has gradually become an indispensable and important means for stable grain production. This may also be the reason why the direct impact of PRE on grain production is gradually weakening and why the EIA is gradually increasing. Future research should focus on related research in this area. The reasons for the insignificant effect of temperature and sunshine duration may be attributed to two points: first, compared with precipitation, the heat resources of the HHH can well meet the growth needs of winter wheat and summer corn, and the yield is less affected by temperature and sunshine duration. At the same time, an increase in temperature and sunshine duration may increase evaporation, thereby offsetting the effect of precipitation. It is also possible to attribute this effect to precipitation rather than temperature and sunshine duration. Second, it may be that the crops in the research area are not subdivided. Different crops have inconsistent requirements for temperature and sunshine duration, which may weaken the effects of temperature and sunshine duration. This will also be the focus of future research. China's food self-sufficiency rate has reached more than 95%. The concept of "with grain in the hand, the heart is unharried" has been made reality. However, in the face of the impact of global warming and COVID-19, as well as the requirements of new urbanization, as a country with one of the largest populations in the world, ensuring the stability of China's food supply is not only related to Even our results suggested that the effect of precipitation on grain production became weaker, we should also be aware of the relationship between rainfall and groundwater; that is, rainfall becomes groundwater through infiltration, providing sufficient water for irrigating crops. With the improvement in irrigation facilities, irrigation has gradually become an indispensable and important means for stable grain production. This may also be the reason why the direct impact of PRE on grain production is gradually weakening and why the EIA is gradually increasing. Future research should focus on related research in this area. The reasons for the insignificant effect of temperature and sunshine duration may be attributed to two points: first, compared with precipitation, the heat resources of the HHH can well meet the growth needs of winter wheat and summer corn, and the yield is less affected by temperature and sunshine duration. At the same time, an increase in temperature and sunshine duration may increase evaporation, thereby offsetting the effect of precipitation. It is also possible to attribute this effect to precipitation rather than temperature and sunshine duration. Second, it may be that the crops in the research area are not subdivided. Different crops have inconsistent requirements for temperature and sunshine duration, which may weaken the effects of temperature and sunshine duration. This will also be the focus of future research. China's food self-sufficiency rate has reached more than 95%. The concept of "with grain in the hand, the heart is unharried" has been made reality. However, in the face of the impact of global warming and COVID-19, as well as the requirements of new urbanization, as a country with one of the largest populations in the world, ensuring the stability of China's food supply is not only related to national security, but is also related to the stability of the world.

national security, but is also related to the stability of the world. Therefore, based on the results of this paper, the following policy suggestions are put forward to increase grain production in the HHH, an important grain production base for food security in China to maintain China's food security. The suggestions include the following: to increase the construction investment for basic farmland infrastructures, such as irrigation facilities; to cultivate and promote good varieties and treatments, and implement soil testing and formula fertilization; to standardize the rural land market; to promote the transfer of rural land in an orderly manner; and to Therefore, based on the results of this paper, the following policy suggestions are put forward to increase grain production in the HHH, an important grain production base for food security in China to maintain China's food security. The suggestions include the following: to increase the construction investment for basic farmland infrastructures, such as irrigation facilities; to cultivate and promote good varieties and treatments, and implement soil testing and formula fertilization; to standardize the rural land market; to promote the transfer of rural land in an orderly manner; and to realize the large-scale management of cultivated land.

realize the large-scale management of cultivated land. However, the impact of relevant agricultural policies issued by the country also needs to be considered in the future. From 1995 to 2018, the overall increase in grain production in HHH and the gradual reduction in spatial differences well reflects the background and national policies of the country in different periods. Since 1995, with the advancement of agricultural technology, the agricultural development of the HHH has been weakened by natural conditions, and the grain production has increased significantly [13]. At the same time, due to the popularization of fertilizers, pesticides, and irrigation, the gap in grain production in various regions in HHH is also narrowing.

Due to rapid urbanization, the conversion of fertile irrigated land to non-agricultural land seems to pose a potential threat to the food security of the HHH, and even to the whole of China [10]. Moreover, farmers are gradually migrating to cities in search of higher incomes due to the urban–rural development gap [44], which has caused the arable land in the HHH to be abandoned, resulting in the expansion of spatial differences in grain production in different regions in the HHH, and also threatening national food security. In order to ensure national food security, the central government proposed the construction of "high-standard basic farmland projects" and "agricultural modernization" to promote the large-scale, intensive, and modernized management of arable land in the HHH, thereby increasing the grain production of the HHH, and promoting the development of sustainable agriculture.

At the same time, we should not ignore the resistance, created by resource endowments, to sustainable agriculture in the HHH. In particular, the precipitation cannot meet the water demand of the crops [8,45], and water shortage is one of the major factors threatening the high and stable production of wheat [20]. The water consumption greatly exceeds the precipitation, and groundwater must be extracted to make up for the deficiency so as to maintain high yields [21]. Some areas in the HHH even appear salinized due to unreasonable irrigation [10], which poses huge challenges for the minimization of environmental impacts and the development of sustainable agriculture [46,47]. Therefore, the government must not only build irrigation facilities, but more importantly, must promote water-saving irrigation technologies, improve water resource utilization efficiency, implement drip irrigation and sprinkler irrigation [46], etc., so as to achieve sustainable agricultural production in the HHH. At the same time, in the long run, in the face of the encroachment of arable land in the promotion of urbanization, the amount of arable land in the future will also face severe challenges [10]. Promoting intensive agricultural production and improving the level of intensive use of agricultural production are also of great significance to ensure future food production [48]. Moreover, these actions can also enhance the ability to respond to natural disasters in the future. However, a sustainable food and agriculture system is one which is environmentally sound, economically viable, socially responsible, nonexploitative, and which serves as the foundation for future generations [49–51]. With the long-term development of intensive agriculture production in the Huang-Huai-Hai Plain, agricultural practices ranging from the development of irrigation projects to the use of agrichemicals have often had negative environmental impacts, such as wildlife kills, pesticide residues in drinking water, soil erosion, groundwater depletion, and salinization [52]. Substituting environmentally sound inputs for those which are damaging is an important step in addressing these problems [49]. In view of this, the Ministry of agriculture of China started the construction of the Key Laboratory of Agricultural Environment in the Huang-Huai-Hai Plain, with the objectives of scientific research, environmental monitoring, detection analysis and technical services, in 2012, aiming to carry out research on regional agricultural pollution prevention by means of agricultural non-point source pollution prevention, the environmental protection of producing areas, and the development and application of environment-friendly inputs. However, for the farmers who are the main body of agricultural production, whether these agricultural technology inputs will increase agricultural production costs and reduce agricultural income will be an important factor affecting the promotion of agricultural technology and the development of sustainable agriculture. Therefore, whether the economic, social and environmental benefits generated by agricultural production in the Huang-Huai-Hai Plain under the influence of agriculture technology can achieve a balance will be the focus of our future research.

#### **6. Conclusions**

In this paper, exploratory spatial data analysis, the gravity center model, and the spatial lag model were used to explore the spatial–temporal variation and influencing factors of the grain production pattern in the HHH from 1995 to 2018. The main conclusions were drawn as follows:

The grain production pattern in the HHH has the characteristics of being non-equilibrated in space and non-stationary in time. The spatial non-equilibrium is reflected in the shift of the grain production center from the southeast to the northwest of Tai'an city. The high-level areas of grain production capacity were mainly distributed in the east and south, while the low-level areas were distributed in the west and north. The non-stationarity of time is reflected in the rising trend in the grain production capacity and the weakening of the non-stationarity of time in the grain production center over time;

The global and local spatial agglomeration characteristics of grain production in the HHH were significant. The global spatial correlation characteristics underwent a "decrease–growth–decrease" change process, and the local spatial correlation characteristics demonstrated a concentrated distribution. Specifically, the hot spots were mainly distributed in the central and eastern regions of the HHH, and the cold spots were distributed in the north and southwest. The global and local spatial autocorrelation characteristics showed that the polarization of grain production in local areas has gradually weakened and the spatial difference has gradually decreased in the HHH, which indicates that its agricultural production has gradually shifted in the direction of sustainable development;

The impact of social–economic factors on grain production was constantly strengthened and the influence of climate factors on grain production was gradually weakened. EIA, AFA, and PCI helped to increase the grain yield per unit of sown area in the HHH; however, the effect of the PRE on grain production became weaker as time went on. We adopted the GWR model to prove that the EIA, AFA, and PCI had clear spatial heterogeneity in the intensity and direction of the local scale. The results showed that the EIA had a larger impact on grain production in the HHH compared with other factors, with the percentage of significance at 75.47%.

**Author Contributions:** C.Z. and X.N. developed the main ideas of the study, gathered the data, performed the model construction and estimation, and wrote the manuscript. R.Z. helped to collect and process data. Z.Z. participated in revising the manuscript and proofreading the article. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China [41771190].

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article Crepis vesicaria* **L. subsp.** *taraxacifolia* **Leaves: Nutritional Profile, Phenolic Composition and Biological Properties**

**Sónia Pedreiro 1,2, Sandrine da Ressurreição 3,4 , Maria Lopes 1,2 , Maria Teresa Cruz 1,5 , Teresa Batista 1,6 , Artur Figueirinha 1,2,\* and Fernando Ramos 1,2**


**Abstract:** *Crepis vesicaria* subsp. *taraxacifolia* (Cv) of Asteraceae family is used as food and in traditional medicine. However there are no studies on its nutritional value, phenolic composition and biological activities. In the present work, a nutritional analysis of Cv leaves was performed and its phenolic content and biological properties evaluated. The nutritional profile was achieved by gas chromatography (GC). A 70% ethanolic extract was prepared and characterized by HLPC-PDA-ESI/MS<sup>n</sup> . The quantification of chicoric acid was determined by HPLC-PDA. Subsequently, it was evaluated its antioxidant activity by DPPH, ABTS and FRAP methods. The anti-inflammatory activity and cellular viability was assessed in Raw 264.7 macrophages. On wet weight basis, carbohydrates were the most abundant macronutrients (9.99%), followed by minerals (2.74%) (mainly K, Ca and Na), protein (1.04%) and lipids (0.69%), with a low energetic contribution (175.19 KJ/100 g). The Cv extract is constituted essentially by phenolic acids as caffeic, ferulic and quinic acid derivatives being the major phenolic constituent chicoric acid (130.5 mg/g extract). The extract exhibited antioxidant activity in DPPH, ABTS and FRAP assays and inhibited the nitric oxide (NO) production induced by LPS (IC<sup>50</sup> = 0.428 ± 0.007 mg/mL) without cytotoxicity at all concentrations tested. Conclusions: Given the nutritional and phenolic profile and antioxidant and anti-inflammatory properties, Cv could be a promising useful source of functional food ingredients.

**Keywords:** *Crepis vesicaria* L. subsp. *taraxacifolia*; nutritional value; phenolic profile; chicoric acid; antioxidant; anti-inflammatory

#### **1. Introduction**

The genus *Crepis* belongs to the Asteraceae family comprising about 200 species, and is widely distributed in the Northern Hemisphere, Africa and also in South East Asia [1]. The aerial parts and roots from the plants of the genus *Crepis* are widely used in foods like salads [2], infusions [3], decoctions [4], omelettes, pasta dough and pan-fried [2]. The plants of this genus are also used in traditional medicine to treat jaundice [5], hepatic disorders [6], cardiovascular diseases [7], cough [3,6], catarrh [3], cold [3], diabetes [4], kidney stones [8], eye diseases [6], abdominal colic [6], depurative, blood cleaning, laxative and as a diuretic [2]. It can be used externally in wound healing, bruises and inflammations [8].

**Citation:** Pedreiro, S.; da Ressurreição, S.; Lopes, M.; Cruz, M.T.; Batista, T.; Figueirinha, A.; Ramos, F. *Crepis vesicaria* L. subsp. *taraxacifolia* Leaves: Nutritional Profile, Phenolic Composition and Biological Properties. *Int. J. Environ. Res. Public Health* **2021**, *18*, 151. https://doi.org/10.3390/ijerph18010151

Received: 14 November 2020 Accepted: 23 December 2020 Published: 28 December 2020

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/ licenses/by/4.0/).

Biological properties and phenolic composition have been evaluated for some species. Aerial part and root extracts of *Crepis foetida* L. subsp. *rhoeadifolia* showed antioxidant activity in DPPH and thiobarbituric acid reactive substances (TBARS) assays. In these extracts phenolic compounds to which antioxidant activity has been attributed were identified, namely chlorogenic acid and luteolin in the aerial parts and chlorogenic acid in the roots [7]. A methanolic extract from the flowers of this species showed that chlorogenic acid was the major phenolic compound. This extract presented high antiproliferative activity in HEPG-2, Caco-2, MCF-7 and MCF-10A cells, antioxidant (in DPPH, ABTS, nitric oxide and superoxide radical scavenging assays), anticholinesterase and antityrosinase activities [9]. Aqueous and ethanolic extracts of *Crepis japonica* L. showed antiproliferative activity against leukemia and sarcoma. Moreover, the ethanolic extract presented antiviral activity against respiratory syncytial vírus (anti-RSV). The phenolic content, including hydrolysable tannins may be responsible for this activity [10]. The 3,4 dicaffeoylquinic acid, 3,5-dicaffeoylquinic acid and luteolin-7-*O*-glucoside were isolated from the *Crepis japonica* ethanolic extract. The first two compounds exhibited significant anti-RSV activity. Moreover these three compounds together showed some antibacterial activity against *Vibrio cholerae* and *Vibrio parahaemolyticus* [11]. Besides these activities in this plant, anti-inflammatory, immunosuppressive, antiallergic, antioxidant, analgesic, central nervous system depressant and nematicidal activities were also reported [12].

It is known that reactive oxygen species (ROS), when unregulated, are related to several pathologies such as inflammation, through NF-kB signaling pathway activation [13]. The NF-kB transcription factor is present in the cytoplasm and is in its inactive state due to its association with the inhibitor complex of nuclear factor kappa B kinase (IKK). When this kinase is activated, NF-kB is released and enters the nucleus by activating the transcription of a variety of genes that participate in inflammatory and immune responses [14], such as interleukines IL-1β, IL-6 and IL-8, tumor necrosis factor-α (TNF-α) [15], prostaglandins [16], chemokine ligand 5 (CCL5), transcription of inducible nitric oxide synthase (iNOS), leading to the production of nitric oxide (NO) [17], and cyclooxygenase-2 (COX-2) [18], among others. Phenolic compounds such as phenolic acids, flavonoids and tannins have been identified as good radical scavengers [17]. The mechanisms by which they act in the radical scavenging are involved in signaling pathways of inflammation activation [13]. Phenolic acids presented anti-inflammatory activity acting mainly at the level of the proteasome, inhibiting this and, also, the activation of the NF-kB, since it maintains the phosphorylation levels of IkBα. These mechanisms were attributed to the action of chlorogenic acid, since this was identified as main phenolic acid together with *p*-coumaric acid derivatives [17]. These derivatives have been found to have the ability to inhibit iNOS-dependent NF-kB and COX-2 expression [19]. Phenolic acids from *Lippia* genus inhibited the carrageenan-induced pro-inflammatory cytokine production, namely IL-1β, IL-33 and TNF-α and consequently suppressed the NF-kB activation [16]. Polyphenols from *Ilex latifolia* Thunb. ethanolic extract showed high antioxidant and anti-inflammatory activities, through decreasing the production of NO, COX-2 and pro-inflammatory cytokines via inhibitions of MAPKs, namely ERK and JNK, and NF-kB activation [20]. In LPS-induced acute lung injury rats model, chlorogenic acid decreased the activity of iNOS and suppressing the activation of NF-kB [16]. Others phenolic compounds decreased TLR-4 upregulation, NOX activation and NF-kB activation in LPS-induced renal inflammation rat model [16]. The antioxidant and anti-inflammatory activities of rice bran (RB) phenolic compounds were evaluated in human umbilical vein endothelial cells (HUVECs) with induced oxidative stress. The RB extract regulated antioxidant genes, namely Nrf2, NQO1, HO1 and NOX4, as well antiinflammatory genes (ICAM1, eNOS, CD39 and CD73). This activities were attributed to synergistic interactions between phenolic acids including *p*-coumaric acid, vanillic acid, caffeic acid, ferulic and syringic acid [21]. These studies support the correlation between oxidative stress and inflammation as well the biological effects of phenolic compounds on these.

There are no studies on phytochemical composition and biological activities of *Crepis vesicaria* L. subsp. *taraxacifolia* (Cv). Commonly known as beaked hawk's beard [2], this plant is used traditionally in foods and the treatment of diverse ailments. The cooking water of Cv young leaves are traditionally used as depurative, blood cleaning, diuretic and laxative [2]. Therefore it's important to study this phytoconstituents and to evaluate its health impact. In the present work, there was evaluated the phenolic and nutritional composition of Cv, as well assessed the antioxidant and anti-inflammatory activities.

#### **2. Material and Methods**

#### *2.1. Plant Material and Extract Preparation*

*Crepis vesicaria* L. subsp. *taraxacifolia* (Cv) plant was collected and identified by J. Paiva (Botany Department, University of Coimbra, Coimbra, Portugal). A voucher specimen (A. Figueirinha, 175) was deposited in the herbarium of the University of Coimbra, Faculty of Pharmacy. The leaves, dried in a circulating air drying oven, were milled in a knife mill (KSM 2, BRAUN, Frankfurt, Germany), avoiding the overheating of the sample, and sieved through a 60 mesh sieve. Subsequently, extracts were prepared from the powdered material with different solvents in a proportion of 1:100 (*w*/*v*). In order to improve the extraction of more active compounds, several extractions of Cv leaves (10 mg of dry plant/mL) were made with ethanol/water in various grades: 10%, 30%, 55%, 70% and 100% EtOH. The results of three independent assays showed that the reduction percentage of DPPH radical for the different extracts at 0.33 mg/mL was: 100% EtOH (28.22 ± 0.1575%) < 15% EtOH (66.92 ± 1.083%) < 30% EtOH (91.83 ± 0.602%) < 55% EtOH (90.96 ± 2.205%) < 70% EtOH extract (91.87 ± 2.066%). Thus, a leaves infusion of Cv (ICv) prepared according to its ethomedicinal uses, and 70% ethanol extracts were obtained. These extracts were filtered under vacuum, concentrated in a rotavapor at 40 ◦C, frozen, freeze-dried and kept at −20 ◦C in the dark until use. In the leaves 70% ethanol Cv extract (Cv 70% EtOH), a yield of 22.85% of dry plant was obtained. Relatively to ICv the yield of dry plant obtained was 30.1%. Both extracts were rich in soluble phenolics.

#### *2.2. Chemical Characterization*

#### 2.2.1. Nutrient Composition Analysis

Proximate composition parameters were measured according to the international standards methods of Official Methods of Analysis of AOAC International [22], except neutral detergent fiber (NDF), acid detergent fiber (ADF) and acid detergent lignin (ADL) for calculation of cellulose, hemicellulose and lenhin, respectively [23]. Moisture evaluation was performed by oven drying sample at 105 ◦C until constant weight. Protein was determined by Kjeldahl method, using a protein conversion factor of 6.25. Lipids were gravimetrically quantified after a continuous extraction process in a Soxhlet apparatus by diethyl ether. Fatty acids were analysed by gas phase chromatography (GC-FID) of fatty acid methyl esters, and the quantification was performed by Supelco standards (Sigma-Aldrich, St Louis, MO, USA). The total dietary fiber, soluble and insoluble dietary fiber contents were determined using the Supelco enzyme kit TDF100A (Sigma-Aldrich). Crude fiber was gravimetrically quantified after chemical digestion and solubilisation of other materials. The fiber residue weight was then corrected for ash content. Ash was obtained by incineration of all organic matter of the sample in a muffle furnace at 550 ◦C. The Nitrogen-free extractives were estimated, considering the following equation: Nitrogen-free extractives = 100 − (moisture + ash + lipids + protein + crude fiber). The total carbohydrates were estimated, considering the following equation: Total carbohydrates = 100 − (moisture + ash + lipids + protein). Quantification and sugars were performed by High Performance Liquid Chromatography with refractive index detection (HPLC-RI). The separation column used was HC-75 Ca2+ <sup>305</sup> <sup>×</sup> 7.8 mm (Hamilton, Energy Way Reno, NV, USA) with ultrapure water with traces of sodium azide mobile phase, at a flow rate of 0.6 mL/min, at 80 ◦C. The quantification was performed by BioUltra standards (Sigma-Aldrich). The available carbohydrates were estimated, based on the following equation:

Available Carbohydrates = 100 − (moisture + ash + lipids + protein + dietary fiber). Quantification of sugars were performed by high performance liquid chromatography (HPLC), using BioUltra standards (Sigma-Aldrich). Energy values are expressed in Kcal and KJ/100 g and were calculated according to Regulation (EU) n◦ 1169/2011 of the European Parliament and of the Council of 25 October 2011. Minerals were determined by flame atomic absorption spectrometry (FAAS), with the exception of cadmium and lead traces, which were determined by graphite furnace atomic absorption spectrometry (GFAAS). Mercury traces was analysed by an AMA254 Mercury Analyzer (Leco, St Joseph, MI, USA) and phosphorus, by spectrophotometry.

#### 2.2.2. Phenolic Profile HPLC-PDA-ESI/MS<sup>n</sup>

The phenolic profile of Cv (ethanol 70% extract) was carried out on a liquid chromatograph (Finnigan Surveyor, THERMO, Waltham, MA, USA) with a Spherisorb ODS-2 column (150 × 2.1 mm i.d.; particle size, 3 µm; Waters Corp., Milford, MA, USA) and a Spherisorb ODS-2 guard cartridge (10 × 4.6 mm i.d.; particle size, 5 µm; Waters Corp. Milford, MA, USA). The separation occurred at 25 ◦C with a mobile phase constituted by 2% aqueous formic acid (v/v) (A) and methanol (B) in a discontinuous gradient of 5–15% B (0–10 min), 15–25% B (10–15 min), 25–50% B (15–40 min), 50–80% B (40–50 min), followed by an isocratic elution (50–60 min), a gradient 80–100% B (60–65 min) and other isocratic elution for 5 min, at a flow rate of 200 µL/min.

The first detection was done with a PDA detector ((Finnigan Surveyor, THERMO, Waltham, MA, USA)) at a wavelength range 200–400 nm, followed by a second detection using an Linear Ion Trap Mass Spectrometer (LIT-MS) (LTQ XL, Thermo Waltham, MA, USA). Mass spectra were obtained in the negative ion mode. The mass spectrometer acquired three consecutive spectra: full mass (*m*/*z* 125–1500), MS<sup>2</sup> of the most abundant ion in the full mass and MS<sup>3</sup> of the most abundant ion in the MS<sup>2</sup> . Source voltage was 4.5 kV and the capillary temperature and voltage were 250 ◦C and −10 V, respectively. The sheath and auxiliary gas used was nitrogen at 20 Finnigan arbitrary units with helium as collision gas with a normalized energy of 45%. XCALIBUR software (Thermo, Waltham, MA, USA) was used for data treatment.

#### 2.2.3. Quantification by HPLC-PDA

Quantification of L-chicoric acid in Cv 70% EtOH was performed in a chromatograph with a photodiode array (PDA) detector (Gilson Electronics SA, Villiers le Bel, France). The analysis were performed on a Spherisorb S5 ODS-2 column (250 × 4.6 mm i.d., 5 µm) (Waters Milford, MA, USA) with a C18 guard cartridge (30 × 4 mm i.d., 5 µm) (Nucleosil, Macherey-Nagel, Düren, Germany), at 24 ◦C. The mobile phase was constituted by methanol 100% (B) and formic acid 5% (A). The elution was made at a flow rate of 1 mL/min. The gradient used was: 5–15% B (0–10 min), 15–25% B (10–15 min), 25–50% B (15–40 min), 50–80% B (40–50 min) followed by an isocratic elution of 80% B (50–60 min), 80–100% B (60–70 min) and finally, an isocratic elution of 100% B (70–85 min). The volume of the sample injected was 100 µL. The UV-Visible spectra acquisition was performed between 200–600 nm and the chromatographic profiles were recorded at the wavelengths 280 and 320 nm. Data treatment was carried out with Unipoint®, version 2.10 software (Gilson, Middleton, WI, USA).

Quantification of the L-chicoric acid was performed using commercial standard dissolved in methanol (10 to 150 µg/mL) as external standard L-chicoric acid (Sigma Aldrich St. Louis, MO, USA). The chicoric acid present in the Cv 70% EtOH extract was quantified by the absorbance recorded in the chromatogram relative to this standard (330 nm). Three independent injections (100 µL) were performed in duplicate for each sample. The least-squares regression model was used to assess the correlation between peak area and concentration. The detection (LOD) and quantification (LOQ) limits were calculated from the calibration curve. The quantification of the chicoric acid in Cv 70% EtOH extract

(identified first by HPLC-PDA-MS<sup>n</sup> ) was made using the standard calibration curve and the peak area of the compound.

#### *2.3. Antioxidant Activity*

#### 2.3.1. 2,2-Diphenyl-1-Picrylhydrazyl Radical Assay (DPPH)

Free radical-scavenging activity of the infusion and ethanol/water Cv extracts was evaluated using the *DPPH* method previously described [24]. Briefly, aliquots of samples (100 µL) were assessed by their reactivity with methanolic solution of 500 µM *DPPH* (500 µL) in the presence of 100 mM acetate buffer, pH 6.0 (100 µL). The reaction mixtures (300 µL) were kept for 30 min at room temperature, in the dark. The decreases in the absorbance were measured at 517 nm, in a Thermo scientific multiskan FC plate reader. The % of reduction of *DPPH* of the Cv extracts were determined by:

$$\% \text{ reduction of DPPH} = 100 - \frac{\text{Abs sample} - \text{Abs control}}{\text{Abs control}} \tag{1}$$

Posteriorly, the obtained values were plotted in a graph % of DPPH reduction vs. concentration in µg/mL. The IC<sup>50</sup> was interpolated in the graph for the correspondent value of 50% reduction.

The results were expressed as Trolox equivalent (TE), defined as the concentration of the extract with antioxidant capacity equivalent to 1 mM of Trolox solution. This value was obtained interpolating the absorbance of 1 mM Trolox in the graph % of DPPH reduction vs. concentration. All the determinations were performed in triplicate.

#### 2.3.2. Ferric Reducing Antioxidant Power Assay (FRAP)

Ferric reducing ability was evaluated with slight modifications [25]. The FRAP reagent was prepared by mixing 300 mM acetate buffer, 10 mL TPTZ (Sigma–Aldrich St. Louis, MO, USA) in 40 mM HCl and 20 mM FeCl3.6H2O (Merck, Darmstadt, Germany) in the proportion of 10:1:1 (*v*/*v*/*v*). The extract (100 µL) was added to 3 mL of the FRAP reagent. An intense blue color complex was formed when ferric tripyridyl triazine (Fe3+ TPTZ) complex was reduced to ferrous (Fe2+) form. The absorbance was measured at 593 nm, against a reagent blank, after incubation at room temperature for 6 min. The results were expressed as trolox equivalent (TE) values obtained using a calibration curve for Trolox (31.25–1000 mM). All the determinations were performed in triplicate.

#### 2.3.3. 2,20 -Azinobis-(3-ethylbenzothiazoline-6-sulfonate) Assay (pH = 7) (ABTS)

The ABTS assay described by [26], consisted in generating the ABTS•+ radical by the oxidation of ABTS (7 mM) with potassium persulphate (2.45 mM) (Merck) in water. After 12–16 h in dark and at room temperature, this solution was diluted with phosphate buffered saline (PBS) at pH 7 to give an absorbance of 0.7 ± 0.02 at 734 nm. The extract (50 µL) was mixed with 2 mL of the ABTS + solution and vortexed for 10 s. After 4 min of reaction, the absorbance was measured at 734 nm. The IC<sup>50</sup> value was interpolated in a graph % of ABTS reduction vs. concentration in µg/mL for the correspondent value of 50% reduction. The results were expressed as TE, obtained interpolating the absorbance of 1 mM trolox in the graph % of ABTS reduction vs. concentration. Three independent experiments in triplicate were performed for each of the assayed extracts.

#### *2.4. Anti-Inflammatory Activity Evaluation*

#### 2.4.1. Nitrite Production by Griess Assay

Raw 264.7, a mouse leukemic monocyte macrophage cell line from American Type Culture Collection (Manassas, VA, USA), and kindly supplied by Dr. Otília Vieira (Center for Neurosciences and Cell Biology, University of Coimbra, Portugal), was cultured in Iscove's Modified Dulbecco's Eagle medium supplemented with 10% non-inactivated fetal bovine serum, 100 U/mL penicillin, and 100 µg/mL streptomycin at 37 ◦C in a humidified atmosphere of 95% air and 5% CO2. The cells were monitored to detect any morphological

change. For the experiments, the cells were plated (0.6 <sup>×</sup> <sup>10</sup><sup>5</sup> cells/well) with culture medium and allowed to stabilize for 12 h. Then the cells were incubated during 24 h at 37 ◦C in culture medium (control) or stimulated with 1 µg/mL of bacteria lipopolysaccharide (LPS) (Sigma) with or without different concentrations of extract (0.1–2.0 mg/mL).

The anti-inflammatory activity was determined by the nitric oxide production, measured indirectly by the accumulation of nitrite in the supernatant through a colorimetric assay with Griess reagent [0.1% (*m*/*v*) of N-(1-naphthyl)-ethylenediamine dihydrocloride and 1% (*m*/*v*) of sulfanilamide with 5% of phosphoric acid] [27]. To perform the assay, it was used 100 µL of the supernatant and 100 µL of Griess's reagent and then stored away from light during 30 min. The absorbance at 550 nm was measured in an automated plate reader (Synergy HT, BioTek Instruments SAS, Colmar, France). Culture medium was used as blank and nitrite concentration was determined from a regression analysis using serial dilutions of sodium nitrite as standard.

#### 2.4.2. Assessment of Cell Viability by Resazurin Assay

In order to evaluate the cytotoxicity it was performed the resazurin assay [28]. After the incubation with the samples, the cells were incubated with 100 µL of a resazurin solution (10 µM in culture medium) during 2 h at 37 ◦C in a humid atmosphere with 5% CO2/95% air. Quantification of resorufin was performed using a plate reader (Synergy HT, BioTek, Instruments SAS, Colmar, France) at 570 nm, with an optical filter for 620 nm.

#### *2.5. Statistical Analysis*

All samples were analysed, at least, in triplicates and the results were expressed as mean ± standard deviation (SD). To calculate the IC<sup>50</sup> values for the anti-inflammatory activity, the linearization of the dose-response curve was performed as described by Chou [29].

The statistical analysis of the cellular viability and anti-inflammatory activity was performed in GraphPad Prism program (version 5.02, GraphPad Software, San Diego, CA, USA). For the comparison between treatment conditions and control it was used two-sided unpaired *t*-test. To evaluate the effect of different treatments to LPS-stimulated cells it was performed One-way ANOVA followed by Bonferroni's test. The limit of significance was set at \*\*\* *p* < 0.001.

#### **3. Results and Discussion**

#### *3.1. Nutrient Composition of C. vesicaria*

The knowledge of the nutritional properties of wild plants is crucial to assess their suitability for human consumption. In this study, the nutritional profile of *Crepis vesicaria* subsp. *taraxacifolia* leaves was analyzed.

#### 3.1.1. Nutritional Analysis of *Crepis vesicaria* subsp. *taraxacifolia* Leaves

The nutritive content of Cv leaves was determined (Table 1). The proximate composition revealed high moisture content, even though all foods contain water; those with a higher content are more prone to the rapid occurrence of microbial spoilage phenomena, enzymatic degradation and other moisture-dependent chemical deterioration reactions. Therefore, precautions should be considered to prevent rapid deterioration during storage, such as drying or freezing.


**Table 1.** Nutritive content of *Crepis vesicaria* subsp. *taraxacifolia* leaves (mean ± SD; *n* = 3).

Total carbohydrates, calculated by difference, were the most abundant macronutrients (9.99 g/100 g wet weight (*w*/*w*)), followed by ash, protein and lipids. Carbohydrates play a major role in human diet. They are the main source of energy, and also help to maintain glycemic homeostasis and gastrointestinal integrity, among other functions. A healthy adult diet should include about 130 g of carbohydrates per day [30]. Cv leaves contain an important amount of carbohydrates, which is in line with what has been reported for other wild Asteraceae plants traditionally consumed in the Mediterranean region, such as *Taraxacum obovatum*, *Chondrilla juncea*, *Sonchus oleraceus*, *Cichorium intybus*, *Scolymus hispanicus* and *Silybum marianum* [31]. An important fraction of the total carbohydrates content in Cv leaves is fiber. In this study, different fiber measurement methods were used and the results showed that the chosen method has an impact on the values observed for different fiber parameters. Weende's crude fiber analysis determines cellulose, lignin and some hemicellulose, pectin, gums and mucilages. The acid detergent lignin (ADL) measures lignin, the acid detergent fiber (ADF) determines cellulose and lignin, and the neutral detergent fiber (NDF) consists mainly in the measurement of cellulose, hemicelluloses and lignin [32]. Regardless of the method, the results reveal that Cv leaves are an interesting dietary fiber source, with insoluble dietary fiber being the major fraction. It is well established that the daily consumption of about 25–30 g of fiber, for an adult, markedly reduces the risk of cardiovascular and digestive diseases [30]. Also, Cv leaves may contain insoluble-bound phenolics present in the cell wall plant components. These insoluble-bound form can contribute for to protection of cardiovascular health [33]. Thus, the use of this plant, either individually or added to other foods, may contribute to a desired increase in fiber intake with the associated health benefits.

With regard to the available carbohydrates, the estimated value was 5.75 g/100 g (*w*/*w*). The total sugars content found was 3.76 g/100 g (*w*/*w*), with maltose as the main sugar (2.47 g/100 g, *w*/*w*), followed by fructose and glucose. Protein makes up 1.04 g/100 g, *w*/*w* of Cv leaves. This value is considerably lower than that reported by Barnett and Crawford [34]. Variations in protein levels may be due to differences between species, environmental and climatic factors, or a mixture of both.

#### 3.1.2. Lipid and Fatty Acids Composition of *Crepis vesicaria* subsp. *taraxacifolia* Leaves

According to Table 2, the lipid content was moderate, 0.69 g/100 g, *w*/*w* (4.78 g/100 g dry weight (dw)), higher than that reported for *C. Juncea*, the highest lipid content Asteraceae (0.79 g/100 g, *w*/*w*) analyzed in the study of García-Herrera [31]. The fatty acid profile of Cv leaves showed a predominance of polyunsaturated fatty acids (PUFA) (402.84 mg/100 g, *w*/*w*), mainly comprised by α-linolenic acid (343.24 mg/100 g). Total

saturated fatty acids (SFA) concentration was 159.82 mg/100 g, *w*/*w*, with the main constituent being palmitic acid (108.75 mg/100 g, *w*/*w*). For a nutritional "good quality", including beneficial effects in terms of cardiovascular risk reduction, the PUFA/SFA ratio should be > 0.45, whilst n-3/n-6 fatty acids ratio should be > 4 [35]. In the present study, the PUFA/SFA ratio was 2.52 and the n-3/n-6 fatty acids ratio was 5.76. The presence of considerable amounts of oleic acid (60.49 mg/100 g, *w*/*w*) should also be highlighted, given the beneficial properties that have been attributed to it in the context of the immunomodulation, prevention and treatment of several pathologies such as cancer, cardiovascular and autoimmune diseases, and metabolic disturbances [36].

**Table 2.** Lipid and fatty acids composition of *Crepis vesicaria* subsp. *taraxacifolia* leaves (mean ± SD; *n* = 3).


#### 3.1.3. Minerals and Heavy Metal Composition of Cv Leaves

Given the results in Table 3, *Crepis vesicaria* L. subsp. *taraxacifolia* leaves exhibited moderate levels of ash (2.74 g/100 g, *w*/*w*). This value is within the recommended range for human consumption and reveals considerable mineral richness, which is corroborated by studies on similar species, such as *C. commutata* and *C. vesicaria* [37]. The mineral fraction is an aspect of greater relevance in the use of edible plants in human nutrition. The inappropriate intake of minerals (macrominerals and trace minerals) is the cause of multiple degenerative and chronic diseases. Calcium (Ca), phosphorous (P), magnesium (Mg), sodium (Na), potassium (K) and iron (Fe) are essential elements and their intake is necessary at mg/kg level to keep the human body healthy. Zinc (Zn), copper (Cu), manganese (Mn), chromium (Cr), and nickel (Ni) are required at trace levels in the diet [38]. Concerning the macrominerals composition of Cv leaves, K, Ca and Na were the most abundant (591.29 mg/100 g, *w*/*w*; 309.93 mg/100 g, *w*/*w*; 76.78 mg/100 g, *w*/*w*, respectively). The macromineral profile found was identical to that reported for *C. vesicaria* (K > Ca > Na > P > Mg), but different from that observed in *C. commutata* (K > P > Na > Mg > Ca) [37]. When compared to other wild Asteraceae, such as *S. hispanicus*, K and Ca contents are lower, but Cv can still be considered a remarkable source of these minerals, better than many conventional vegetables [31]. Zn and Mn were found as major trace minerals (0.86 mg/100 g, *w*/*w* and 0.83 mg/100 g, *w*/*w*, respectively). The most abundant trace minerals in *C. vesicaria* and *C. commutate* were Fe and Mn. *Crepis* spp. seem to be a good source of Mn. The contaminants cadmium (Cd), lead (Pb) and mercury (Hg) were detected. These toxic metallic elements can induce damage to multiple organs and have carcinogenic effects [39]. Pb levels are below the maximum values legislated, 0.30 mg/kg, *w*/*w*. However, the concentration of Cd coincides with the maximum level of contamination that is considered safe, 0.2 mg/kg, *w*/*w* [40]. Overall, the results indicate that, when located in polluted areas, these plants can accumulate toxic metals in concentrations that may represent a risk to the consumer's health.


**Table 3.** Minerals and heavy metal composition of *Crepis vesicaria* subsp. *taraxacifolia* leaves (mean ± SD; *n* = 3).

#### *3.2. Screening for Antioxidant/Scavenging Activity*

The ability of ROS to activate the inflammation signaling pathway, through activation of pro-inflammatory cytokines is well known. The literature describes that colorimetric methods to assess antioxidant activity like DPPH and ABTS, are a good tool to select the extracts more promising [41]. Also, it was reported that phenolic extracts bearing higher radical scavenging towards DPPH and ABTS, present higher inhibition of NF-kB activation mediated by ROS [41]. Given the correlation between antiradical activity and inhibition of the NF-kB signaling pathway, the antioxidant activity of the extracts was screened using the DPPH and ABTS colorimetric methods. Based on the results obtained, it was chosen the extract that demonstrated the greatest activity in these tests.

The infusion (10 mg of dry plant/mL) was screened for antioxidant activity as it is the form of use in traditional medicine. However, the percentage of reduction observed was 20.27%. As the Cv 70% EtOH extract was the most active extract it was lyophilized (previously described in Material and Methods) and characterized relatively to its phenolic profile, and antioxidant and anti-inflammatory activities.

Regarding antioxidant activity, the infusion presented an IC<sup>50</sup> of 103.22 ± 5.61 µg/mL and a TEAC of 441.980 ± 0.058 mg/mL. The IC<sup>50</sup> of Cv 70% EtOH was 26.20 ± 1.86 µg/mL and TEAC of 111.980 ± 0.041 mg/mL meaning that this extract is more active than the infusion. Subsequently, the antioxidant activity of the 70% ethanol extract was assessed by FRAP and ABTS methods (Table 4). The ABTS and DPPH methods are based on electron and H atom transfer while the FRAP is based on electron transfer reaction [42]. Attending to the results, the Cv extract present reducing power besides their ability in scavenging free radicals. The results shown that the Cv 70% EtOH extract has a good radical-scavenging activity and antioxidant activity.

**Table 4.** Antioxidant activity of ethanolic extract (Cv 70% EtOH) from *Crepis vesicaria* L. subsp. *taraxacifolia*.


\* TE (Trolox Equivalent): Amount of the samples (µg/mL) that has the same anti-radical activity of Trolox 1 mM. The results are expressed as mean ± SD of three independent experiments.

*taraxacifolia*.

#### *3.3. Phenolic Profile of 70% Ethanolic Extract from Crepis vesicaria subsp. taraxacifolia 3.3. Phenolic Profile of 70% Ethanolic Extract from Crepis vesicaria subsp. taraxacifolia*

DPPH● 26.20 ± 1.86 111.980 ± 0.041

(pH = 7) 18.92 ± 2.24 21.670 ± 0.012

Trolox 1 mM. The results are expressed as mean ± SD of three independent experiments.

**Table 4.** Antioxidant activity of ethanolic extract (Cv 70% EtOH) from *Crepis vesicaria* L. subsp.

**IC<sup>50</sup> (µg/mL) TE \***

*Int. J. Environ. Res. Public Health* **2021**, *18*, x 10 of 16

Based on the given results relatively to the antioxidant activity, the 70% EtOH from Cv extract is the most active. Therefore, the phenolic profile by HPLC-PDA-MS<sup>n</sup> of this extract was assessed (Figure 1). Based on the given results relatively to the antioxidant activity, the 70% EtOH from Cv extract is the most active. Therefore, the phenolic profile by HPLC-PDA-MS<sup>n</sup> of this extract was assessed (Figure 1).

\* TE (Trolox Equivalent): Amount of the samples (µg/mL) that has the same anti-radical activity of

Figure 1. HPLC-PDA-ESI/MS<sup>n</sup> profile of 70% ethanol extract from Cv recorded at 280 mn. It was used the gradient 2 described in Material and Methods section. (The chromatogram of the extract is not shown up to 40 min as no further compounds were eluted after this time period. Peaks 1–6 identification is showed in Table 5). **Figure 1.** HPLC-PDA-ESI/MS<sup>n</sup> profile of 70% ethanol extract from Cv recorded at 280 mn. It was used the gradient 2 described in Material and Methods section. (The chromatogram of the extract is not shown up to 40 min as no further compounds were eluted after this time period. Peaks 1–6 identification is showed in Table 5).

According to the absorption spectra, the 70% EtOH Cv extract is mainly composed of phenolic acids, generally presenting a shoulder at 295 nm and a maximum wavelength of 330 nm (Table 5), indicating to be caffeic or ferulic acids derivatives [43]. According to the absorption spectra, the 70% EtOH Cv extract is mainly composed of phenolic acids, generally presenting a shoulder at 295 nm and a maximum wavelength of 330 nm (Table 5), indicating to be caffeic or ferulic acids derivatives [43].


**Table 5.** Compounds identified in Cv 70% ethanol extract by HPLC-PDA-ESI/MS <sup>n</sup> . **Table 5.** Compounds identified in Cv 70% ethanol extract by HPLC-PDA-ESI/MS <sup>n</sup> .

Identification based on the UV-Vis spectra, molecular weight and fragmentation patterns, which are according to authors referred.

In an attempt to identify the compounds of this extract, HPLC-PDA-ESI/MS<sup>n</sup> was performed. The results (Table 5) showed that the extract consisted mainly of phenolic acids, namely caffeic and ferulic acid derivatives as well as chicoric acid isomers. The chicoric acid was identified as the main compound of the Cv 70% ethanol extract.

Compound **1**. MS analysis showed a molecular ion [M-H]- at *m*/*z* 179 and a fragmentation pattern typical of caffeic acid. MS<sup>2</sup> data presented fragments at *m*/*z* 135 indicating a

decarboxylated caffeic acid moiety [(M-H-CO2] - . The compound 1 was tentatively identified as caffeic acid [44].

Compound **2**. This compound presents a molecular ion [M-H]- at *m*/*z* 191. MS<sup>2</sup> most abundant fragments are *m*/*z* 173 indicating a dehydrated quinic acid moiety [M-H-H2O]- . This compound was tentatively identified as quinic acid [44].

Compounds **3**, **4** and **5**. The molecular ion [M-H]- occurs at *m*/*z* 473. The MS<sup>2</sup> presents a fragment at *m*/*z* 311, indicating the presence of deprotonated caftaric acid [M-H-C13H12O9] - and *m*/*z* 293 corresponding to the neutral loss of caffeic acid [M-H-C9H8O4] - . MS<sup>3</sup> profiles have a fragment at *m*/*z* 149 corresponding to the tartaric acid and at *m*/*z* 179 corresponding to a deprotonated molecule of caffeic acid [M-H-C9H8O4] - . Based on this fragmentation pattern, these compounds were tentatively identified as chicoric acid isomers [47]. Accordingly with the literature, the most abundant chicoric acid isomer is L-chicoric acid [47]. The quantification by HPLC-PDA of this isomer was performed using a standard. The peak of the L-chicoric acid standard has approximately the same retention time of peak 5. Therefore, peak 5 was tentatively identified as L-chicoric acid.

Compound **6**. This compound has a molecular ion [M-H]- at *m*/*z* 487. MS<sup>2</sup> showed fragments at *m*/*z* 325 (loss of 162 Da) indicating the loss of a hexose [M-H-C6H12O6] - and at 307 that indicates the subsequent loss of water [M-H-C6H12O6-H2O]- . The MS<sup>3</sup> fragment *m*/*z* 193 indicates the presence of ferulic acid [48], that can be probably associated to a hexosylpentosyl residue. All the data suggest that compound **6** was tentatively identified as feruloyl hexosylpentoside [46].

#### *3.4. Quantification of Chicoric Acid*

The major constituents in *Crepis vesicaria* subsp. *taraxacifolia* were chicoric acid derivatives and the evaluated activities were attributed to these compounds. Therefore, the Lchicoric acid was quantified by HPLC-PDA. The equation of calibration curve of L-chicoric acid was y = 4011236.2307 × −36474316.1324 (r<sup>2</sup> = 0.99). Based on this equation, the concentration of L-chicoric acid in the Cv 70% EtOH extract was 130.5 ± 4.2 mg/g extract of Cv 70% EtOH. The LOD and LOQ were 19.74 ± 3.33 mg/g extract and 44.58 ± 2.96 mg/g extract, respectively.

#### *3.5. Assessment of Cell Viability of the Cv 70% EtOH Extract Int. J. Environ. Res. Public Health* **2021**, *18*, x 12 of 16

The citotoxicity of the Cv 70% EtOH extract in macrophages was evaluated. The results (Figure 2) showed that none of the tested concentrations were cytotoxic.

**Figure 2.** Effect of *Crepis vesicaria* subsp. *taraxacifolia* ethanolic extract on macrophages cell viability (RAW 264.7 cells). Each result represents the mean ± SD (minimum of three independent assays, performed in duplicate). The statistical tests were performed with *p* < 0.05 compared to control. **Figure 2.** Effect of *Crepis vesicaria* subsp. *taraxacifolia* ethanolic extract on macrophages cell viability (RAW 264.7 cells). Each result represents the mean ± SD (minimum of three independent assays, performed in duplicate). The statistical tests were performed with *p* < 0.05 compared to control.

The cytotoxicity of chicoric acid in macrophages (Raw 264.7 cells) has also been was

The results showed that Cv 70% ethanol extract inhibited NO production in a dosedependent manner (Figure 3) and the IC<sup>50</sup> was 0.428 ± 0.00669 mg/mL. Cv extract is little

are few cell viability studies in normal cells with Cv extracts. However, some researchers studied the effects of a methanol extract of Cv flowers on tumor (HEPG-2, Caco-2 and MCF-7) and non-tumor (MCF-10A) cells. The Cv extract was not cytotoxic to the nontumor cell line and cytotoxic to the tumor lines and therefore had some selectivity over

*3.6. Antioxidant and Anti-Inflammatory Activity of the Cv 70% Ethanol Extract*

active in inflammation compared to the results obtained in antioxidant activity.

tumor cells [9].

**Control**

**0**

**50**

**100**

**Nitrite Production (%)**

**150**

**LPS**

**0.1 mg/mL**

**0.25 mg/mL**

**0.5 mg/mL**

assays). \*\*\* *p* < 0.001 compared with the LPS group.

\*\*\* \*\*\*

**0.75 mg/mL**

**[Cv 70% EtOH] (mg/mL)**

\*\*\*\*\*\*

**1 mg/mL**

**1.25 mg/mL**

**1.5 mg/mL**

\*\*\*\*\*\*

**2 mg/mL**

Figure 3**.** Effect of *Crepis vesicaria* subsp*. taraxacifolia* ethanolic extract on NO production in macrophages (RAW 264.7 cells). Each result represents the mean ± SD (minimum of three independent

\*\*\*

\*\*\*

**Control**

**0**

**50**

**100**

**Cellular viability (%)**

**150**

**0.1 mg/mL**

**0.25 mg/mL**

**0.5 mg/mL**

**0.75 mg/mL**

**[Cv 70% EtOH] (mg/mL)**

**1 mg/mL**

**1.25 mg/mL**

**1.5 mg/mL**

**2 mg/mL**

**Figure 2.** Effect of *Crepis vesicaria* subsp. *taraxacifolia* ethanolic extract on macrophages cell viability

The cytotoxicity of chicoric acid in macrophages (Raw 264.7 cells) has also been was tested. The studies performed have shown that this compound is not cytotoxic [49]. There are few cell viability studies in normal cells with Cv extracts. However, some researchers studied the effects of a methanol extract of Cv flowers on tumor (HEPG-2, Caco-2 and MCF-7) and non-tumor (MCF-10A) cells. The Cv extract was not cytotoxic to the non-tumor cell line and cytotoxic to the tumor lines and therefore had some selectivity over tumor cells [9]. tested. The studies performed have shown that this compound is not cytotoxic [49]. There are few cell viability studies in normal cells with Cv extracts. However, some researchers studied the effects of a methanol extract of Cv flowers on tumor (HEPG-2, Caco-2 and MCF-7) and non-tumor (MCF-10A) cells. The Cv extract was not cytotoxic to the nontumor cell line and cytotoxic to the tumor lines and therefore had some selectivity over tumor cells [9].

The cytotoxicity of chicoric acid in macrophages (Raw 264.7 cells) has also been was

#### *3.6. Antioxidant and Anti-Inflammatory Activity of the Cv 70% Ethanol Extract 3.6. Antioxidant and Anti-Inflammatory Activity of the Cv 70% Ethanol Extract* The results showed that Cv 70% ethanol extract inhibited NO production in a dose-

*Int. J. Environ. Res. Public Health* **2021**, *18*, x 12 of 16

The results showed that Cv 70% ethanol extract inhibited NO production in a dosedependent manner (Figure 3) and the IC<sup>50</sup> was 0.428 ± 0.00669 mg/mL. Cv extract is little active in inflammation compared to the results obtained in antioxidant activity. dependent manner (Figure 3) and the IC<sup>50</sup> was 0.428 ± 0.00669 mg/mL. Cv extract is little active in inflammation compared to the results obtained in antioxidant activity.

Figure 3**.** Effect of *Crepis vesicaria* subsp*. taraxacifolia* ethanolic extract on NO production in macrophages (RAW 264.7 cells). Each result represents the mean ± SD (minimum of three independent assays). \*\*\* *p* < 0.001 compared with the LPS group. **Figure 3.** Effect of *Crepis vesicaria* subsp. *taraxacifolia* ethanolic extract on NO production in macrophages (RAW 264.7 cells). Each result represents the mean ± SD (minimum of three independent assays). \*\*\* *p* < 0.001 compared with the LPS group.

Reactive oxygen species are involved in various pathologies including inflammation. Cv extract was shown to have a high antioxidant activity. According to the characterization of the extract by HPLC-PDA-ESI/MS<sup>n</sup> , the extract has phenolic acids namely caffeic and ferulic acid derivatives as well as chicoric acid isomers. The mechanisms involved in the antioxidant activity by which phenolic compounds act are based on: ability to chelate metals, such as copper and iron, that participate in the Fenton reaction generating hydroxyl radicals; interrupt signaling pathways triggered by free radicals; interfere with enzyme activity [50]. It is known that the antioxidant activity of phenolic compounds is directly related to the number of hydroxyl groups. Chicoric acid was identified as the major compound present in the extract. This compound has two caffeic acid moieties which are responsible for the high activity observed relatively to caffeic acid [51]. Some authors relate high molecular weight phenolic compounds, such as chicoric acid, with antioxidant activity [52]. Therefore, the observed antioxidant activity by H transfer and electron transfer reaction can be mostly attributed to chicoric acid.

Other plants of the genus *Crepis* are reported to have various biological activities including anti-inflammatory activity [6]. Chicoric acid has been reported to have that activity inhibiting activated immune cells, nitric oxide synthase and cyclooxygenase-2 through its inhibitory effects on nuclear factor NF-κB and TNF-α [53–55]. However, the results of Cv 70% EtOH extract weren't satisfactory in the anti-inflammatory activity. This fact can be due to the extract matrix [53] or antagonistic interactions between the matrix compound [56]. Moreover, this extract have chicoric acid isomers, and the observed

activity can be due to this isomers. Further studies are needed to understand the inherent mechanisms of these compounds in anti-inflammatory activity.

#### **4. Conclusions**

The 70% ethanol *Crepis vesicaria* subsp. *taraxacifolia* leaves extract presents antioxidant and anti-inflammatory activities. Besides its biological activities, the Cv leaves extract demonstrated to have a high content of lipids and fatty acids. Given the observed antioxidant activity, Cv extract may be used as a functional food to prevent oxidative stress and its associated pathologies. In fact, Cv leaves displayed an interesting nutritional composition, with a low energetic contribution of 41.84 kcal/100 g, *w*/*w*. Some potential uses for this plant's leaves may be the development of new additives for human and/or animal consumption or food supplements that contribute to a balanced diet. In a context of food insecurity, their incorporation as an ingredient in recipes may also be of interest to increase their nutritional and functional value. From another perspective, to increase the consumption of this abundant and under-exploited plant, it is important to investigate its nutritional value and antioxidant and anti-inflammatory properties, but also to ensure that its consumption is safe, i.e., without neglecting the risk of contamination by toxic metals.

**Author Contributions:** S.P. contributed to the methodology, formal analysis and original draft preparation. A.F., T.B., M.T.C. and F.R. contributed to the project's administration, conceptualization and review and editing. S.d.R. and M.L. participated in the nutritional analysis and its interpretation. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by PT national funds (FCT/MCTES, Fundação para a Ciência e Tecnologia and Ministério da Ciência, Tecnologia e Ensino Superior) through grant UIDB/50006/2020 and by Programa de Cooperación Interreg V-A España–Portugal (POCTEP) 2014–2020 (project 0377\_IBERPHENOL\_6\_E).

**Institutional Review Board Statement:** Not applicable. This work didn't performed in vivo studies in animals or humans.

**Informed Consent Statement:** Not applicable. This work didn't performed human's assays.

**Data Availability Statement:** All data is available based on "MDPI Research Data Policies" at https://www.mdpi.com/ethics.

**Acknowledgments:** We are grateful to Laboratory of Mass Spectrometry (LEM) of UC Node integrated in the National Mass Spectrometry Network (RNEM) of Portugal, for MS analyses.

**Conflicts of Interest:** The authors state no conflict of interest.

#### **Abbreviations**

GC-FID: gas chromatography-flame ionization detector; HPLC-PDA-ESI-MS<sup>n</sup> : high performance liquid chromatography-photodiode array-electrospray ionization-mass spectrometry; HPLC-PDA: high performance liquid chromatography-photodiode array; DPPH: 2,2-diphenyl-1-picrylhydrazyl radical; ABTS: 2,20 -azinobis-(3-ethylbenzothiazoline-6-sulfonate); FRAP: ferric reduction antioxidant power.

#### **References**


### *Article* **Food Insecurity among Low-Income Food Handlers: A Nationwide Study in Brazilian Community Restaurants**

**Ingrid C. Fideles <sup>1</sup> , Rita de Cassia Coelho de Almeida Akutsu 2,\* , Rosemary da Rocha Fonseca Barroso <sup>1</sup> , Jamacy Costa-Souza <sup>1</sup> , Renata Puppin Zandonadi <sup>2</sup> , António Raposo 3,\* and Raquel Braz Assunção Botelho <sup>2</sup>**


**Abstract:** This study aims to evaluate food insecurity (FI) among Brazilian Community restaurant food handlers and its associated factors. This cross-sectional study was performed with a representative sample of 471 food handlers working in community restaurants (CR) from all Brazilian regions. Participants are mostly female (62.2%), ≤40 years old (67.7%), with a partner (52.0%), and with up to eight years of education (54.1%). Predictors of participants' socioeconomic status and CR geographic location are associated with the household food insecurity categories (*p* < 0.05). The predictors of socioeconomic conditions are associated with mild and moderate/severe FI category. Workers with less education are twice as likely to belong to the category with the highest FI severity. Lower per capita household income increased the chances of belonging to the mild insecurity category by 86%. It more than doubled the chance to be in the category of moderate/severe insecurity. Predictors of health status, lifestyle, and work are not associated with any multinomial outcome categories. However, working in the South, Southeast, or Midwest regions of Brazilian decreased the chances of belonging to one of the FI categories, with significance only for the mild category. Variables that show an association for this population are per capita household income for the different levels of FI and the CR region for mild FI. A high prevalence of FI in this population points to the need for more studies with low-income workers to prevent FI and its health consequences.

**Keywords:** Brazil; community restaurants; food handlers; food insecurity; low-income

#### **1. Introduction**

Food insecurity is a global public health problem, affecting more than 2 billion people worldwide [1]. Food security is a basic need covering access to safe, sufficient, and proper nutritional food [2,3]. Among adults, food insecurity is related to the high prevalence of chronic health problems that compromise this population's quality of life and longevity. It is strongly linked to economic vulnerability [1,4]. In Brazil, the estimation is that more than 14 million households have some food insecurity level [5]; additionally, people eat more outside their homes because they have less time to prepare meals. In Brazil, to the low-income population, eating out tends to represent an intake of cheaper and fast snacks, usually representing inadequate nutrient intake that influences their health outcomes increasing the risk of chronic diseases and nutrient deficits [2,6–8].

In Brazil, people with no or limited access to adequate, safe, and nutritional food are characterized in food insecurity situations (FIS) [5]. FIS can generate fear of the inability to obtain food (influencing quantity or quality of food choice) or generate hunger in the most severe cases due to the scarcity of food at home [5]. The FIS determinants can present social,

**Citation:** Fideles, I.C.; Akutsu, R.d.C.C.d.A.; Barroso, R.d.R.F.; Costa-Souza, J.; Zandonadi, R.P.; Raposo, A.; Botelho, R.B.A. Food Insecurity among Low-Income Food Handlers: A Nationwide Study in Brazilian Community Restaurants. *Int. J. Environ. Res. Public Health* **2021**, *18*, 1160. https://doi.org/10.3390/ ijerph18031160

Academic Editors: Paul B. Tchounwou and Lorrene D. Ritchie Received: 24 November 2020 Accepted: 26 January 2021 Published: 28 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

economic, and political nature, affecting populations' health. The living and working conditions of individuals and population groups are related to their health situation [9]. Among the factors associated with FIS, mainly in the moderate and severe categories of FI, income, and weight excess have been highlighted [10]. Women in developed countries have higher chances of being overweight and obese when in FIS [11].

Based on this information, the Brazilian Government developed the community restaurant (CR) program as a strategy to face the FIS. CRs were created to offer safe, cheap, and healthy meals to the low-income population [12,13]. The CR Program is one of the programs integrated into the "Fome Zero" network of actions, a social inclusion policy established in 2003 [13]. From the beginning, the Government planned to increase the number of CR distributed throughout the Brazilian territory located in regions with more significant numbers of low-income people, reaching 135 CRs in 2020 [13,14]. The production and distribution of CR involve professionals with different education and income levels [15], highlighting the food handlers that directly produce meals [13]. Among food handlers, it is common to have weight excess [15–18], low education levels, and low-income [15,19], making them more susceptible to FIS [10] even though they work in places that produce a great amount of food. Godoy et al. [20] studied food insecurity among CR customers. Despite a high percentage of FI (40.6% for males and 43.8% for females), there were no significant associations between FI and Body Mass Index or body fat percentage. There were significant associations between FI and household income and educational levels. There was also an association with Brazilian's regions among females, FI being worse in the North and Northeast regions.

Few studies evaluated FIS among workers in Brazil [21–23] and in the world [24–27]. However, only three of these studies were performed with professionals working in food service. A study showed a high prevalence of FIS perception among these professionals in Canada compared to other professionals [24]. The two studies conducted in Brazil only evaluated food handlers in a specific region without a nationwide Brazilian representative sample of food handlers [21,22]. Studies with the population of food handlers demonstrate their susceptibility to food insecurity [21,22], to excess weight [22,28], to management in foodservice with a greater focus on the costs involved with the production of meals than with the health of workers [29]. Even though the studies found in Brazil were both developed in the format of case studies [30–32], geographically delimited to institutional spaces, such as universities, Brazilian states or municipalities [21–23,28,31–33], food insecurity is multifactorial and, in CR, we hypothesize that socioeconomic and demographic factors of food handlers are associated with food insecurity. Since the program tends to increase the number of CR in Brazilian territory, more CRs open over the country, possibly resulting in more opportunities to decrease FI among customers and more work opportunities for food handlers. However, these food handlers need to be closely investigated because of their well-being, quality of life, and food security to guarantee motivation to work and meals safely produced. Therefore, the objective was to evaluate nationwide the food insecurity among Brazilian CR food handlers and its associated factors. As it has an exploratory character, the study sought to answer the following question: What are the factors associated with food insecurity among food handlers in Brazilian CR? We expect to provide data allowing to promote interventions to reduce food insecurity among this vulnerable group that works inside foodservices but does not overcome FI inside their families.

#### **2. Materials and Methods**

#### *2.1. Design, Settings, and Participants*

This cross-sectional study was performed in Brazilian CR (focused on the low-income Brazilian population offering meals from Monday to Friday). It was conducted according to the Declaration of Helsinki guidelines and approved by the University of Brasília Research Ethics Committee (#037210). Written informed consent was obtained from all participants.

The basis of the sample calculation was the official list of all existing CR throughout the Brazilian territory at the moment of data collection [13]. With CR nationwide selection,

food handlers in different Brazilian regions and similar work conditions were evaluated. Researchers were allowed to enter all of the CR with nationwide permission. The restaurant inclusion criteria were to be part of the Brazilian CR program, sign the Institutional Acknowledgement Agreement by the CR responsible, and offer daily more than 500 meals. With more than 500 meals, CR allowed researchers to evaluate many workers, helping to achieve a representative sample. Restaurants that provide less than 500 daily meals are considered small and present fewer food handlers than medium or large restaurants that offer more than 500 daily meals [34]. There were 65 existing CR in Brazil, and all of them met the inclusion criteria.

From the 65 CR, a sampling plan was calculated considering a level of significance (α) of 5% [35] using the "survey select" of the SAS 9.1.3 program (SAS Institute Inc., Cary, NC, USA). Therefore, a minimum of 31 CR existing in each of the five Brazilian regions (North, Northeast, Midwest, South, and Southeast) was randomly selected to be part of the study. The final sample of 36 restaurants was used, respecting the stratification criteria by the Brazilian region. The distribution of randomly drawn CRs was proportional to the number of CR in each region. The researchers visited 36 CR and included them in the sample (4 in North, 10 in Northeast, 1 in Midwest, 15 in Southeast, and 6 in South). All of the working food handlers in the 36 drawn CR were invited to participate in this study (*n* = 1062). Data collection was carried out over four consecutive days to cover the various work schedules to guarantee access to the largest possible number of handlers in each CR. Therefore, from the total of 1062 handlers working on the selected CRs, 970 met the inclusion criteria (e.g., not being pregnant, and workers on vacation or not working due to medical issues). From them, 471 (48.6%) agreed to participate and completed the study.

Some handlers refused to participate in the research because they did not want to stop working or worried about employability. Even though researchers explained that participation would not influence their work, some decided not to participate in the study. A total of 383 individuals was necessary to be a representative sample stratified in the five Brazilian regions according to the number of food handlers in each region. Participants were not compensated for the participation. They were just informed about the importance of the study for their category. The study used a 95% confidence interval and an error of 4%, respecting the minimum handlers' sample per Brazilian region [35]. Therefore, it was necessary to achieve a minimum of 11% of the handlers from the North, 28% from the Northeast, 3% from the Midwest, 41% from the Southeast, and 17% from the South.

#### *2.2. Data Collection*

Trained researchers performed data collection using standardized instruments to identify socio-demographic characteristics and the Brazilian Food Insecurity Scale (EBIA/BFIS) [36]. The socio-demographic variables were gender, age, educational level, per capita household income, marital status, smoking status, participation in a governmental program, and how many years or months the individual has worked in the CR. The presence of diagnosed non-communicable diseases (NCDs) (by a physician) was self-reported. The information was recorded in a specific form showing the presence or absence of one or more than one of the following NCDs: Systemic arterial hypertension (SAH), type 2 diabetes mellitus (DM), and others (cancer, dyslipidemia, cardiovascular diseases, respiratory diseases, depression). The self-reported NCD data was used because population studies widely use this method for its convenience and economy [37,38].

The individuals who agreed to participate signed the consent form after receiving information about the research. Before participants' lunch in a reserved room, weight and height were measured with a Plenna® (São Paulo, Brazil) weighing scale (150 kg) and a stadiometer (220 cm). Participants had to take off their shoes and coats. After that, the body mass index (BMI) was calculated [39]. The anthropometric status based on body mass index classification [39] was dichotomized on weight excess (0. No; 1. Yes). The cut-off point used to indicate excess weight was a BMI value <sup>≥</sup> 25 kg/m<sup>2</sup> , which covers both the category of overweight and obesity.

#### *2.3. Dependent Variable*

The study's outcome, household food security situation (HFSS), was obtained through EBIA/BFIS adapted and validated for the Brazilian population [40,41]. The EBIA/BFIS seeks to assess the perception and experience of household residents' hunger in the three months before the instrument application. The positive answers to the 14 questions of the instrument categorize the level of food security/insecurity (considering the age of the residents) in food security (0 points); mild food insecurity (1–5 points in the presence of residents <18 years old, or 1–3 points in the absence of residents <18 years old); moderate food insecurity (6–9 points in the presence of residents <18 years old, or 4–5 in the absence of residents <18 years old) and, severe food insecurity (10–14 points in the presence of residents <18 years old, or 6–8 points in the absence of residents <18 years). Based on the methodology adopted by Panigassi et al. [42], for this study, the HFSS categories were grouped into food security, mild food insecurity, and moderate/severe food insecurity.

#### *2.4. Statistical Analysis*

The data were analyzed with the STATA 15.0® (StataCorp LP, College Station, TX, USA), using frequencies to describe the categorical variables and Pearson's chi-square test (χ 2 ) to identify associations (*p*-value < 0.05).

The outcome with three categories of HFSS, multinomial (polytomous) logistic regression models were used. They are applicable for outcomes with three or more levels [43]. For this study, the category "food security" was chosen as a reference. In this type of analysis, the logistic model will have two logit functions: the ratio between Y = 1 and Y = 0 and the ratio between Y = 2 and Y = 0, with Y = 0 as the referent. The numbers 0, 1, and 2 represent the food security status classification for statistical analysis—(0) Food Security—the household has regular and permanent access to food in sufficient quantity and quality without compromising access to other needs; (1) Mild food insecurity—at this level there is uncertainty regarding access to food in the future, with a change in the quality of food, but without compromising the amount of food; (2) includes the levels of moderate food insecurity and severe, in which there is already a quantitative reduction in access to food, causing changes in the dietary pattern of residents at home. These categories were grouped to enable comparison with population data given that national food insecurity assessment studies [44] and international [45] present their analyses considering the grouped on these two levels.

Bivariate multinomial regression models, with HFSS as the dependent variable, were applied to all predictor variables. This stage helped to understand the initial associations of the determinant factors for the HFSS of this population and the magnitude of each predictor's effect by calculating the gross Odds Ratio (OR) and their respective 95% Confidence Intervals (CI). The predictor variables were included in the multivariate model based on two assumptions. The first assumption was the theoretical basis underlying the HFSS. The second one was the statistical decision based on the bivariate multinomial model result with a value of *p* < 0.20 [16]. Thus, the multivariate multinomial model was composed of the predictors that met the assumptions and were adjusted together, without considering hierarchical order or determination level.

The backward stepwise procedure was used for the selection in the final model. The Likelihood Ratio Test (LRT) was used to test hypotheses about the significance of the predictor variables, that is, to evaluate the effect at all levels of the outcome simultaneously, which affects the number of parameters tested and the degrees of freedom associated with the test [46]. Thus, the comparison of the observed and the expected values using the likelihood function was based on the expression:

$$D = -2lnL\_{reduced} - \left(-2lnL\_{full}\right) \sim X^2$$

The full model corresponds to the complete multivariate model and the reduced model to the model without a corresponding predictor, following a chi-square distribution with degrees of freedom equal to the number of set parameters (defined as zero under the null hypothesis—H0). Statistically significant results (*p* < 0.05) reject H0 and indicate that the predictor is significant for the model, being maintained in the final model.

#### **3. Results**

The socioeconomic and demographic characteristics, health status, HFSS, and aspects related to food handlers' work in CR are in Table 1. The sample was composed mainly of women (62.2%), aged ≤40 years old (67.7%), with a partner (52.0%), and with up to eight years of education (54.1%). In most households, there were up to two people employed (78.8%), less than six rooms (64.4%), per capita household income below half a minimum wage (40%), and 19.3% of the workers in government programs. Among the health risk behaviors or conditions, overweight and alcohol use was more prevalent (60.8% and 46%, respectively) when compared to the smoking habit (17.2%) and the presence of NCDs (19.1%). More than half of these professionals were working in the CR for ≥12 months (57.3%) (Table 1). For most of the factors in Table 1, the median value was used to split the sample and perform the analysis.

Predictors of socioeconomic status (education, per capita household income, and participation in governmental programs) and geographic location of the CR were associated with the HFSS categories by the chi-square test (*p* < 0.05) (Table 1).

There was no significant association between mild FI and the demographic variables in the bivariate multinomial analysis. In contrast, the predictors of socioeconomic conditions were associated with both the mild and the moderate/severe FI category. Workers with less education were twice as likely to belong to the category with the highest severity of food insecurity (OR: 2.17; 95% CI: 1.25–3.77). Having lower per capita household income increased the chances of belonging to the category of mild insecurity by 86% (OR: 1.86; 95% CI: 1.20–2.88) and more than doubled the chance to be in the category of moderate/severe Insecurity (OR: 3.8; 95% CI: 2.18–6.60) in addition to participating in governmental programs (OR: 3.17; 95% CI: 1.75–5.74). Living in households with fewer rooms was also associated with the most severe food insecurity category (OR: 1.87; 95% CI: 1.05–3.34) (Table 2).

Predictors of health status, lifestyle, and work were not associated with any of the multinomial outcomes' categories at a 5% significance level. However, working in the South, Southeast, or Midwest regions of Brazilian decreased the chances of belonging to one of the Food Insecurity categories between 21% to 42%, with significance only for the mild category (OR: 0.58; 95% CI: 0.38–0.88) (Table 2).

The bivariate multinomial regression models confirmed the associations of education, per capita income, participation in a governmental program, and the Brazilian region with the categories of HFSS, with an increase in the number of rooms among the predictors. Besides, variables that showed association with any of the categories of HFSS (*p*-value < 0.20) were also included in the multivariate multinomial model (gender, number of people working at home, and the presence of NCDs) (Table 2).

*Int. J. Environ. Res. Public Health* **2021**, *18*, 1160



*Int. J. Environ. Res. Public Health* **2021**, *18*, 1160



**Table 2.** Bivariate multinomial logit model with odds ratio (OR) estimation between household food security situation (reference group: food security) and demographic, socioeconomic, health, lifestyle predictors, and those related to work in a Community Restaurant in all the Brazilian regions.


OR: Odds ratio; SE: standard error. \* *Chi*-Square test statistic *p*-value. 1 Minimum wage (MW): 175.90 USD per month.

Table 3 shows the results of the multivariate model containing eight predictors selected in the previous steps. They were associated with at least one of the categories of HFSS with a value of *p* < 0.20. In the complete crude model, per capita income (Mild Insecurity— OR: 1.64; 95% CI: 1.01–2.66/Moderate/Severe Insecurity—OR: 2.7; 95% CI: 1.48–4.94); participation in a governmental program (Moderate/Severe Insecurity—OR: 2.02; 95% CI: 1.04–3.92) and CR region (Mild insecurity—OR: 0.5; 95% CI: 0.31–0.80) maintained an association with HFSS, increasing the chances of a CR worker belonging to one of the risk categories for food insecurity when in unfavorable socioeconomic situations, and compared to those with better socioeconomic conditions.

A stepwise backward method was applied to the multinomial model to adjust the final HFSS model. Calculating the odds ratios and standard errors and the model's goodness adjustment indexes guided the model's choice that best-suited data. Reduced models with smaller Akaike information criterion (AIC) than the previous model indicated the removal of the tested variable, as well as the *p*-value (>0.05) of the LRT (Table 4). The variables per capita household income, number of rooms, and number of people working (at the same home), in addition to the CR Brazilian region, remained in the model (Table 5).

Despite being maintained in the final model, considering the LRT results, the number of rooms and the number of people working did not maintain the association with food insecurity in the adjusted model. When tested, they were not confirmed (range of estimates < 10%). The multivariate model results confirmed the statistical significance of per capita household income as a predictor of HFSS. Food handlers with per capita household income equal to or less than half the minimum wage were 1.77 times more likely to be in a situation of mild food insecurity (OR: 1.77; CI95: 1.12–2.80) and around three times more likely to be in a situation of moderate/severe food insecurity (OR: 3.52; CI95: 2.00–6.19) when compared to individuals belonging to the food security category. In the adjusted model, the association with the CR Brazilian region was maintained in mild food insecurity. Therefore, the food handlers working in CR located in the Midwest/Southeast/South regions were 47% less likely to experience mild food insecurity (OR: 0.53; CI95: 0.34–0.83) compared to individuals in food security (Table 5).


**Table 3.** Multivariate multinomial logit model for predictors of Household Food Security Situation (reference group: food security) in food handlers of community restaurants.

OR: Odds ratio; SE: standard error.



**Table 4.** Comparison of the complete and incomplete model performed using the Likelihood Ratio Test (LRT)—stepwise backward method.

LRT: Likelihood Ratio Test; ll = Log likelihood; gl = degrees of freedom; AIC = Akaike information criterion; BIC = Bayesian information criterion.

**Table 5.** An adjusted multivariate model of the Household Food Security Situation (HFSS) (reference group: food security) determinants among food handlers in Brazilian community restaurants.


Model adjusted for income, number of rooms, and workers in the household and region of the CR. 95% CI: 95% confidence interval; OR: Odds Ratio; SE: standard error. \* significance level *p* < 0.05.

#### **4. Discussion**

Food insecurity is a public health problem, affecting individuals in all parts of the world. The latest survey conducted in 2018 by the Food and Agriculture Organization of the United Nations (FAO) showed that about two billion people suffered from moderate or severe food insecurity, representing 26.4% of the world population [1]. In Latin America, including Brazil, 188 million people had the perception of moderate to severe food insecurity, representing 30.9% of the population [1]. In Brazil, the last populational study that assessed food insecurity in households was the National Household Sample Survey (PNAD/NHSS), carried out in 2013 and showed a prevalence of 22.6% of food insecurity in the population [36].

In Brazil, the foodservice market employs about 250 thousand food handlers [47], who work in several foodservice segments, such as this study's locus (the CR). The findings can be extended to others in the foodservice segment, such as the food handlers working in hospitals, industrial or commercial food services. The activities of workers in the Brazilian foodservice sector revolve around a common goal, the production of meals for groups, whether healthy or sick [47], and, in the case of CR handlers, they provide meals for groups in social vulnerability [13]. The results of this study can be extended to the population of food handlers, given the peculiarities of this sector that include: (a) the composition of teams with professionals from different levels of educational background and with different and complementary activities [31]; (b) by workers from the different foodservice segments with a high prevalence of excess weight [30,32,48–50]; (c) by workers without professional training to carry out the activities, the in-service learning process takes place, which also leads to high turnover in the sector [31]; (d) inadequate work conditions, characterized by the requirement of long hours, rhythm and intense efforts with work overload, tasks such as repetitive movements, leading to absenteeism and high turnover of workers and the prevalence of work-related diseases [31]. This scenario reinforces the importance of the study and its scope.

Our results showed a total prevalence of 44.4% of food insecurity, much higher than the general population [51]. Moderate/severe food insecurity was almost twice the prevalence of the Brazilian population. In this study, socio-demographic characteristics showed that most food handlers were females in the age group ≤40 years old, similar to previous studies carried out with food handlers [15,16,18,22].

Two studies on food handlers from CR, the first in the city of Rio de Janeiro (*n* = 273 individuals from 7 CR) [21] and the second in the city of Belo Horizonte (*n* = 180 individuals from 4 CR) [22] showed a prevalence of total food insecurity of 53.7% and 24%, respectively. This prevalence, as in our study, is higher than food insecurity for the Brazilian population. It is worth mentioning that in these studies [21,22], most food handlers perceived themselves in mild food insecurity, according to the EBIA/BFIS classification [36]. There is a concern with access to food in the future, and it may impact the quality of the diet, choosing to buy cheaper foods [52].

In our study, an association was found between income and food insecurity, maintained in the regression model. Of the study participants, 39.9% had a per capita income at home equal to or less than half the minimum wage, classified as low-income families, according to the Brazilian Government [53]. Moreover, income showed a more significant impact on households with moderate to severe insecurity. Low-income families were about three times more likely to have food insecurity than workers with income above half the minimum wage. Other variables related to household income were associated with food insecurity as the number of rooms in the house and the number of individuals working, also maintained in the final regression model. Another variable related to household income that showed an association with food insecurity was the participation in a governmental program. Almost 40% of participants presented the requirements to receive governmental benefits [53], but only 19.3% of the participants received any benefit. Godoy et al. [20] also showed differences in food insecurity prevalence and income and education for CR customers. Among males, food insecurity prevalence was 64.2% when per capita income

was below <sup>1</sup> <sup>4</sup> minimum wage and 54.8% from <sup>1</sup> 4 to half of the minimum wage. Females showed similar results, with 61.5% in food insecurity when per capita income was lower than <sup>1</sup> <sup>4</sup> minimum wage, and 54.8% when income was between <sup>1</sup> 4 to half of the minimum wage. CR customers with less than eight years of education also presented higher food insecurity for males and females (51% and 56.3%, respectively) [20].

In a previous study with the Brazilian population, 54.8% of households with moderate to severe food insecurity received per capita income up to half a minimum wage [51]. In the study carried out with food handlers in the CR in Rio de Janeiro, 59.4% of workers reported having another job to increase their income [21]. Per capita household income showed an inverse association with food insecurity in a study developed with food handlers in CR in Belo Horizonte/Brazil [22].

In a study that assessed food insecurity among workers in Canada, multivariate regression revealed that increased income independently decreases the chances of food insecurity in Canadian families, as well as having more workers in the family or at home. Income was considered a significant factor for the chances of food insecurity at home [24]. Some studies show an association between FI lower per capita income at home [53–55].

Schooling was associated with food insecurity at home in the bivariate analysis (*p* = 0.02), not remaining in the final regression model. In our study, 54.1% of the participants had an education level equal to or less than eight years, in agreement with other studies conducted with this professional category, identifying the highest percentage of individuals with educational level up to eight years of formal study (equivalent to elementary school) [19,21,34]. In general, education is associated with food insecurity among the general Brazilian population [51], in which a higher level of education relates to the lower prevalence of moderate or severe insecurity. Since food handlers in Brazil receive low salaries and present fewer years of formal education, they are more related to food insecurity even though working with food production.

In the study by Falcão et al. [21], an association was found between education and food insecurity in food handlers' homes. In the studied sample, the risk of those who had up to nine years of formal study was almost three times more likely to find themselves in a food insecurity situation than those with higher education.

Another variable that showed association with food insecurity was the region of the CR. Food handlers from the Midwest-South regions perceived themselves in household food security. The National Household Sample Survey (PNAD/NHSS) [47] identified a higher prevalence of food insecurity in the North and Northeast regions (30%) than the Midwest/Southeast/South regions (less than 20%). Bezerra et al. [54] also demonstrated higher FI in the North (40.3%) and the Northeast (46.1%) with positive and moderate correlation between FI and the percentage of the extremely poor population. Gubert et al. [55] demonstrated a reduction in FI prevalence in all Brazilian states from 2004 and 2013, higher in the Northeast region. Santos et al. [56] evaluated the tendency and associated factors of food insecurity in Brazil, analyzing the NHSS from 2004, 2009, and 2013. These researchers verified that even though there was a reduction in food insecurity levels with time, there was an increase in the association force between food insecurity and regions North and Northeast of Brazil. Residents from these regions had three to four more chances of having moderate to severe food insecurity. With CR customers, only for females, there was an association between region and food insecurity, with higher prevalence in the north and northeast regions [20].

This study did not show an association between food security and excess weight, contrary to other studies [57–60] that showed an association between food insecurity and excess weight. Shamah-Levy et al. [60] associated food insecurity with undernutrition and obesity, discussing that both situations come from poor eating habits. In children, Kac et al. [58] discussed a less consistent association with overweight, but Franklin et al. [57] studying women presented a higher prevalence of excess weight with food insecurity. Godoy et al. [61] did not find a significant association between food insecurity and BMI or body fat percentage for CR's male or female customers. Even though obesity is included in

the NCDs, excess weight was not included as NCD for this study's associations. There was an association between food insecurity and NCDs; however, it was not maintained in the model.

One of this study's contributions is the confirmation of the dichotomy between FI and the performance in the food field. In the case of this study, places designed to promote access to food (CR) to prevent the occurrence of food insecurity also constitute a workspace for people in FI. Besides, the geographic location of the CRs act as protectors, confirming that regional inequalities in Brazil affect the same population group differently and confirming their claim that income is associated with FI. It is important to emphasize that the mean income from the North and Northeast regions is almost half of the income from workers from the South and Southeast regions in Brazil. This inequality is also observed in this food handler group and needs to be discussed inside a governmental program.

#### **5. Conclusions**

In Brazil, few studies assess workers' food insecurity. Food handlers play an essential role in promoting food security as responsible for producing and supplying meals in hygienic and nutritional conditions suitable for the population using this state equipment. These workers have characteristics, confirmed in the study, as low education and low income, which places them both as actors and target population of public policies aimed at food insecurity. This study is the first to assess food handlers nationwide, and a high prevalence of food insecurity was observed, higher than that of the general population. For moderate/severe food insecurity, it was almost twice the Brazilian population. Based on the study design, it is not possible to establish a causal link between the variables, but the study shows paths for further studies. The variables that showed an association for this population were the household income per capita for the different levels of food insecurity and the CR region for mild food insecurity. The high prevalence of food insecurity in this population and the return of the increase in the prevalence of food insecurity in the Brazilian and worldwide population point to the need for further studies on the theme with categories of workers, mainly low-income, which allow identifying factors related to the work environment and developing policies and interventions aimed at preventing food insecurity and its consequences for health.

**Author Contributions:** Conceptualization, I.C.F.; R.d.C.C.d.A.A.; R.d.R.F.B.; J.C.-S.; methodology, I.C.F.; R.d.C.C.d.A.A.; R.d.R.F.B.; J.C.-S.; R.P.Z.; R.B.A.B.; validation, I.C.F.; R.d.C.C.d.A.A.; R.d.R.F.B.; J.C.-S.; formal analysis, I.C.F.; R.d.C.C.d.A.A.; R.d.R.F.B.; J.C.-S.; R.P.Z.; R.B.A.B.; investigation, I.C.F.; R.d.C.C.d.A.A.; R.P.Z.; R.B.A.B.; data curation, I.C.F.; R.d.C.C.d.A.A.; writing—original draft preparation, I.C.F.; R.d.C.C.d.A.A.; R.P.Z.; R.B.A.B.; writing—review and editing, I.C.F.; R.d.C.C.d.A.A.; R.P.Z.; A.R.; R.B.A.B.; visualization, A.R.; supervision, R.d.C.C.d.A.A.; project administration, R.d.C.C.d.A.A.; R.B.A.B.; funding acquisition, A.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of University of Brasília (protocol code #037210).

**Informed Consent Statement:** Written informed consent was obtained from all participants.

**Data Availability Statement:** Data sharing is not applicable to this article.

**Acknowledgments:** The authors acknowledge the support of "Ministério do Desenvolvimento Social".

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

