For the objects containing patterns of two dots and a line segment, the chance for humans to see faces is high because the two dots and line segments appear to the pattern of a human face with eyes and a mouth, which is the most common pareidolia [14]. However, neuropathological and psychiatric conditions can have a negative impact on the ability to identify facial expressions of sensation [14]. Human eyes can distinguish brightness differences in a range of 1–100, which means human eyes' dynamic contract is high. The human eye is a brilliant detecting system that will undergo both physical and chemical processes from the view and send it to the brain. Though the outcome of the brain–eye system is perception, the brain does not always expound the information of the object seen in the real perception [11]. The most suspect human perceptual apparatus is sight because human eyes do not see a centralized field, while the saccadic movement of the eyes stores the image from detail to detail and is perhaps used later for data recall by forming a whole image [19]. The rhizome structure in the brain has many chemical reaction saccades for the first time of seeing, and it will leave many impressions unprocessed, disordered, and unknown. Since the human brain is a complex structure, the percipient information received in real-time will be eliminated, expanded, rearranged, and codified to form a common and logical layout of the external world. According to Wertheimer and Riezler [20], the Gestalt theory states that visual perception is the effect of the relationship among the objects observed instead of the simple add-up of the elements seen. As a result, the misperception of the shape and color of an object may be caused by the complicated and fast processing work of the brain and lead to illusion and pareidolia.

 Fusiform gyri shows higher activation when a facial expression is shown on the subject [22]. The face-sensitive neurons in nonhuman primates' inferior temporal cortex show the same selectivity for face-like object configurations [23]. Even a newborn baby can undergo face recognition and perform across view change [24]. According to Hoback [3], a 3-month-old baby already starts learning to identify the mother's face, while when the baby grows to four to nine months, the baby can distinguish several facial expressions such as happiness, fear, anger, surprise, or sadness. Moreover, detecting faces from the environment was the survival intuition of humans to ensure vigilance towards danger in the surrounding area [25].

#### *1.2. Pareidolia in Religion*

A German neurologist and psychiatrist, Klaus Conrad (1905–1961), promoted that pareidolia has a bearing on an even wider range of illusory phenomena, including the discernment of religious themes such as the faces of the Virgin Mary and Jesus [8]. Paranormal believers frequently have perceptual illusions in ambiguous visual stimuli [1]. Even though both paranormal believers and non-believers have the same ability to detect face perception, non-paranormal believers have less liberal response bias than paranormal believers [26], which is probably because paranormal believers perceive ambiguous stimuli as face-like patterns more easily. There is an interrelationship between illusory agency detection and paranormal belief in the studies regarding schizotypy and schizophrenia [13].

Both schizotypy and schizophrenia can cause pareidolia. The schizotypal personality sometimes comes with delusional beliefs, abnormal perceptual experiences, and magical thoughts [27]. According to Galdos et al. [28], a person with a schizotypal personality can sense meaningful patterns in random meaningless noise, while a person with schizophrenia is more able to link unrelated events. By using the moving triangle task, the person with schizotypal can gain meaning from geometrical images that move randomly [13]. Furthermore, some natural landscape elements are sacralized in the three mythological characters: Pan-Gu, Fu-Xi, and Shen-Nong. All these three mythological characters belong to the Chinese culture's formative period. The three humanized mountains were sacralized to identify these three outstanding Pareidolian Images in their origin. In the case of the commercial aircraft that destroyed New York City's World Trade Center on 9/11/01, many protestants saw demons and devil faces in the smoke, and they believed that the crisis destroyed the current seat of power of Satan living in the American financial institutions [29]. Pareidolia extends the illusionary connections between landscape and the mythological characters.

#### *1.3. Present Study*

Pareidolia exists in the built environment, such as landscapes, buildings, and furniture. Many pareidolian faces occur as ornaments or adornments in the architectural history of different cultures [14]. Some house envelopes consist of a leading pattern that can be justified as a human face. The facial expressions perceived in houses can trigger an emotional reaction; thus, the techniques of using pareidolia to produce positive emotions are applied in house design [3]. Experimental studies on the perception of illusory agency discovered that the tendency to sense human faces among random noise is high, which is more than 40%, especially in low information content [17]. There are interrelationships between pareidolian faces and aesthetic value [30]. Though there are many kinds of stimuli able to trigger the emotion of observers, a human face is the most effective one [5]. Examples of pareidolia in the landscape are faces in trees, faces on mountains, faces on stones, the man on the moon, and faces in clouds [18]. Pareidolia sometimes appears on house facades, windows of buildings, doors of buildings, electrical outlets, cupboards, and chairs [14]. Previous research validated the relationship between demographics and delusional disorder using questionnaire surveys and pilot experiments [3,14]. Mental anomalies may be the source causes of many tragedies and incidents [31]. In addition to job pressures, lifestyles and other aspects may affect the mental state to varying degrees. In fact, everyone will be in different mental states at different times of normal life. In general, people are in blurred thinking status when they wake up, while high brain activity during the day will leave people exhausted at night, and the brain activity time is negatively correlated with the human spirit [32]. Negative spirit status increases the possibility of pareidolia [33]. Therefore, this research compares built environmental pareidolian phenomena at different spirit state time points (6 a.m., 12 p.m., and 2 a.m.) with the purpose of determining people's sensitivity and reactions towards pareidolia in a built environment.

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

We applied a mix of qualitative and quantitative approaches in this study. We combined a pilot experiment and a questionnaire survey to triangulate the data collected in this study. The pilot experiment is an explanatory and descriptive method to categorize the phenomena of pareidolia in the built environment, which was then used to structure the questionnaire survey and provide readers with a better understanding of pareidolia. The qualitative approach relies on evidence rather than frequency to illuminate issues and uncover possible explanations.

#### *2.1. Participants*

#### 2.1.1. Pilot Experiment

According to Liu [18], 20 participants (10 females and 10 males, aged from 18 to 60) were randomly recruited from society for the pilot experiment and were offered financial remuneration for participation. The participants mainly consisted of five age groups: 11–20, 21–30, 31–40, 41–50, and 51–60. Four participants per age group. All participants had normal or corrected vision.

#### 2.1.2. Questionnaire Survey

The anonymous questionnaire survey was distributed to 500 people ranging from 11 to 60 years old through the post. In this study, the power was set to 0.8, and the threshold for significance (α) was set to 0.05. According to the power analysis, if we wanted an 80% probability of correctly rejecting the null hypothesis, we needed a sample size of at least 217. A total of 228 valid forms were returned and assessed, representing a 45.6% response rate. The demographic details of the respondents are shown in Table 1. The percentage of females at 59.6% was slightly higher than that of males at 40.4%. The questionnaire survey targeted five age groups: 11–20 years old (6.6%), 21–30 years old (75.4%), 31–40 years old (6.6%), 51–60 years old (6.1%), and 41–50 years old (5.3%). Approximately 72.8% of respondents had a bachelor's degree, followed by 10.1% with secondary education, 9.2% with a diploma/advanced diploma, 5.3% with primary education, 2.2% with postgraduate degree/professional level, and 0.4% without formal education. In total, 49.6% of respondents considered their current job stressful, and 50.4% did not. There were 62.7% single respondents and 37.3% married. Moreover, 55.7% of respondents normally went to bed after 12 a.m., and 34.6% of respondents did not have any habits such as smoking, alcohol abuse, staying up late, or chronic diseases.

**Table 1.** Demographic details of respondents.

#### *2.2. Stimuli*

In total, 26 pareidolian images were used for this pilot experiment. The images were obtained from building photos taken on-site and were presented on paper with a picture size of 7.5 × 7.5 cm (H × W). A pareidolian face consists of a facial outline, left and right eyes, and a mouth [30]; thus, the images used in the experiment had similar characteristics.

#### *2.3. Experimental Setup and Procedure*

#### 2.3.1. Pilot Experiment

In the pilot experiment, 20 participants were required to observe the 26 images provided by the researcher during the morning, midday, and midnight and share their feedback. Each image was observed as an individual experiment by 20 participants. Their feedback was then analyzed using content analysis and structured through tabulation.

#### 2.3.2. Questionnaire Survey

Sensitivity in this study refers to a change in pareidolia's likelihood when one factor changes. The questionnaire survey determined sensitivity and reactions towards pareidolia in a built environment. There were 69 questions in the questionnaire form. Section A focused on the respondents' background, such as age, gender, occupation status, education level, lifestyle, and stress level ("Gender": male, female; "Age group": 11–20, 21–30, 31–40, 41–50, 51–60; "Educational Level": no formal education, primary education, secondary education, certificate/diploma/advanced diploma level, degree, post-graduate degree/professional level; "Current job stress level": stressful, not stressful; "Current status": single, married; "Lifestyle": smoking, alcohol abuse, sleep after 12 a.m., long-term medicine user). Section B focused on their sensitivity and reactions towards pareidolia in different spaces in the building. The questions offered a 5-point Likert-type scale for respondents to rank.  Investigating people's emotional feedback determines what type of architectural design is more likely to be acceptable rather than radical when people generate pareidolia. The specific content of the questionnaire is presented in the Appendix A.

Two software were applied in the analysis: Statistical Package for the Social Sciences (SPSS) Version 22 and Smart PLS 2.0. Structural Equation Modeling (SEM) was used to test the theoretically supported model.

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

#### *3.1. Findings in Pilot Experiments*

Data collected from the pilot experiments were used to compare the participants' sensitivity and reactions toward pareidolian images at different time points (6 a.m., 12 p.m., and 2 a.m.). There were 26 experiments conducted among the 20 participants, and the results are tabulated in Table 2. In general, the participants had very similar sensitivity and reactions towards these pareidolian images. Most participants could identify human faces at 6:00 a.m. in the morning and 2:00 a.m. at night, but not at 12:00 p.m. midday. Most participants felt scared when looking at the pareidolian images at 2 a.m.

**Table 2.** Tabulation of experiments results.

#### **Table 2.** *Cont.*

#### **Table 2.** *Cont.*

According to the feedback from pilot experiments, face identification in the built environments was easier at midnight and in the morning. Especially at midnight, many participants described the faces as scared, whereas they only considered those images as normal pareidolian faces at other time points (6:00 a.m. and 12:00 p.m.). Therefore, the results suggested that pareidolia was more likely to occur when people were exhausted. In other words, when people are fatigued, it will be easier to identify pareidolian faces. Face identification was not significant at 12 p.m., probably because of a higher level of relaxation in brain activity; thus, it did not reach exhaustion at that time. Results from partial images' feedback showed people reacted differently between morning and midnight when identifying a face in Pareidolian Image 2, 10, 11, and 13. The people's reactions to images oscillated from "Funny" in the morning to "Scary" at midnight. The majority of pareidolian images reactions were "Scary", "Neutral", and "Creepy". As a result, the pilot experiments revealed that people were more likely to identify pareidolian faces at midnight and in the morning than at midday, and people reacted more negatively to pareidolian images at night.

#### *3.2. Results of Questionnaire Survey*

The Cronbach's Alpha of the questionnaire survey was 0.899, showing a high internal consistency and high reliability. In the study, factor analysis was first used to screen the images and eliminate the image results that could not be well classified by factors. The factor analysis' eigenvalue was set at 0.9, and the factor loadings of the original 21 images with an absolute value greater than 0.3 were considered acceptable. There were two criteria for screening images: first, one image had a large loading on two or more factors (factor loading > 0.4), and the difference between the absolute value of the factor loading was less than 0.3; and second, only one image under a factor had factor loading. Only one image could be removed per screening and required reanalysis. The purpose of factor analysis is to integrate highly correlated factors to form new factors. The final KMO value was 0.938 (KMO value > 0.6), and the significance of the Bartlett sphericity test was 0.000 (*p* < 0.05), indicating that the final questionnaire had excellent validity, and the last retained image was suitable for factor analysis. We used factor analysis to remove pareidolian images 1, 4, 5, 6, 7, 8, 9, 11, 13 and 21; the purpose of each image was to examine the effect of demographic information on pareidolia; thus, we could integrate the deleted images into one factor through factor analysis. The final new factor extracted the results from the majority of the images; thus, it could well represent the overall results of the returned questionnaire. We named this factor the pareidolia questionnaire result and used it for multivariate analysis. The image factor loadings in the pareidolia questionnaire result are shown in Table 3.

**Table 3.** The image factor loadings in the pareidolia questionnaire result.

Table 4 show the results of effect analysis between the demographic information in the pareidolia questionnaire results. Among them, age, gender, and lifestyle as a single factor showed a significant impact on pareidolia. Many two–factor between–subject effects also significantly influenced pareidolia identification, including Age\*Gender (*p* = 0.036 < 0.05), Age\*Lifestyle (*p* = 0.000 < 0.05), Gender\*Lifestyle (*p* = 0.000 < 0.05) and Job Stress\*Lifestyle (*p* = 0.000 < 0.05). It was worth noting that the three–factor and four–factor interaction effects did not show a significant influence on pareidolia identification.

**Table 4.** Tests of Between–Subjects Effects.

Figure 1 show the effect of significant two–factor between-subject effects on pareidolia identification. Two–factor estimated marginal means can not only show between–subject effects but also reflect the influence of a single factor on pareidolia. Figure 1a show the differences in pareidolia identification between males and females in different age groups. From the overall trend, females were more likely than males to identify pareidolian faces. Especially in the 21–40 age group, females performed much better than males in pareidolian face identification. However, gender had no significant influence on face identification in the 41–60 age group. Compared with gender, age had a significant impact on pareidolia. The performance of pareidolian face identification improved with increasing age. It seemed that 40 years old was a watershed since the 41–50 age group displayed a significant difference from the 31–40 age group. When the age exceeded 40 years old, the performance of pareidolian face identification showed a steady trend.

**Figure 1.** Two–factor Estimated Marginal Means of pareidolia questionnaire result. (**a**) Gender and age group; (**b**) Gender and lifestyle; (**c**) Job stress and lifestyle; (**d**) Age group and lifestyle.

Figure 1b show the differences in pareidolia identification between males and females with different lifestyles. Lifestyle significantly influenced pareidolia identification. The two– factor between-subject effects of lifestyle with age, gender, and job stress all demonstrated significant influence. Figure 1b,c both showed that long-term medicine users were more likely to identify pareidolia. Topers and smokers behaved similarly, but both influenced pareidolia identification to a large extent. Sleeping late had little influence on pareidolia recognition. However, people without these lifestyles seemed to be less prone to pareidolia. Figure 1b show the differences in pareidolia between males and females with different lifestyles. Male smokers and long-term medicine users were more likely to identify human faces than females. However, females who abused alcohol, slept late and did not have these lifestyles were more likely to identify human faces than males. This also confirmed that females were more likely to identify human faces than males in daily life. Figure 1c show the between–subject effects between job stress and lifestyle. There was little difference in pareidolia identification between long-term medicine users, smokers, and people without these lifestyles, whether they had job stress. Among late sleepers, people with job stress were more likely to identify human faces. Conversely, topers without job stress appeared to easily identify human faces.

Age and lifestyle were two factors that significantly affected pareidolia; the betweensubject effect remained significant for pareidolia identification. From Figure 1d, people in the 41–60 age group had significantly higher pareidolia identification between different lifestyles than those in the 11–30 age group. Other lifestyle differences were revealed to not be significant in the 41–60 age group, except for long-term medicine users. However, in the 11–40 age group, long-term medicine users had a significant impact on pareidolia. Topers

and smokers had a large degree of sensitivity to identifying human faces in images. Other lifestyles did not show a significant influence on pareidolia.

#### *3.3. PLS–SEM Analysis*

This study used the PLS–SEM model to investigate the causal relationships between pareidolian images, pareidolia sensitivity, demographics, and pareidolian reactions. Since PLS modeling focuses on prediction and makes no assumptions regarding data distribution, residuals, or parameters of the observable variables, only generic linear regression requirements must be met. PLS modeling is more reliable than LISREL in analyzing the causal relationship between pareidolia factors. Figure 2 show the initial PLS–SEM model of this study.

**Figure 2.** Initial PLS–SEM model of demographic data and sensitivity to pareidolia in the built environment.

To verify the rationality of the observed variables among the four latent variables, we conducted a reliability analysis on the latent variables. The reliability test was regarded passed if the Cronbach's Alpha value and CR value in the internal consistency reliability test were both greater than 0.70. Table 5 show the initial latent variable reliability test results. The CR values of "Demographic" and "Can identify phenomenon of pareidolia" were both greater than 0.7, but the Cronbach's Alpha value failed the test, indicating that the observed variables need to be modified or removed. Observed variables with variable loadings less than 0.5 were considered to lack correlation; thus, we removed "Working" and "Education" from "Demographic" and "Pareidolia 2" and "Pareidolia 21" from "Can identify the phenomenon of pareidolia", respectively. The reliability test and correlation test of the revised PLS–SEM model both met the standard, indicating that the selected data had good reliability. Table 6 show the convergent validity test results of the PLS–SEM model. The AVE values of revised models were all greater than 0.5, suggesting that the revised model had better construct validity.

According to the PLS–SEM model path analysis in Figure 3, "Human face identification" (0.675), "Can identify phenomenon of pareidolia" (0.271), and "Demographics" (0.374) all had a direct impact on "Reaction", and "Human face identification" (0.675) witnessed the most significant influence on "Reaction"; this may infer that people found it easier to identify pareidolian faces in the built environment, resulting in stronger emotional

reactions. "Demographics" showed no significant influence on both "Human face identification" (0.163) and "Can identify the phenomenon of pareidolia" (0.256). As demographic factors such as age and lifestyle changed, "Human face identification" and "Can identify the phenomenon of pareidolia" also experienced some variations. "Can identify the phenomenon of pareidolia" had a moderate impact on "Human face identification"; this may indicate that people were more likely to identify human faces in a built environment under an illusion. As a result, "Human face identification" was the most significant factor in influencing people's reactions to the built environment. In addition, "Can identify the phenomenon of pareidolia" and "Demographics" both had a greater influence on pareidolian reaction through the indirect influence of "Human face identification".

**Table 5.** Initial latent variable reliability test results.

**Table 6.** Convergent validity test results of revised model.

**Figure 3.** SEM tests between demographic data and sensitivity to pareidolia in the built environment.

#### *3.4. Discussions*

Consistent with the results of many studies, gender influenced the sensitivity and reactions towards pareidolia in a built environment, which might be explained by Pavlova et al. [34], who claimed that females are more likely to detect faces in arrangements of food on a plate. Age had a significant influence on pareidolia sensitivity. According to the single factor analysis, respondents' sensitivity to identify faces increased with age in the 11–40 age group, and they were more likely to have pareidolia than those over 40. This might indicate that the age of 40 is a watershed for pareidolia. However, there were significant differences between males and females in the 21–40 age group but non-significant differences in other age groups. Job stress seemed to have no significant impact on sensitivity to pareidolia in the built environment, which was not fully in line with [11]. Although the between-subject effects between job stress and lifestyle demonstrated a significant influence on pareidolia identification, the experimental results revealed that this influence mainly relied on the lifestyle factor. This is probably because pareidolia is an automatic phenomenon that occurs rapidly, insensibly, and compulsorily, and people do not realize their perception is misled by pareidolia [30]. Different lifestyles had some impact on the sensitivity and reactions towards pareidolia. This could be explained by Chalup et al. [14], who claimed that neuropathological and psychiatric conditions could affect the sensitivity and reactions to identify facial expressions of sensation. In the current research, smokers, topers, and long-term medicine users tended to be more sensitive to pareidolian images in the built environment. Combining gender factors, long-term medicine users and smokers in males were more sensitive than females in face identification. Additionally, females who abused alcohol, slept late, and did not have these lifestyles were more likely to identify faces than males. This finding was consistent with previous findings that females are more likely to identify faces in pareidolian images [34] since sleeping late and having none of these lifestyles is common of most people nowadays. In addition, the between-subject effects of age and lifestyle both showed a significant influence on pareidolia. The 41–60 age group were more likely to identify faces than those in the 11–40 year age group with different lifestyles. However, long-term medicine users in different age groups were more likely to identify faces in pareidolian images. Among the respondents who slept late and did not have these lifestyles, the 41–60 age group differed significantly from those in the 31–40 age group in pareidolian face identification.

Through the path analysis of the PLS-SEM model, "Human face identification" significantly influenced people's reactions to pareidolia in the built environment, and demographic information also showed a moderate influence on pareidolian reactions. Pareidolian human faces appearing in buildings may cause positive or negative effects on human life. Some may lead to happiness and relaxation, but more lead to feelings of creepiness and depression. The facial expressions perceived in house designs through pareidolia can trigger emotional reactions based on the time of day [3].

#### **4. Conclusions**

We can categorize Pareidolia in a built environment into unintentional and intentional. The sensitivity and reactions towards pareidolia in a built environment are spontaneous and influenced by the age, gender, and lifestyles of observers. The majority of the respondents chose to change the design layout when they detected pareidolia in the built environment, including those living in the built environment with intentional pareidolian design elements. Residents rarely like any pareidolian faces in buildings, whether intentionally or unintentionally. The difference in gender influences the sensitivity and reactions towards pareidolia in a built environment because females have a stronger capability in face identification. The age of the respondents has a significant influence on pareidolia, and the age of 40 seems to be a watershed. Job stress has little impact on sensitivity to pareidolia in a built environment, but the between-subject effects of job stress and lifestyle have a considerable impact on pareidolia due to the dominance of lifestyle. The sensitivity to pareidolia varies depending on one's lifestyle. Smokers, topers, and long-term medicine users are more sensitive to pareidolian images in the built environment. Similar results are

observed even when gender and age factors are taken into account. The identification of the pareidolian face in the built environment significantly influences pareidolian reactions and sensitivity. In addition, the pareidolian human faces appearing in buildings may cause positive or negative effects on human life. Some may lead to happiness and relaxation, but more lead to creepiness and depression. The facial expressions perceived in house design through pareidolia can evoke the emotion of the observers. As a limitation, a finding that remained unexpected is that most of the pareidolian images in the built environment are much more easily detected during early morning and at midnight but remain much less able to be perceived at midday. As a result, further and more in-depth studies are required, especially multi-disciplinary studies, to address this issue satisfactorily.

**Author Contributions:** Conceptualization, C.W. and Y.M.; data analysis, L.Y.; resources, L.C.W. and C.G.; writing—original draft preparation, C.W.; writing—review and editing, Y.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by [Fujian Provincial Department of Science and Technology] grant number [2018J05121], [Fujian Provincial Department of Science and Technology] grant number [2021I0014], [Quanzhou Government] grant number [605-50X17234], [Huaqiao University] grant number [17BS201], and [Huaqiao University] grant number [605-50Y14007]. The APC was funded by [Fujian Provincial Department of Science and Technology].

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** All data, models, and code generated or used during the study appear in the submitted article.

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

#### **Appendix A. Pareidolia Questionnaire Survey**

You are invited to participate in a study examining the sensitiveness and reactions to pareidolia in a built environment. Pareidolia is a phenomenon of seeing faces in the environment, buildings and objects that surround us. Participation in this study is completely voluntary, your response will be completely confidential and only used for research purpose. In return for helping us with this, we would be very happy to give you the copies of the results of the study. For questions regarding this study, please contact Liangcheng YU at<EMAIL_ADDRESS>
#### *Appendix A.1. Background Information*

Please TICK the following items that most closely correspond to your background.

#### *Appendix A.2. Response to Pareidolia in the Built Environment*

This section is designed to ask you about your sensitiveness and reactions toward pareidolia in a built environment. There are no right or wrong answers, hence, you don't need to spend too long on each question, just go with your first impulse. Please TICK the item that that best describes your feeling and reactions.

Thank You for Your Participation.

#### **References**

### *Article* **Co-Movement between Carbon Prices and Energy Prices in Time and Frequency Domains: A Wavelet-Based Analysis for Beijing Carbon Emission Trading System**

**Rundong Luo, Yan Li \*, Zhicheng Wang and Mengjiao Sun**

Business School, Shandong University, Weihai 264209, China<EMAIL_ADDRESS>(R.L<EMAIL_ADDRESS>(Z.W<EMAIL_ADDRESS>(M.S.)

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** This study aims to investigate the co-movement and lead–lag relationship between carbon prices and energy prices in the time–frequency domain in the carbon emission trading system (ETS) of Beijing. Based on wavelet analysis method, this study examines the weekly data on oil and natural gas prices and carbon prices in Beijing ETS from its establishment in November 2013 to April 2019. Empirical results show the following important findings: (1) Carbon and natural gas prices are mainly negatively correlated, with natural gas prices occupying a leading position in the 12–20 weeks frequency band, indicating that the increase (decrease) of natural gas price will lead to the decrease (increase) of carbon price; (2) carbon and oil prices show an unstable dependence relationship, and their leadership position in the market constantly changes. The partial wavelet coherency and partial phase differences vary greatly in different time–frequency domains, indicating that there is no stable coherency between oil prices and carbon prices. The estimation results prove the existence of coherency between the carbon and energy prices in the Beijing ETS. The research findings of this paper provide quantifiable references for investors to achieve risk control in asset allocation and investment portfolio in the ETS market.

**Keywords:** Beijing carbon emission trading system; energy prices; wavelet analysis; time–frequency domain

#### **1. Introduction**

Climate change should be addressed through coordination and consummation of energy policies in the global society [1,2]. The Kyoto Protocol, which was adopted in December 1997, stipulates that by 2008–2012, total emissions of major greenhouse gases in developed countries should have been reduced by 5.2% relative to the baseline in 1990. Since the Kyoto Protocol came into effect, a carbon emission trading system (ETS) has been implemented by a growing number of countries and regions [3]. The global carbon market has emerged since the implementation of the European Union's ETS in 2005, with both carbon trading volume and turnover increasing rapidly [4].

The rationale for ETS is theoretically straightforward: Operating the carbon market encourages enterprises with low emission-reduction costs to exceed their emission reductions and sell their remaining carbon quotas or greenhouse gas emission reductions to companies with high emission-reduction costs through trading. Such a method helps companies with high emission-reduction costs to achieve the set emission-reduction targets and effectively reduces emission-reduction costs.

In the carbon market, greenhouse gas emission credits have become a tradable commodity that can be traded between regions and countries [5]. Carbon prices are jointly formed by the supply-and-demand relationship of greenhouse gas emission quotas in the carbon market [6]. The ETS embraces multiple attributes of environment, market, finance, and policy [7]. The price level itself can have varying degrees of effects on im-

**Citation:** Luo, R.; Li, Y.; Wang, Z.; Sun, M. Co-Movement between Carbon Prices and Energy Prices in Time and Frequency Domains: A Wavelet-Based Analysis for Beijing Carbon Emission Trading System. *Int. J. Environ. Res. Public Health* **2022**, *19*, 5217. https://doi.org/10.3390/ ijerph19095217

Academic Editors: Francesco Nocera, Roberto Alonso González Lezcano and Rosa Giuseppina Caponetto

Received: 28 March 2022 Accepted: 24 April 2022 Published: 25 April 2022

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

**Copyright:** © 2022 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/).

proving environmental quality, reducing energy demand, and promoting macroeconomic growth [8].

China has committed to reducing its carbon emissions from peak levels by 2030, and seven pilot carbon ETSs in Beijing, Tianjin, Shanghai, Chongqing, Hubei Province, Guangdong Province, and Shenzhen have been launched [9]. The objective is to compel companies to carry out low-carbon technology innovation through market-oriented means so that low-cost emission reduction can be achieved. The carbon reduction in each pilot has achieved initial results since the pilot initiation of the Chinese ETS. *China's Policies and Actions for Addressing Climate Change (2019)* states that the country has fulfilled its international commitment to reducing carbon emissions intensity by 40% to 45% by 2020 compared with the levels in 2005. However, China's carbon market remains troubled by inadequate development and a complex policy environment for energy saving and emission reduction [10]. In this context, improving the operating efficiency of the ETS and using market means to achieve carbon dioxide emission reduction has become a centrally important research question [11].

In this study, carbon and energy prices are dynamically analyzed by using wavelet analysis method based on the data in the Beijing carbon ETS from November 2013 to April 2019. Results reveal the coherency between the carbon and energy prices. The lead–lag relationship between the carbon and energy prices is analyzed through the lag effect of the response between changes in the carbon and energy prices.

The innovation of this study lies in the following aspects: (1) wavelet analysis, which can deal with nonstationarity in economic time series, is adopted for accurate analysis of the structural changes in carbon and energy prices in a nonstationary time series; (2) local coherency analysis is performed from different frequency scales (short and long term), and then the dynamic coherency between carbon and energy prices is studied; and (3) wavelet partial phase difference is used to analyze the lead–lag relationship amongst changes in carbon prices, natural gas prices, and oil prices.

This study adds quantifiable evidence on the relationship between carbon and energy prices in China's ETS and contributes to deep understanding of the influencing factors of carbon prices and the interaction between China's ETS and energy markets. Thus, this study provides a valuable reference for investors in asset allocation and investment portfolio for risk control.

The rest of this paper is organized as follows: Section 2 describes Beijing's carbon market, Section 3 reviews the literature, Section 4 introduces the methodology and data, Section 5 discusses empirical results, and Section 6 presents the conclusion and recommendations.

#### **2. Beijing Carbon Market**

The Beijing carbon ETS was established on 28 November 2013. According to data provided by the Beijing Carbon Exchange, by the end of 2018, the total amount of Beijing's carbon market quota reached 51 million tons, the total exchange volume was 29.07 million tons, the total trading turnover reached CNY 1.049 billion, and the average annual exchange price fluctuated around CNY 50/ton. Among the pilot carbon markets in China, the Beijing carbon market ranks first in volume, turnover, and trading activity. According to the *Beijing Municipal Commission of Development and Reform and the Beijing Environment Exchange (2018)*, the Beijing carbon market covers approximately 1000 companies, including electric power, heat, cement, petrochemicals, other industrial companies, service industries, urban rail transit, and public electric passenger cars.

The experience of carbon ETS pilots, as pioneers, plays an important role in establishing a national carbon emission market [12,13]. As one of the pilot regions of China's carbon market, Beijing ETS has achieved good performance and results in reducing emissions. According to the *Annual Report of Beijing Carbon Market 2018*, from 2013 to 2018, the energy consumption per CNY 10,000 of GDP and carbon dioxide emissions in Beijing dropped by 22.5% and 28.2%, respectively, while the city's energy efficiency ranks first in China. This finding indicates that initial results have been achieved in Beijing through the adoption

of the market means to promote energy conservation and emission reduction. In terms of market construction, a relatively complete policy and regulation system of "1 + 1 + N" has formed in the Beijing carbon exchange market. Specifically, the policies include the legislation set forth during the Beijing Municipal People's Congress on 27 December 2013 (1), local government regulations presented on 28 May 2014 (1), and the specific rules and regulations (N) of the Municipal Commission of Development and Reform and relevant departments in recent years. Meanwhile, a market supervision system based on the maximum position limit and price warning has been formed, which can be an important reference and inspiration for the construction of the national carbon market. The establishment of China's carbon market is a major impetus to further strengthen the country's ecological development and promote energy conservation and emission reduction measures.

Beijing ETS was chosen for this study for three main reasons. Firstly, the Beijing carbon market is more active, with trading entities actively trading in the market, its turnover is increasing, and the trading system is operating safely and stably [14]. The price of carbon allowances in Beijing has become an important reference for other pilot cities [15]. Secondly, Beijing carbon market information is disclosed more timely, the market carbon price can reflect market information in a more timely manner, and the construction of the carbon market system is at the forefront of the pilot carbon market [16]. Finally, Beijing's energy use efficiency is among the highest of all Chinese provinces and cities, indicating that the use of market instruments to promote energy conservation and emission reduction has begun to bear fruit in Beijing [17]. Based on these considerations, this study selects the Beijing carbon market, which has performed well in the carbon market pilots, as the subject of this study.

#### **3. Literature Review**

A growing body of literature has begun to examine carbon price fluctuations and their influencing factors [18]. The contributing factors of the carbon market can be grouped into internal factors and external ones. The internal factors mainly include the market trading behavior, the adjustment of carbon trading policy, and the uncertainty about economic policy. For example, Wang et al. [19] studied the impact of different types of trading behaviors in the EU ETS on the carbon market. It was found that "compliant transactions" and the inventories of the financial intermediaries have an impact on carbon price trends, while "non-compliant transactions" have an impact on carbon price volatility. Other researchers have confirmed that the EU and China's carbon trading policies will have an impact on the carbon market [20,21]. Li et al. [22] found that trade policy uncertainty and monetary policy uncertainty will have a positive impact on China's carbon market prices, while exchange rate policy uncertainty has a negative impact. External factors mainly refer to macroeconomic and weather variables. For example, Bento and Gianfrate [23] proved that macro-level determinants, such as gross domestic product (GDP), the presence of a national carbon price mechanism, and industry attributes, are correlated with the carbon market, and that countries with higher per-capita GDP are usually subject to higher carbon market price levels. Hintermann [24] found that in the EU carbon market, temperature and precipitation affect ETS quota prices in a nonlinear manner. In addition, some studies argue that the performance of financial market, energy market, and foreign exchange market also exert a significant effect on the volatility of carbon prices [25–27]. Sun et al. [28] studied the causal relationship between China's carbon price and the four major energy-intensive stock indexes and found that there is a weak two-way causal relationship between China's carbon market and the stock market.

Research on the coherency between the carbon market and the energy market has motivated an in-depth study of the influencing factors and effects to form a theoretical underpinning for the development and improvement of the carbon market and its efficiency [29–31]. Adekoya [32] found that the prices of energy sources such as crude oil, natural gas, and coal can help investors predict the price of carbon allowances in Europe. Energy prices, such as those of oil and natural gas, are considered amongst the important

influencing factors of the carbon price [33]. Zhang and Sun [34] proved that a one-way and time-varying spillover effect exists between the European ETS and the coal and natural gas markets, and no significant spillover effect can be observed between the European carbon market and the crude oil market. However, Ji et al. [35] studied the coherency between the returns and volatility of the carbon energy system under the EU carbon emissions exchange system. Their findings show that Brent crude oil prices play an important role in affecting carbon price fluctuation and risks. Wu, Wang, and Tian [29] proved that the volatility spillover between the coal and carbon emission markets is the strongest by using the recursive graph method and recursive quantitative analysis method. Hammoudeh et al. [36] studied the effect of changes in United States crude oil, natural gas, and other energy prices on carbon prices. Jiang and Chen [37] found that the outbreak of COVID-19 will expand the spillover effect between carbon and energy markets. To sum up, there has been growing evidence of the coherency between the carbon market and the energy market.

For China's carbon market, few studies have focused on its coherency with the energy markets. Zeng, Nan, Liu, and Chen [15] studied the dynamic relationship between carbon prices in the Beijing carbon market and economic development and energy prices through the structural vector autoregression (SVAR) model. Applying the adaptive lasso method, Guo [38] found that the carbon price in the Shenzhen ETS was mainly affected by the EUR exchange rate and domestic oil prices. Wang et al. [39] simulated the changes in China's energy system and the trend of carbon emissions by establishing a hybrid energy model to minimize costs. Dong et al. [40] studied the relationship between China's regional carbon emission intensity and energy structure by using a static spatial econometric model and a panel co-integration model. Wagner et al. [41] investigated the relationship between renewable energy costs and carbon prices in China. Using the vector error correction model, Zhou and Li [42] studied the dynamic relationship amongst energy prices, macroeconomic indicators, air quality, and carbon prices in the Hubei carbon ETS. By using a multi-time series model and an ARCH model, Li and Lei [43] found that the industrial revenue, energy prices, government intervention, and number of participating companies in the carbon ETS in Hubei affected carbon prices considerably.

Although the above studies shed light on our research, they still have shortcomings. In terms of research perspective, although the coherency between carbon and energy prices in the Chinese carbon market has been investigated, studies are mostly conducted from a static perspective without reflecting the dynamic coherency and leader–follower (lead– lag) relationship between variables in different time and frequency domains. In terms of research method, empirical research on the relationship between carbon and energy prices is mainly conducted by using geographically weighted regression, SVAR, Bayesian SVAR, and other traditional approaches [15,44,45]. These research methods have certain weakness and difficulty in accurately identifying the effective information of carbon prices and energy prices from the perspective of time and frequency domains. Multivariate wavelet analysis can effectively compensate for the inadequacies of traditional methods [46]. Wavelet analysis tools can reflect the changes in variables in the time domain and the short-term and long-term coherency between variables in the frequency domain [47]. Wavelet analysis has been used in research on carbon emissions and carbon assets [48–50]. Therefore, a wavelet analysis tool is a powerful alternative to confirm the coherency between carbon prices and energy prices in various time and frequency domains. Based on this, we use the wavelet analysis method to conduct dynamic analysis of carbon price and energy price and explore the correlation between carbon price and energy price. At the same time, the lead–lag relationship between the carbon and energy prices is analyzed through the lag effect of the response between changes in the carbon and energy prices.

#### **4. Methods and Data**

#### *4.1. Research Methods*

Wavelet analysis is a mathematical analysis method that was gradually developed in the middle and late 1980s. This method inherits the localized idea of Fourier transformation and overcomes its shortcomings in that the window size is not changed along with frequency [51]. Wavelet theory has become one of the most commonly used tools in the field of signal processing [52]. Its advantages are the following: (1) Analysis is conducted by combining the windows with different frequency fluctuations and different time scales, that is, combining time and frequency domains [53,54]. (2) The method has good localization properties and can find information about structural characteristics hidden in the data that cannot be recognized by other signal methods. (3) Nonstationary time series can be handled and analyzed well. (4)In the long wavelength band, the wavenumber resolution of the overall spectrum of the wavelet transform is higher and the error is lower [55]. On the basis of these excellent characteristics, wavelet analysis is widely used in signal processing, physics, geography, and other fields, and it began to be used in the economic and financial fields in the 1990s [56]. Continuous wavelet transform and wavelet analysis tools, such as wavelet power spectrum, wavelet coherency coefficient, and wavelet phase difference, are introduced in the present paper. In accordance with the research needs of this study, wavelet analysis tools can be extended to cases with three time variables.

#### 4.1.1. Continuous Wavelet Transform

Given a time series *x*(*t*), the expression of its continuous wavelet transform function is

$$\mathcal{W}\_{\mathbf{x}}(\mathbf{r},s) = \int\_{-\infty}^{\infty} x(t) \overline{\psi}\_{\mathbf{r},s}(t) dt. \tag{1}$$

In Equation (1), *ψτ*,*s*(*t*) indicates a complex conjugation function of the daughter wavelet function *ψτ*,*s*,

Where *ψτ*,*s*(*t*) is a sequence of functions obtained after scaling and translation of mother wavelet function *ψ*(*t*) [55,57], which can be shown as

$$
\psi\_{\mathbf{r},s}(t) = \frac{1}{\sqrt{|s|}} \psi\left(\frac{t-\tau}{s}\right), s, \tau \in \mathbb{R}, s \neq 0. \tag{2}
$$

where *s* and *τ* are the scaling and translation parameters, respectively. When |*s*| > 1, the mother wavelet is stretched and transformed to obtain a high-frequency part with a short duration in the sample time series. When |*s*| < 1, the mother wavelet is compressed and transformed to obtain a low-frequency part with a long duration in the sample time series. The translation transformation is the movement of the mother wavelet function on the time axis. Different values of *τ* represent different window positions.

The mother wavelet function *ψ*(*t*) needs to satisfy the following conditions, which means that the function fluctuates up and down around the time axis:

$$0 < \int\_{-\infty}^{+\infty} \frac{\psi(\omega)}{|\omega|} d\omega < \infty. \tag{3}$$

$$\int\_{-\infty}^{+\infty} \psi(t)dt = 0.\tag{4}$$

$$\int\_{-\infty}^{+\infty} \psi^2(t)dt = 0.\tag{5}$$

The mother wavelet function has numerous forms. In this study, Morlet wavelet is used as the mother wavelet function because it is more commonly used in the analysis of financial and economic time series. The Morlet wavelet is expressed as follows:

$$
\psi\_{\omega\_0}(t) = \pi^{-\frac{1}{4}} e^{i\omega\_0 t} e^{-\frac{t^2}{T}}.\tag{6}
$$

The value of *ω*<sup>0</sup> deserves special attention. A high value leads to poor localization property of the Morlet wavelet in the time domain, while a small value leads to the poor localization property of the Morlet wavelet in the frequency domain. To effectively balance the effect of these two aspects, *ω*<sup>0</sup> = 6 is usually considered.

#### 4.1.2. Wavelet Power Spectrum

When the wavelet power spectrum is expressed as the volatility of a single time series under the time and frequency domains, it is called the self-wavelet power spectrum and is defined as |*Wx*(*τ*,*s*)| 2 . The volatility in the self-wavelet power spectrum can be indicated by the warm and cool colours in the picture (Figures 1–3). The blue part corresponds to the low-power region, indicating that the volatility is low at this time and frequency domain; the red part represents the high-power region, indicating that the volatility is high at this time and frequency domain.

When the wavelet power spectrum is expressed as the coherency of two time series in the time and frequency domains, it is called cross-wavelet power spectrum and defined as *Wxy*(*τ*,*s*) <sup>2</sup> <sup>=</sup> <sup>|</sup>*Wx*(*τ*,*s*)<sup>|</sup> 2 *Wy*(*τ*,*s*) 2 . Similarly, the strength of the coherency between the two time series can be indicated by the warm and cool colours in the cross-wavelet power spectrum. The blue part corresponds to the low-power region, indicating that the correlation of two time series is low at this time and frequency domain. The red part represents the high-power region, indicating that the coherency of two time series is high at this time and frequency domain.

**Figure 2.** Multiple-wavelet power spectrum between Beijing carbon price and energy (natural gas and oil) prices.

**Figure 3.** Partial wavelet coherency, phase difference, and wavelet gain. **Note:** (**a**) partial wavelet coherency between Beijing carbon price and natural gas price (**a.1**) or oil price (**a.2**); (**b**) partial wavelet phase difference between Beijing carbon price and natural gas price (**b.1**) or oil price (**b.2**); and (**c**) partial wavelet gain between Beijing carbon price and natural gas price (**c.1**) or oil price (**c.2**).

#### 4.1.3. Wavelet Coherency Coefficient and Wavelet Phase Difference

Although the cross-wavelet power spectrum can obtain the correlation between the two time series in the time and frequency domain, knowing which sequence plays a leading role is difficult, that is, the lead–lag relationship between the two time series is unclear. Therefore, the introduction of wavelet coherency coefficients and wavelet phase differences is required.

The wavelet coherency coefficient, which can measure the coherency and size of two time series, is defined as

$$R\_{xy} = \frac{\left| s\left(\mathcal{W}\_{xy}\right) \right|}{\left[ s\left( \left| \mathcal{W}\_x \right|^2 \right) s \left( \left| \mathcal{W}\_y \right|^2 \right) \right]^{1/2}}. \tag{7}$$

In Equation (7), *s* is the smoothing operator achieved by convolution in both time and scale. To further obtain the lead–lag relationship of the time series, we define the wavelet phase difference as

$$\phi\_{\rm yx} = \tan^{-1} \left( \frac{\heartsuit \left( W\_{\rm xy} \right)}{\Re \left( W\_{\rm xy} \right)} \right). \tag{8}$$

In Equation (8), *Wxy* and *Wxy* represent the real and imaginary parts of crosswavelet power spectrum *Wxy*. The different values of wavelet phase difference represent the lead–lag relationship between two time series in specific time and frequency domains. If *φyx* = 0, then *x*(*t*) and *y*(*t*) of the time series are completely positively correlated. If *φyx* = ±*π*, then the time series are completely negatively correlated. If *φyx* ranges between 0 and *π*/2, then the series are positively correlated and *y*(*t*) occupies a leading position compared with *x*(*t*). If *φyx* ranges between −*π*/2 and 0, then the series are negatively correlated and *x*(*t*) occupies a leading position compared with *y*(*t*). If *φyx* ranges between *π*/2 and *π*, then the series are negatively correlated and *x*(*t*) occupies a leading position compared with *y*(*t*). If *φyx* ranges between −*π* and −*π*/2, then the series are negatively correlated and *y*(*t*) occupies a leading position compared with *x*(*t*).

Wavelet analysis can be extended to cases with multiple time series. This study focuses on three time series. The corresponding wavelet power spectrum is called multiple-wavelet power spectrum, which can be used to analyze the coherency between multiple time series. The warm and cool colors indicate the coherency between the series. Corresponding partial wavelet coherency and partial phase difference indicate research on the coherency and lead–lag relationship between the other two time series, with one time series controlled. Due to the limited space, the paper of Aguiar can be referred to for the specific content of multiple time series [58].

#### *4.2. Data*

Three variables are considered in the estimation model: carbon exchange price in the Beijing carbon ETS (BEA), natural gas price, and oil price. Carbon dioxide gas emissions originate mainly from the combustion of fossil energy. Different enterprises have different needs for various types of fossil energy as a result of changes in fossil energy prices. When the price of fossil energy is low, a company increasingly uses fossil energy. This condition leads to greater carbon dioxide emissions, which raises the demand for carbon subsidies and further increases carbon prices, and vice versa. Therefore, oil and gas prices are selected as two variables that affect carbon prices. The present study uses weekly data. Carbon price is the transaction price data on the Beijing carbon market from its establishment until 11 April 2019; these data are obtained from the China Carbon Emissions Exchange Network. Due to the high correlation between the Chinese energy market and the international energy market [59] and the lack of relevant data in China, natural gas prices adopt the natural gas spot price of Port Henry, Louisiana, and the oil price adopts the Russian ESPO crude oil spot price.

#### **5. Empirical Results**

#### *5.1. Empirical Findings of Self-Wavelet Power Spectrum*

Self-wavelet power spectrums of carbon price, natural gas price, and oil price in the Beijing carbon market are given in the following. A preliminary analysis on the characteristics of the three time series for carbon price, natural gas price, and oil price from the perspective of volatility shows the actual oscillations of each time series. On the left side of Figure 1, the monthly rate of return of the BEA price and trend chart of the monthly increase rate of natural gas and oil prices are plotted. On the right side of Figure 1, the self-wavelet power spectrum is plotted, with a time span from the end of 2013 to the 14th week of 2019.

In terms of the energy spectrum, the abscissa represents time and the ordinate represents frequency, which is represented by a time period in weeks. Frequency and period can be converted by using the following formula: frequency = 1/period. At each time and frequency, the color of the wavelet power spectrum shows the fluctuation degree of the time series. The color ranges from dark blue to deep red. The deep blue color indicates weaker volatility, while the deep red color indicates stronger volatility. The part inside the red area, that is, the part indicated by the white line, is the time–frequency domain with the most volatility. The black and grey circles represent significance levels of 5% and 10%, respectively. The black conic line indicates the cone of influence, which identifies the wavelet power spectrum region and existence of edge effects. Aguiar's wavelets toolkit [58] allows current research using wavelet analysis to proceed smoothly. The toolkit can be run with MATLAB, which facilitates the research work of scholars. We thank them for their development and sharing of wavelets tools.

The results of self-wavelet power spectrum are presented in Figure 1, which helps identify the volatility of carbon prices and natural gas and oil prices in the Beijing carbon market under different timescales. Figure 1a.1 shows that the carbon price in the Beijing carbon ETS exhibits a relatively stable fluctuation overall, while strong fluctuations occurred in several areas, which is consistent with the red region in the power spectrum in Figure 1b.1. Days with strong fluctuations are usually concentrated from mid-2014, mid-2015 to mid-2016, and from early 2018 to 2019, which represents the high-frequency domain (less than 17 weeks). In general, the carbon price in the Beijing carbon ETS has been relatively stable. The local fluctuations may be due to the following reasons: in mid-2014, the Beijing carbon ETS conducted its first annual compliance work since its establishment. During this period, emission-controlling companies that have limited knowledge of the carbon market decided to participate so that they could complete the year's emission-reduction tasks. Furthermore, the secondary market relaxed the conditions for institutional investor participation and introduced the participation of individual investors, making carbon trading active at this stage and thereby leading to large fluctuations in carbon prices. From mid-2015 to mid-2016, the Beijing carbon ETS had more abundant exchange products and initially exerted its central market function, which resulted in considerable activity. In 2018, the Beijing carbon ETS committed to promoting the formation of a multilevel carbon financial market system that coordinates onsite and offsite co-movement, which makes the carbon market active and prices volatile.

The price of natural gas always exhibits obvious cyclical fluctuations, and large fluctuations occur in the fourth quarter of each year. This condition occurs mainly because the increased heating demand, inventory data, and cold air lead to fluctuations in gas prices as the weather gradually cools down [60]. This condition corresponds to the time interval in which the red region appears in the power spectrum, and volatility is higher under the high-frequency domain (periods under 17 weeks), as shown in Figure 1b.2.

In comparison, oil prices are more susceptible to political and economic factors [61]; thus, the oil price is more volatile. Figure 1a.3 shows that the oil price was on a downward trend before 2015, whereas it rapidly increased after 2015. Since then, the price continuously declined and increased at a great amplitude. In terms of oil prices, the more volatile dates are mainly from 2015 to 2017, which are high-frequency domain with a period less than

24 weeks. This condition may have occurred because of the recovery of the world economy in 2015, which effectively promoted the recovery of the global crude oil market. The oil price increased sharply and then fluctuated sharply in the next two years. Oil prices are susceptible to demand and competition, trade prospects, and politics [62], which are also important reasons for oil price fluctuations.

#### *5.2. Empirical Findings of Multiple-Wavelet Power Spectrum*

Multiple-wavelet power spectrum is helpful in identifying the coherency of multiple time series. Figure 2 shows an estimate of the multiple coherency between the carbon price and the other two time series (natural gas price and oil price), indicating the time and frequency domains when the coherency is the strongest.

Figure 2 shows that from the latter half of 2015 to the end of 2018, the carbon ETS is related to the natural gas and oil markets. In addition, the span of each black circle line is less than half a year, which indicates that no significant long-term coherency exists between the carbon and energy markets. Furthermore, in the high-frequency domain (12–20 weeks), the carbon price and the natural gas and oil price in the Beijing carbon ETS in the two periods of the first half of 2017 and 2018 exerted a significant coherency. This coherency can be reflected in the red area distribution of the multiple power spectrum. In the low-frequency domain (34–58 weeks), two circles exist at the same time point. In general, the relationship between the carbon price and the natural gas and oil price in Beijing's carbon ETS is hardly affected by long-term and continuous changes, whereas it is more affected by short-term price shocks. The influence on carbon prices and natural gas and oil prices intensifies over time, as proven by the fact that the significantly coherent time–frequency domains in 2017 and 2018 are more concentrated.

Multiple-wavelet coherency presents evidence of coherency between carbon and energy prices (natural gas and oil) at different frequency and time domains, illustrating the importance of the relationship between carbon and energy prices (natural gas and oil). However, the relationship between the two variables cannot be distinguished by this information alone. Therefore, further analysis should be based on the partial wavelet coherency and partial phase difference.

#### *5.3. Empirical Results of Partial Wavelet Coherency and Partial Phase Differences*

This study employed partial wavelet coherency and partial phase difference to identify the lead–lag relationship of carbon and natural gas and oil prices in the Beijing carbon ETS as well as evaluate the dynamic relationship between these prices in the time–frequency domain.

Figure 3 describes partial wavelet coherency, partial wavelet phase difference, and partial wavelet gain between carbon price and each energy price (after another energy price is controlled), thereby reflecting the effect of one variable on another. The middle and right sides of Figure 3 are selected to display the wavelet phase difference and wavelet gain of two frequency bands, which are, respectively, 12–20 and 34–58 weeks, to intuitively reflect the coherency information. This finding is consistent with the coherent and statistically significant frequency bands in Figure 2.

In some partial wavelet coherency diagrams, the coherency color ranges from blue (low coherency) to red (high coherency). For wavelet phase difference, the horizontal axis represents the time scale and the vertical axis represents the phase difference, ranging from −π to +π. The research results prove several findings.

Firstly, the partial coherency (Figure 3a.1) and phase difference (Figure 3b.1) results between the carbon and natural gas prices are obtained from the top part of Figure 3. Three important regions with high coherency overlapped with the highly coherent regions estimated in Figure 2. In the low-frequency domain (34–58 weeks) from the second half of 2015 to the first quarter of 2016, the phase difference ranges between −*π*/2 and 0, indicating that carbon and natural gas prices have a positive coherency. Fluctuations in natural gas prices lead to fluctuations in carbon prices, that is, natural gas is in a leading

position. Another two regions are in the high-frequency domain (12–20 weeks), including the first half of 2017 and the second half of 2018. The phase difference ranges between *π*/2 and *π*, indicating that carbon and natural gas prices have a negative coherency. An increase in natural gas prices lead to a decrease in carbon prices. The preceding research results indicate the following: (1) In the short term, as in the first half of 2017 and the second half of 2018, the carbon price of the Beijing carbon ETS is negatively related to the price of natural gas, and changes in natural gas prices lead to changes in carbon prices. (2) Over time, the dependence between carbon and natural gas prices intensifies due to the continuous improvement of the carbon price formation mechanism of the Beijing carbon ETS. Partial wavelet gain also shows that, compared with the first half of 2017, the price gains fluctuated more in the second half of 2018, indicating that natural gas and carbon prices are greatly affected by external shocks. Prior to this, the relationship between carbon prices and natural gas prices has also been studied. Duan et al. [63] found that natural gas prices can have an opposite effect on carbon prices through two channels, the aggregated carbon demand effect and the fuel switching mechanism. Hammoudeh, Nguyen, and Sousa [45] also demonstrated that an increase in the price of natural gas reduces the purchase demand of firms and reduces carbon emissions, thus causing the price of carbon to fall. However, Batten et al. [64] show that the effect of natural gas prices on carbon prices is not significant. Our findings confirm the conclusions of the first two studies.

Secondly, the partial coherency (Figure 3a.2) and phase difference (Figure 3b.2) results between the carbon and oil prices are obtained from the bottom part of Figure 3. The regional distribution of this coherency is scattered, indicating that only a short-term co-movement exists between the carbon and oil markets. Therefore, the co-movement relationship between the two is more affected by short-term shocks than by long-term persistence. To further explain the relationship between carbon prices and oil prices, three important regions are selected for further analysis in the high-frequency (period of 12–20 weeks) and low-frequency (period of 34–58 weeks) domains. These three regions have a high coherency at 10% significance level, as shown in Figures 2 and 3. The first region lies in the low-frequency domain from the fourth quarter of 2015 to the first quarter of 2016, whose phase difference ranges between −*π*/2 and 0. This result indicates that carbon and natural gas prices have a positive coherency. Oil is in a leading position in that fluctuations in oil prices lead to fluctuations in carbon prices. In the first quarter of 2018, the second region lies in the high-frequency domain, where the phase difference ranges between *π*/2 and *π*. This result indicates that carbon and natural gas prices have a positive coherency. Oil is in a leading position in that fluctuations in oil prices lead to fluctuations in carbon prices. In the fourth quarter of 2018, the third region lies in the high-frequency domain, where phase difference ranges between −*π* and −*π*/2. This finding indicates that carbon and oil prices have a negative coherency. Carbon is in a leading position in that fluctuations in carbon prices lead to fluctuations in oil prices. The analysis results based on the aforementioned three highly coherent regions show that the carbon and oil prices have an unstable dependence relationship, and their leadership position in the market constantly changes. This result is consistent with the complexity of oil price changes and the fact that the carbon price formation mechanism of the Beijing ETS is incomplete. In addition, partial wavelet gain graph shows that partial gain is stable in both frequency bands, and the value is close to 0.1, which is approximately twice the partial gain between carbon and natural gas prices. The relationship between oil and carbon prices has also been explored by previous researchers. Duan, Ren, Shi, Mishra, and Yan [63] found that oil prices can positively affect the price of natural gas through aggregated carbon demand effects. Hammoudeh, Nguyen, and Sousa [45] found that an increase in crude oil price promotes carbon price in the short run but suppresses it in the long run. In the short run, our findings are consistent with previous studies. However, in the longer term, we believe that the dominant relationship between the price of carbon and the price of oil is changing and there is no clear relationship.

#### **6. Conclusions and Policy Recommendations**

Understanding the influencing factors of carbon price changes is crucial for studying carbon market prices. In this study, empirical analysis of the coherency and lead–lag relationship between carbon prices and energy prices (natural gas and oil) of the Beijing carbon ETS in the time and frequency domains is performed by using wavelet analysis tools, including partial wavelet coherency and partial phase difference. The conclusions help carbon market participants make appropriate investment and portfolio decisions based on the relationship between the carbon and energy markets. Moreover, the conclusions also have important policy implications for market regulators. For example, risk control and management of market exchanges should be strengthened, and reasonable exchange rules, such as price restrictions and risk reserves, should be formulated to ensure the security and stability of market transactions.

The analysis reveals the following: (1) Multiple-wavelet coherency results show that the carbon price and energy price series are in the high-frequency domain around 2017 and from 2018 to early 2019 (less than 20 weeks), indicating a significantly strong coherency. However, the low-frequency domain (more than 34 weeks) shows a significant coherency in the local period from the second half of 2015 to the middle of 2016. Therefore, the coherency between carbon prices and energy price series has strong coherency in the short term and weak coherency in the long term. The estimation results indicate the correlation between the carbon market and the energy market, especially in the short term. (2) Partial wavelet coherency results of carbon and natural gas prices reflect that these prices have a significantly high coherency in the short term (specifically, the first half of 2017 and the second half of 2018). In the high-frequency band, the phase difference ranges between *π*/2 and *π* more than half the time, indicating that carbon and natural gas prices have a negative coherency, in that changes in natural gas prices lead to changes in carbon prices. (3) The results of partial coherency and partial phase difference between the carbon price and the oil price series vary greatly in different time–frequency domains. The lead–lag relationship between the two in the high-frequency and low-frequency bands also changes constantly. Therefore, the carbon and oil prices show an unstable dependence relationship. In conclusion, these prices in the Beijing carbon market have a close relationship.

Based on the research results, the following recommendations are presented:

fluctuations have emerged. Therefore, regulators should carry out effective risk control by establishing a sound market supervision mechanism so that the carbon emissions trading market can effectively achieve energy saving and emission reduction.

(4) Governments need to ensure price stability and secure supply of energy. Our findings show that energy markets are closely linked to carbon markets, and that stability in energy markets helps to ensure the effective functioning of carbon markets. Therefore, in order to fulfil the role of carbon markets in environmental sustainability, it is necessary to ensure a stable and secure supply of energy and to continue to promote the use of renewable energy.

Future research can be further explored in terms of both the research object and the factors influencing the carbon market. In terms of the research object, the official opening of China's national carbon emissions trading market was in July 2021. The relationship between the national carbon market price and the energy price can be further explored by taking the national carbon market as the research object. In terms of carbon market influencing factors, researchers could also consider the linkages between other energy markets and the carbon market, such as the renewable energy market. This could provide more reference for government policy-making and investors' investment decisions.

**Author Contributions:** Conceptualization, R.L. and Y.L.; methodology, R.L., Y.L. and Z.W.; software, Z.W. and M.S.; validation, Z.W. and M.S.; writing—original draft, Z.W. and M.S.; writing—review and editing, Y.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Science Foundation of China (NO. 71904009), Beijing Social Science Fund NO. 20GLC054, MOE (Ministry of Education in China) Project of Humanities and Social Sciences (NO. 18YJC840041), the Postdoctoral Science Foundation of China NO. 2019M652415, Young Scholars Program of Shandong University at Weihai, and Scientific Research Fund for Graduate Students of Business School of Shandong University.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Thank Sun Zuoren for the polishing service provided during the modification process.

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

#### **References**

### *Article* **Health-Related Quality of Life According to Sociodemographic Characteristics in the South Korean Population**

**Chan-Hee Park 1,†, Eunhee Park 2,3,†, Hyun-Min Oh 1, Su-Jin Lee 2, Sun-Hee Park <sup>4</sup> and Tae-Du Jung 2,3,\***

**Abstract:** Health-related quality of life (HRQoL) concerns satisfaction with life and happiness with regard to physical, mental, and social factors. RAND-36 is a publicly available, self-administered questionnaire that examines eight health dimensions. This study evaluated the HRQoL of the South Korean population using the RAND-36 questionnaire and compared HRQoL across sociodemographic characteristics. From May 2015 to May 2019, South Koreans who visited public places aged 19–80 years were recruited and the RAND-36 questionnaire was administered. Overall, 1002 participants were recruited (mean age 45.34 years, 52% men). Men scored better than women in both physical and mental health (*p* < 0.05). There were significant differences in bodily pain (*p* < 0.05), general health perception (*p* < 0.05), and energy/fatigue (*p* < 0.05) dimensions according to the participants' health condition. The HRQoL of South Koreans was lower than average in most dimensions compared with other countries. As the first study to assess this, its data can be used in future studies that apply RAND-36 to evaluate the HRQoL of diseased individuals, as they can compare their findings with those of our study population.

**Keywords:** health-related quality of life; RAND-36; South Korea

#### **1. Introduction**

The concept of health continues to change as medicine and the medical sciences develop. Health is defined as a dynamic state of human welfare characterized by physical, mental, and social potential that meets life-related needs as well as the absence of illness or infirmity. Health-related quality of life (HRQoL) relates to individuals' life satisfaction or happiness with regard to their health status and social and cultural context [1]. This HRQoL mainly addresses public health issues in high-income and rapidly aging countries. HRQoL is typically measured using generic or disease-specific tools [2]. The decision to use a generic or disease-specific tool is determined by several factors, such as the reason for the measurement, the efficiency of the applicable tools, and the population in question. Generic measures are not specific to any disease, age, or treatment group and can be used to compare patients in different situations and with the general population [2,3].

General population data play an important role in determining whether groups or individuals' scores are higher or lower than the average based on country of residence, age, and sex. Published general population data are available for several countries around the world, including the United States [4], Sweden [5], China [1,6], Hong Kong [7], Taiwan [8], and Japan [9]. However, to our knowledge, although HRQoL data of the general population

**Citation:** Park, C.-H.; Park, E.; Oh, H.-M.; Lee, S.-J.; Park, S.-H.; Jung, T.-D. Health-Related Quality of Life According to Sociodemographic Characteristics in the South Korean Population. *Int. J. Environ. Res. Public Health* **2022**, *19*, 5223. https:// doi.org/10.3390/ijerph19095223

Academic Editors: Roberto Alonso González Lezcano, Francesco Nocera and Rosa Giuseppina Caponetto

Received: 22 March 2022 Accepted: 24 April 2022 Published: 25 April 2022

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

**Copyright:** © 2022 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/).

475

are very important, there are no suitable general population data on HRQoL in South Korea. In particular, in a revolutionary dynamic situation such as the coronavirus pandemic period, there is a need for research on the HRQoL of the general population to compare between general populations and coronavirus survivors, or to confirm the change in HRQoL before and after coronavirus treatments. It is hoped that this study may be helpful for follow-up studies by presenting results prior to the coronavirus pandemic period, which will be the subject of comparison for future studies.

Many generic questionnaires have been developed to evaluate HRQoL; one of the most widely used instruments is the Short Form-36 Health Survey (SF-36) [1]. The SF-36 was created in the Medical Outcome Study (MOS), a four-year study that specifically focused on outcomes of care [10]. Since its creation, many variations and derivatives of the SF-36 have been developed. In particular, working with the RAND Corporation in 1992, Ware and Sherbourne published their version of SF-36, known as Ware-36, which focuses on eight health domains [11,12]. In 1993, Hays, Sherbourne, and Mazel released the RAND-36 [13]. Unlike Ware-36, the RAND-36 questionnaire and scoring methods are publicly available on the RAND Corporation's website. A longitudinal study performed by MOS showed that the difference in scores between the Ware-36 and RAND-36 scales was subtle, without any significant differences [13]. Specifically, the two surveys differed only in terms of bodily pain and general health dimension subscales. However, according to Hays et al., bodily pain and general health dimension scores had a correlation coefficient of 0.99 [13].

Considering this, we aimed to obtain HRQoL-related data for South Korea using the RAND-36 survey and identify the sociodemographic factors that influence HRQoL.

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

#### *2.1. Data Collection*

This retrospective cross-sectional study was conducted in South Korea to evaluate HRQoL in five urban and suburban areas. A total of 1002 South Korean citizens aged 19–80 years were recruited to visit public places where areas had foot traffic, including community centers and train stations that attract large populations. This study used a faceto-face questionnaire in Korean. Informed consent was obtained from all of the participants involved in this study. The inclusion criteria were those 19 years of age or older who are capable of independent daily living. From May 2015 to May 2019, a total of 1002 South Koreans were recruited. The participants were asked to fill out a questionnaire that included the RAND-36 health survey and items designed to obtain their general information, such as age, sex, region of residence, education status, occupational status, and marital status. In addition, respondents were asked whether they had been diagnosed with the following diseases: hypertension, angina, myocardial infarction (MI), diabetes mellitus (DM), or cancer and whether they had any of the following health problems: allergies, back pain, visual impairment, skin problems, chronic lung problems, hearing problems, functional impairment in a leg or arm, or joint pain.

Responses to items concerning diseases and current health problems were classified into four categories: (1) no diseases and current health problems, (2) past or current diseases only, (3) current health problems only, and (4) diseases and current health problems [14]. The respondents took approximately 10–15 min to complete the questionnaire, immediately checked the questionnaire after completion, and if data were missing, they were asked to enter the data in question. If there were missing data even after this, the questionnaires were excluded from the statistical analysis. Ethical approval was obtained from the Institutional Review Board of Kyungpook National University Chilgok Hospital (KNUCH IRB no. 2018-11-020).

#### *2.2. RAND-36*

The RAND-36 is a standardized questionnaire, perhaps the most widely used healthrelated quality of life (HRQoL) survey tool in the world today [15]. It includes eight health dimensions: physical function (PF), role limitations due to physical health problems (RP), bodily pain (BP), general health perception (GHP), energy/fatigue (EF), social function (SF), role limitations due to emotional problems (RE), and emotional well-being (EWB). According to the instructions provided on the RAND Corporation website, the eight health dimensions were calculated using 36 items (https://www.rand.org/health-care/surveys\_ tools/mos/36-item-short-form.html, accessed on 5 May 2015).

The items evaluating role limitations (physical and emotional) are answered using "yes" or "no" responses, whereas the other items are scored using scales ranging from 3 to 6 points. According to RAND-36 scoring algorithms, the raw scores for each item should be summed and converted to a 0–100 scale, with 100 representing the best possible health state [13]. The RAND-36 dimensions can be classified into two categories: physical health composites (PHCs), which relate to physical function, and mental health composites (MHCs), which relate to mental and emotional well-being. PHCs were obtained by summing the scores for the PF, RP, BP, and GHP dimensions, whereas MHC was determined by summing the scores for the EF, SF, RE, and EWB dimensions [13].

#### *2.3. Statistical Analysis*

Descriptive statistics are presented as mean values and standard deviations (SDs) and included percentages for absolute counts regarding age, sex, education level, marital status, and health condition. The internal consistency of the RAND-36 questionnaire was evaluated using Cronbach's α coefficient; a Cronbach's α of 0.7 or greater is generally considered to indicate adequate internal consistency [16]. In addition, construct validity was evaluated using correlation analysis, and factor loadings greater than 0.50 were considered to indicate validity [17]. Finally, analysis of covariance was used to assess the independent effects of education, marital status, and health condition on the RAND-36 scale by controlling for age and sex effects. The data were processed and analyzed using a statistical software program (IBM SPSS Statistics 25.0, IBM, New York, NY, USA).

#### **3. Results**

#### *3.1. Sample Characteristic*

A questionnaire survey was administered between May 2015 and May 2019. Data from 1002 participants (525 men and 477 women) were analyzed. All subjects were aged between 19 and 80 years, with an average age of 45.34 ± 12.00 years. Respondents' characteristics are listed in Table 1.

**Table 1.** The demographic characteristics of the respondents (n = 1002).

Of the 1002 subjects, 26% had a past or current disease, with hypertension being the most common (160; 63%), followed by MI (16; 6%), angina (10; 4%), DM (65; 26%), and cancer (61; 24%). Approximately 30% had a current health problem, with chronic back pain (115; 39%) and chronic allergies (89; 30%) being the most common.

#### *3.2. Data Quality*

The internal reliability of the RAND-36 was assessed using Cronbach's α coefficient. This indicated that six of the eight dimensions had good internal consistency, with the exception of EF and EWB. Overall, of the eight dimensions, the lowest Cronbach's α coefficient was found for the EWB (0.366), whereas the SF scored the highest (0.989; Table 2).

**Table 2.** Reliability estimates (Cronbach's α) and correlation for the RAND-36 dimensions.

PF, physical functioning; RP, role limitation, physical; BP, bodily pain; GHP, general health perception; EF, energy/fatigue; SF, social function; RE, role limitation, emotional; EWB, emotional well-being; PHC, physical health composites; MHC, mental health composites.

Spearman's correlation analysis was used to analyze the correlation between the dimensions and items. PF, RP, BP, and GHP were more closely related to the PHC scale, whereas EF, SF, RE, and EWB were more closely related to MHC. Moreover, the correlation between PHC and its constituent dimensions ranged from 0.577 (PF) to 0.710 (GHP), whereas the correlation between MHC and its component dimensions ranged from 0.580 (RE) to 0.875 (EWB; Table 2).

#### *3.3. HRQoL Data*

The HRQoL data stratified by age and sex are presented in Table 3. For each dimension, when all age groups' scores were considered cumulatively, men scored higher than women (*p* < 0.05). In this regard, the difference in the RE dimension was greatest at 4.35, and the difference in the EWB dimension was smallest at 1.96 points on a scale ranging from 0 to 100. Comparing specific age groups across sexes, females had lower scores than males in almost all dimensions.

Among males, subjects aged 50–59 years showed the lowest scores in PF, GHP, and EWB, and those aged 40–49 years showed the lowest scores in RP, BP, and RE.

Among females, subjects over 60 years old returned the lowest scores in PF, RP, BP, GHP, and EF, whereas subjects aged 20–29 years showed the lowest scores in the EWB dimension. For the MHC scale, subjects aged 30–39 years returned the highest scores, whereas subjects aged 40–49 years showed the lowest scores. For the PHC scales, the lowest scores were provided by participants over 60 years of age; the higher the age, the lower the score.

**Table 3.** Mean RAND-36 dimensions scores (SD) by gender and age groups (higher scores indicate better health).

PF, physical functioning; RP, role limitation, physical; BP, bodily pain; GHP, general health perception; EF, energy/fatigue; SF, social function; RE, role limitation, emotional; EWB, emotional well-being; PHC, physical health composites; MHC, mental health composites.

When age and sex were adjusted for, marital status was found to have no significant effect on any dimension (Table 4). However, when comparing raw scores, widow/widower subjects (15; 1%) showed the highest scores in all dimensions except for PF, in which divorced/separated subjects (11; 1%) showed the highest score. However, divorced/separated participants had the lowest scores in all dimensions except PF, RP, and EWB. Education level had no significant effect on any of the dimensions. In all dimensions, subjects who only had elementary school education (45; 4%) showed the lowest scores; however, there was no corresponding pattern among subjects with high education levels (i.e., they did not return the highest scores in all dimensions).

For BP, GHP, and EF, after controlling for the effects of age and sex, we found statistical differences in health conditions (Table 4). As expected, subjects without diseases and current health problems had the highest scores, whereas respondents with no diseases but current health problems returned the lowest scores. Health conditions had no significant effect on the MHC scale but had a significant effect on the PHC scale, which seems to reflect the results for BP and GHP.

**Table 4.** RAND-36 dimensions scores in relation to marital status, education level, and health condition (raw scores and controlled for age and gender by ANCOVAs (=in bold)).

PF, physical functioning; RP, role limitation, physical; BP, bodily pain; GHP, general health perception; EF, energy/fatigue; SF, social function; RE, role limitation, emotional; EWB, emotional well-being; PHC, physical health composites; MHC, mental health composites. \* *p*-Value < 0.05, *p*-Value apply for adjusted scores for age and gender.

#### **4. Discussion**

This study observed HRQoL-related data for South Korea using the RAND-36 health questionnaire. Such data can provide a basis for future studies seeking to analyze and evaluate the quality of life of South Koreans with regard to age and sex. Additionally, the HRQoL data can be used as comparative data for future studies that seek to use RAND-36 to evaluate the HRQoL of patients with specific diseases.

To perform this, internal consistency was evaluated using Cronbach's α coefficient, and we consequently determined that six of RAND-36's eight dimensions had good internal consistency, with EF and EWB being the exceptions. Specifically, the Cronbach's α coefficients for the EWB and EF were 0.366 and 0.369, respectively. This low coefficient can be attributed to the characteristics of South Koreans. Traditionally, South Koreans are reluctant to reveal the state of their mental health [18,19]. For example, among South Koreans, many mental health patients choose to keep their treatment and condition secret. Socially, the country has a culture that regards the diagnoses of mental illnesses as different from those of physical illnesses. For items 24, 25, and 28, which are part of the EWB dimension and concern the negative aspects of mental health, the respondents may have selected responses that portrayed their mental health more positively, regardless of their actual condition. In contrast, for items with positive content, such as items 26 and 30, the respondents may have selected responses reflecting their reality. For EF, participants may have given more positive responses to items 29 and 31 and may have answered positive items, such as 23 and 27, more accurately. This should be understood as a characteristic of regions with significant differences in the diagnosis rate of depression, including East Asia, or regions, such as South Korea, with low access to mental health services for various

reasons (social awareness, institutional disadvantage, treatment costs, drug side effects). These characteristics should be considered when interpreting the study results.

Our results showed that the RAND-36 questionnaire has reliable construct validity, which is consistent with the findings of previous surveys [20]. Our correlation analysis indicated that each of the eight dimensions was highly correlated with the related composite scale. Specifically, we found that PF, RP, BP, and GHP had higher correlations with PHC, whereas EF, SF, RE, and EWB had higher correlations with MHC, which is consistent with previous studies [11,13,20]. Thus, using RAND-36 to evaluate the HRQoL of the general South Korean population was determined to be reasonable.

Previous studies have shown that men generally scored higher than women in terms of quality of life index by age group. In particular, the higher the age group, the larger the gap, which shows that there are time and economic reasons behind the decline in quality of life as people age [21,22]. It turns out that time factors are at play in middle age and economic factors are at play in old age.

In this study, quality of life was divided into PHC and MHC groups. For the PHC dimensions, male respondents scored higher than female respondents. For females, physical scores tended to decrease with age, which can be attributed to age-related biological changes, and there was no change in their physical role or position despite being older. In males, physical scores decreased up to the 50–59 years age group and then increased, with better health status observed among older individuals. This may be due to sociological changes that occur in men when they retire in their early 60s. At that time, they can pay more attention to their health because of the economic and time allowance.

Regarding MHC dimensions, as with the PHC scales, males scored higher than females. Men had the poorest mental health at 30–39 years, when they were likely to be in the most stressful situations. At this age, they must find jobs and struggle to survive in employment. They also establish new homes and support themselves, both economically and mentally. The finding that older individuals had the best mental health was attributed to the abovementioned sociological changes that occur at retirement. Sprangers and Schwartz [23] suggested that men over 60 years of age scored higher on the PHC and MHC scales as a result of the theoretical model of response shift.  At this age, middle-aged women may have a very poor quality of life, not only because of family problems but also because of children's educational problems.

Education level and marital status had no significant effects on any dimension. In South Korea, education up to middle school is compulsory, but South Koreans are generally zealous about educating their children, with a high school enrollment rate of 99% and university enrollment rate of 69.76%, which is the highest among OECD countries [24,25]. Among the participants in this survey, 91% had graduated from high school. The fact that marital status did not have a significant effect was attributed to the Confucian culture that persists in South Korea. In South Korea, divorce and separation are considered sinful and worthy of social punishments. In addition, there are many cases where couples refrain from divorce or separation to avoid negatively affecting their children. Reflecting on these Confucian characteristics, 97% of the respondents in this study were single or married/cohabitant. These social and cultural characteristics significantly affected education level and marital status.

This study had some limitations. First of all, we did not conduct sampling before recruiting our study population. Random sampling is the best method for evaluating the quality of life in participants. However, this was not done in our study, and there is a possibility of sampling bias. Second, we recruited participants and interviewed them face-to-face in the foot traffic areas that attract large populations, so it took a long time, namely four years. This could cause a source of bias. Third, among the total respondents, the number of respondents in the widow/widower and divorced/separated groups in terms of marital status was very small, and the number of respondents in the elementary and middle school education groups in terms of education level was very small. This

was due to socioeconomic and cultural differences, which may have affected the statistical analysis. To recruit respondents from various areas, questionnaires were distributed to public and commercial service organizations, community centers, and train stations that attract large populations. Among the respondents recruited, some lived in suburban areas; however, as most of the areas were located downtown, we may have recruited more urban-based than suburban-based people, which may have affected the results. Although the interviewers who instructed the participants to complete the questionnaire received identical instructions, there could have been differences in the explanations provided by the interviewers, which might have influenced the results. It is difficult to evaluate this quantitatively; however, it can still be considered a limitation of this study.

#### **5. Conclusions**

In conclusion, the HRQoL data for the South Korean population showed that the scores were somewhat lower than those of other countries or cities. In particular, South Korea's EF and EWB dimensions, which are part of the MHC, were the sixth lowest among the countries and cities analyzed. There were no significant differences in HRQoL scores with regard to education level and marital status in South Korea, but there were significant differences in BP, GHP, and EF with regard to health conditions. These HRQoL data for our study population can be used in future studies that seek to use the RAND-36 questionnaire to evaluate the HRQoL of patients with specific diseases and compare the results with those of our study population. In addition, this is a kind of preliminary study in a study targeting a general population that has been subjected to appropriate random sampling in the future, and it is expected that it will be helpful in interpreting the results.

**Author Contributions:** Conceptualization, S.-H.P. and T.-D.J.; methodology, S.-J.L.; formal analysis, H.-M.O.; investigation, S.-J.L.; writing—original draft preparation, E.P.; writing—review and editing, C.-H.P.; supervision, S.-H.P. and T.-D.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** There is no funding to report for this submission.

**Institutional Review Board Statement:** The study was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board of Kyungpook National University Chilgok Hospitial (No. 2018-11-020).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** Available upon reasonable request.

**Acknowledgments:** We thank all the study respondents for their participation in our survey. We thank the anonymous reviewer for their helpful comments, which improved our manuscript.

**Conflicts of Interest:** There are no conflict of interest for any of the authors.

#### **References**

### *Article* **An Analysis of the Colony Structure of Prokaryotes in the Jialing River Waters in Chongqing**

**Maolan Zhang 1,†, Guoming Zeng 1,†, Dong Liang 1, Yiran Xu 1, Yan Li 2, Xin Huang 1, Yonggang Ma 2, Fei Wang 1, Chenhui Liao 1, Cheng Tang 1, Hong Li 3,\*, Yunzhu Pan 4,\* and Da Sun 5,\***

**Abstract:** At present, research on the influence of human activities (especially urbanization) on the microbial diversity, structural composition, and spatial distribution of rivers is limited. In this paper, to explore the prokaryotic community structure and the relationship between the community and environmental factors in the Jialing River Basin of Chongqing, so as to provide a basis for monitoring microorganisms in the watershed. The V3–V4 region of the 16 S rRNA gene was analyzed by high-throughput sequencing and the microbial community of the waters of the Jialing River was analyzed for the diversity and composition of the prokaryotic community as well as the species difference of four samples and correlations with environmental factors. The main results of this study were as follows: (1) The diversity index showed that there were significant differences in the biodiversity among the four regions. At the genus level, *Limnohabitans*, *unclassified\_f\_Comamonadaceae*, and *Hgcl\_clade* were the main dominant flora with a high abundance and evenness. (2) A Kruskal– Wallis H test was used to analyze the differences of species composition among the communities and the following conclusions were drawn: each group contained a relatively high abundance of *Limnohabitans*; the Shapingba District had a higher abundance of *Limnohabitans*, the Hechuan District had a wide range of *unclassified\_f\_Comamonadaceae*, and the Beibei District had a higher *Hgcl\_clade*. (3) Through the determination of the physical and chemical indicators of the water—namely, total nitrogen, total phosphorus, chemical oxygen demand, chlorophyll A, and an analysis by an RDA diagram, the results demonstrated that the distribution of microbial colonies was significantly affected by the environmental factors of the water. Chemical oxygen demand and ammonia nitrogen had a significant influence on the distribution of the colonies. Different biological colonies were also affected by different environmental factors.

**Keywords:** prokaryotes; high throughput sequencing; biological information analysis; community structure

#### **1. Introduction**

Microbial communities in aquatic ecosystems play an important role in degrading organic matter and participating in nutrient cycling, but they are sensitive to the external

**Citation:** Zhang, M.; Zeng, G.; Liang, D.; Xu, Y.; Li, Y.; Huang, X.; Ma, Y.; Wang, F.; Liao, C.; Tang, C.; et al. An Analysis of the Colony Structure of Prokaryotes in the Jialing River Waters in Chongqing. *Int. J. Environ. Res. Public Health* **2022**, *19*, 5525. https://doi.org/10.3390/ ijerph19095525

Academic Editors: Roberto Alonso González Lezcano, Francesco Nocera, Rosa Giuseppina Caponetto and Paul B. Tchounwou

Received: 29 March 2022 Accepted: 29 April 2022 Published: 2 May 2022

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

**Copyright:** © 2022 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/).

environment [1,2]. At present, the research on bacterial community diversity and its relationship with environmental factors is growing [3]. However, most studies focus on coastal areas, lakes, and oceans [4,5], and there are few studies on rivers. Due to the increasingly serious pollution of urban rivers, the response of aquatic biota in urban rivers has gradually attracted researcher's attention.

Rivers are an important part of the global water cycle and one of the important water sources for human industrial and agricultural production as well as for life; they can provide an early warning to the overall quality of the water environment [6]. The Jialing River Basin in Chongqing is an ecological security barrier in the southwest of the country. The Jialing River is an important water source for industrial and agricultural production as well as domestic water and fishery in the urban area of Chongqing. It is also an important water source conservation area in the upper reaches of the Yangtze River. The water quality of the Jialing River affects the economic and social development of Chongqing [7]. Therefore, it is significant to study the characteristics of flora in the Jialing River Basin.

Microorganisms are an important part of ecosystems and are highly adaptable and prone to variation. Their structures and functions follow the changes of environmental factors. Microorganisms can maintain the overall stability of the ecosystem [8,9] and play a vital role in water quality control [10]. They have significant ecological functions in nutrient cycling, pollution reduction, soil protection, climate regulation, and water conservation [11]. As a key part of the composition of the biosphere, bacteria play an extraordinary role in the flow of matter and energy. The compositional diversity, distribution characteristics and inter-community structure of bacterial communities in lake ecosystems are the basis and prerequisite for understanding the structure and maintenance mechanism of lake ecosystems, and provide an important basis for basic theoretical research and environmental monitoring of lake microbial ecology.

In this paper, the research method based on high-throughput sequencing was used to conduct sampling at four sampling points in the Jialing River Basin, the Hechuan District, the Beibei District, the Shapingba District, and the Jiangbei District. The total microbial DNA was extracted from the samples, a PCR library was constructed, and the V3–V4 of the fragment 16S rDNA fragment was used as the research object. The structural diversity of bacterial flora in Jialing River water was analyzed by high-throughput sequencing, their relationship with environmental factors was analyzed by redundancy analysis (RDA) [12], and the biological community was analyzed by bioinformatics and metagenomics methods. To explore the prokaryotic community structure and the relationship between the community and environmental factors in the upper, middle, and lower reaches of the Jialing River Basin, so as to provide a theoretical basis for microbial monitoring in the entire river basin in Chongqing.

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

#### *2.1. Experimental Materials*

A bacterial genomic DNA extraction kit (Baimaike Biotechnology Co., Ltd., Beijing, China) and a PCR kit (Baimaike Biotechnology Co., Ltd., Beijing, China) were used in addition to sulfuric acid, nitric acid, potassium persulfate, phenolphthalein, potassium dihydrogen phosphate, sodium hydroxide, potassium iodide, potassium dichromate, potassium dihydrogen phosphate, and acetone (Meiji Biotechnology Co., Ltd., Shanghai, China). All chemical reagents of absolute ethanol, boric acid, hydrochloric acid, and potassium phthalate (Meiji Biotechnology Co., Ltd., Shanghai, China) are of analytical grade. A vertical pressure steam sterilizer (YXQ-LS-50A), UV-Vis spectrophotometer (TU-1901, T68), and small Benchtop Refrigerated Centrifuge (Centrifμge 5418R) were also used in the experiments.

#### *2.2. Sample Collection*

The sample collection took place between July–September 2020, as this was the highwater season of the Jialing River. To analyze the characteristics of microbial colonies in the Chongqing area of the Jialing River, four sampling points (as shown in Figure 1) were

selected for sampling. The four sampling points were the Nanping Bridge of the Jialing River in the Hechuan District (upstream, indicated as H), the Beibei Jialing River Bridge (upper and middle reaches, indicated as B), the Shapingba Shimen Bridge (middle and lower reaches, indicated as S), and the Jiangbei QianSimen Bridge (downstream, indicated as J). Each sampling point is divided into upper, middle, and lower vertical levels. Due to the large variation in the river water depth, the upper sampling depth (height from the river water surface) was approximately 0.5 m, the middle layer was approximately 3–5 m, and the lower layer was approximately 5–10 m. The tertiary water was then mixed and repeated three times for each sampling point. A total of 1 L of the water sample was filtered through a 0.22 μm microporous membrane and stored at −80 ◦C. This was then sent to a sequencing company for 16S rDNA transcriptome sequencing.

**Figure 1.** Sample collection site, 1:100,000.

#### *2.3. Determination of Physical and Chemical Indicators of Water Bodies*

Four sampling points were in the Jialing River Basin, the Hechuan District, the Beibei District, the Shapingba District, and the Jiangbei District. The reason for choosing Shapingba District and Jiangbei District is that they belong to the upper end of the Jialing River and have a large population base. Beibei District belongs to the middle section of the Jiangling River and has a medium population. Hechuan District belongs to the lower end of the Jialing River, which is a suburb and has a small population. The determination of physical and chemical indicators of the samples obtained, water quality total nitrogen (GB-11894-1989), total phosphorus (GB-11893-1989), ammonia nitrogen (GB-7479-87), chemical oxygen demand (GB-11914-1989), and chlorophyll A (HJ-897-2017) were performed in accordance with the national standards.

#### *2.4. Bacterial Genomic DNA Extraction*

The total bacterial DNA extraction was performed according to the Bacterial Genomic DNA Extraction Kit Manual in Table 1. Electrophoresis of 1% agarose gel was then used to detect the amount of DNA extracted, and NanoDrop2000 (micro-ultraviolet spectrophotometer) was used to determine the concentration and purity of DNA. A total of 100 μL of buffer B1 was added to 0.5 mL of a bacterial solution, and then 100 μL of buffer B2 were added. This mixture was well-shaken, placed in a centrifuge, and centrifuged at

12,000 rpm/min for 2 min. After centrifugation, 100 μL of the supernatant were pipetted into another clean centrifuge tube for use as a template.

**Table 1.** DNA extraction kit product content.

#### *2.5. PCR Amplification*

Bacterial 16S rDNA PCR operation steps, establishment of a PCR reaction system, and 20 μL system in Table 2. After adding TransStart FastPfu buffer, dNTPs, upstream primers, downstream primers, TransStart FastPfu DNA polymerase, and template DNA were added to produce a 20 μL mixture. All reagents were collected at the bottom of the tube by brief centrifugation. The 16S rRNA gene V3-V4 variable region was PCR amplified using 338F(5 - ACTCCTACGGGAGGCAGCAG-3 ) and 806R (5 -GGACTACHVGGGTWTCTAAT-3 ) in triplicate per sample.

**Table 2.** PCR reaction system.

#### *2.6. Data Statistical Analysis*

The bases were filtered with a tail quality value below 20 of the reads, a 50 bp window was set, and the reads were filtered below 50 bp after quality control. The reads containing N bases were then removed, according to the overlap relationship between the PE reads. The paired reads were then merged into one sequence; the minimum overlap length was 10 bp, and the maximum mismatch ratio allowed in the overlap region of the spliced sequence was 0.2. The samples were distinguished according to the barcodes and primers at the beginning and end of the sequence, and the sequence direction was adjusted. The number of mismatches allowed by a barcode was 0, and the maximum number of primer mismatches was 2. Using UPARSE software (version 7.1) [13], the sequences were OTU clustered with a similarity of 97% [14,15], and chimeras were eliminated. The RDP classifier (version 2.2) [13] was used to classify and annotate each sequence, and then compared with the Silva 16S rRNA database (v138) with an alignment threshold of 70%. Solexa was then used, which is sequencing-by-synthesis (SBS) technology that uses DNA clustering, bridge PCR, and reversible blocking methods to read the base sequences in DNA in turn. Optical equipment is used to read and record the fluorescent signal; finally, the computer converts the analysis results into sequencing information.

The data obtained in the experiment were analyzed and processed using the Origin 8.0 software, OriginLab, (OriginLab, Redmond, WA, USA) processing system, and Windows Excel and Word (2003, 2010 editions) office software, (Microsoft, OriginLab, Redmond, WA, USA.) All determinations were carried out in triplicate, and the mean values were presented. The data are reported as the averages of three separate experiments ± SD (*n* = 3).

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

#### *3.1. Water Quality Measurement Results*

As important freshwater habitats for microorganisms, urban rivers are mainly affected by the terrestrial environment and human activities. When pollutants flow into the river, basic parameters such as pH, temperature, dissolved oxygen, salinity, and nitrogen and phosphorus content in the water body change, which is of great significance for regulating the complex microbial community and further reshapes the bacterial community structure [16–18]. For the microbial community structure at the water–sediment interface found in the Chaka Salt Lake (China) study, the apparent changes in the lake water salinity were associated with a decrease in lake water [19]; studies in the Ebro River (Iberian Peninsula, Spain) revealed that electrical conductivity, temperature, and dissolved inorganic nitrogen are the main environmental factors affecting the composition of phytoplankton communities [20]. The water quality indicators of the sampled water bodies were measured; the results are shown in Table 3. The sampling areas had obvious differences, particularly the Hechuan District, which had a high ammonia nitrogen content, a high total nitrogen content, and a low chemical oxygen demand. Furthermore, by monitoring the water quality of the four sampling points, it was found that the water quality of the four points was polluted to a certain extent. The indicators of nitrogen and phosphorus exceeded the normal values of the water body; in particular, the indicators of algae and chlorophyll a exceeded the average value of the river body. Therefore, it is speculated that the four sampling points are prone to the risk of water eutrophication in summer and autumn. Changes in the river environment can affect the balance between alien bacteria introduced by surface runoff and indigenous communities, predation or competition in the food chain, and directly or indirectly affect the function of bacterial communities [21]. Therefore, understanding river bacterial community diversity and spatial distribution has profound implications for monitoring ecosystem health.

**Table 3.** Measurement results of water quality indicators.

#### *3.2. Alpha/Beat Diversity Analysis*

Alpha analysis is a composition analysis—the similarity of the species composition of each sample or group could be reflected. Determining bacterial metabolic potential through river bacterial diversity is an extremely important but under-recognized scientific issue. Changes in the bacterial structure of rivers can further lead to changes in metabolic and environmental information processing pathways [22,23]. From the upper to the lower reaches of the Danube, the relative abundance of its typical freshwater flora gradually increased due to a decrease in microbial richness and evenness [24]. The microbial diversity of the Yenisei peaked in the middle section, which may be associated with high flow and rapid turbulence [25]. As seen in Figure 2, the Shannon indices of the samples were all between 3.2 and 4.5 and the species composition richness and diversity of the four groups were relatively high. PCA analysis showed that the four groups could clearly be distinguished, indicating that the four groups had significant differences in the composition of the biological community. The R value was relatively close to 1, indicating that the difference between the groups was slightly greater than the difference within the group. There were obvious differences in the community composition structure of the four groups and the community composition in the different regions was different. The reason for

this may be that Jiangbei District has the largest industrial area, agricultural area, and population, which has a greater impact on water bodies, so the diversity of bacterial species is the highest; Beibei District and Hechuan District are rich in ecological communities, so the microorganisms in the riparian zone may be due to rainwater washed into the lake, so the microbial community is relatively rich, while the Shapingba area has the weakest community diversity because the industrial area, agricultural area, and ecological community zone are relatively less than the above-mentioned sampling points. Therefore, identifying the bacterial diversity of rivers can help to determine the impact of the metabolic potential of river bacterial communities.

**Figure 2.** Alpha Shannon index diagram/PCA analysis diagram. (H represents the Hechuan District, S represents the Shapingba District, J represents the Jiangbei District, and B represents the Beibei District).

#### *3.3. Community Composition Analysis*

#### 3.3.1. Venn Diagram Analysis

A Venn diagram permits a macroscopic reflection of the number of species in the community through an intersection, which can reflect the composition of species. As seen in Figure 3, the Hechuan District and Beibei District shared 308 species, the Hechuan District and Shapingba District shared 217 species, the Hechuan District and Jiangbei District trict shared 319 species, the Beibei District and Shapingba District shared 243 species, the Beibei District and Jiangbei District shared 232 species, and the Shapingba District and Jiangbei District share 248 species. The Beibei District had the greatest number of distinct species (109 species), and the Shapingba District had the smallest number of distinct species (9 species). This may be due to riparian influences leading to the proliferation of exotic bacteria (i.e., soil and groundwater) and indigenous flora (i.e., increased bacterial richness). For example, studies on the Mississippi River (US) [26], Thames River (UK) [27], and Amazon River [28] have found that the "heterologous input process" or "mass effect" in the riparian zone further indicates that the diffusion rate of organic matter exceeds the local bacterial extinction rate, making the exogenous input bacteria [29–31]. The diversity is much higher than the river itself.

**Figure 3.** Venn diagram analysis of community distribution.

#### 3.3.2. Community Bar Chart Analysis

Through a bar chart analysis, it could be seen that the measured water samples had a significantly high bacterial diversity. Figure 4 showed the bar map at the phylum level. The community composition structure of the Beibei District and Jiangbei District was relatively similar. The higher relative abundance in the Hechuan District were *unclassifiedf -Comamonadaceae* [32], *unclassified-f -Rhodobacteraceae* [33], and the *Hgcl\_clade* [34]. The dominant flora was the *Hgcl\_clade*, and its abundance in the three samples in this group was 12.35, 12.25, and 15.28%, respectively. The relatively high relative abundance levels in the Beibei District were *unclassified\_f\_Comamonadaceae*, the *Hgcl\_clade*, and *Acinetobacter*. The dominant flora was *Hgcl-clade*, whose abundances in the three samples were 12.98, 10.78, and 11.57%, respectively. In the Shapingba District, the relative abundances form *Limnohabitans*, *Acinetobacter*, and *Vogesella.* The dominant flora was *Limnohabitans*, whose abundances in the three samples were 13.17, 14.65, and 17.51%, respectively. The most abundant species in the Jiangbei area were *Limnohabitans*, *unclassified-f-Comamonadaceae*, and the *CL500-29-marine-group*. The dominant flora was *Limnohabitans*, whose abundances in the three samples were 11.18, 11.56, and 11.81%, respectively.

**Figure 4.** Analysis of the bacterial colony structure at the genus level of the four groups. (H represents the Hechuan District, S represents the Shapingba District, J represents the Jiangbei District, and B represents the Beibei District).

#### 3.3.3. Heatmap Analysis of Community Composition

The heatmap displayed the species composition information of the community by reflecting the number of species through a color ladder. As can be seen in Figure 5, the community composition of the three samples collected from the Jiangbei District and the Beibei District was relatively uniform; one of the three samples collected from the Hechuan District and the Shapingba District had a few differences. The community composition of Hechuan and Beibei was the most similar. *Limnohabitans* [35], *comamonadaceae*, *hgcl\_clade*, and *CL500-29-marine-group* had a higher abundance in the four samples; in the Shapingba Distric, it also contained high levels of *Acinetobacter* and *Vogesella*, and Hechuan District also contains high levels of *Aquabacterium*. Bacilli are the predominant flora in urban river water and sediments, which are associated with the gut flora of many mammals and humans, and are effective fecal surrogate indicators. The higher abundance of bacilli in the downstream may be due to the high concentrations of nitrate, ammonia, feces, etc. in medical waste liquid and domestic wastewater [36]. The results of this study showed that the concentration of pollutants gradually increased along the flow direction of the river due to the discharge of domestic and industrial wastewater, which seriously affected the abundance changes of the bacterial community in the river.

#### 3.3.4. Pie Chart

We statistically analyzed the species abundance of the samples at a certain taxonomic level and present the analysis results in the form of pie charts to present the composition of the species in different groups.

Figure 6a reflects the species composition of the Beibei District, which can clearly reflect that *unclassified\_f\_Comamonadaceae*, *Hgcl\_clade*, and the *CL500-29 marine group* were the dominant species in the Beibei District with a relatively high relative abundance that accounted for the total community of 28.8%. Figure 6b reflects the species composition of the Shapingba District, which clearly reflects that *Limnohabitans*, *Vogesella*, and *Acinetobacter* were the dominant species in the Shapingba District. The relative abundance of *Limnohabitans* was relatively high, accounting for 15.12% of the total community. Figure 6c reflects the species composition of the Jiangbei District, which clearly reflects that *Limnohabitans*, *unclassified\_f\_Comamonadaceae*, and the *CL500-29 marine group* were the dominant species in the Jiangbei District. The relative abundance of *Limnohabitans* was high, accounting for 11.52% of the total community. Figure 6d reflects the species composition of Hechuan District, which clearly reflect that *unclassified\_f\_Comamonadaceae*, the *Hgcl\_clade*, and the *CL500-29*

*marine group* were the dominant species in the Hechuan District. The relative abundance of *unclassified\_f\_Comamonadaceae* was high, accounting for 13.30% of the total community. Combining Venn diagrams, community Bar diagrams, Heatmap diagrams of community composition, and pie diagram test results can help us explore the species composition of the communities at the four sampling points, including the number of species and species abundance, by comparing samples or groups. The differences in species composition can reflect the composition information between them more clearly, and provide detailed parameters for monitoring and judging the pollution of the entire Jialing River.

**Figure 6.** (**a**) Community composition analysis map of the Beibei District, (**b**) community composition analysis of the Shapingba District, (**c**) community composition analysis of the Jiangbei District, and (**d**) community composition analysis of the Hechuan District.

#### *3.4. Analysis of Species Differences*

#### 3.4.1. Kruskal–Wallis H Test Analysis

In this study, the Kruskal–Wallis H test was used to analyze the species difference of the samples. As seen in Figure 7, each group contained the genus *Limnohabitans*, the genus *unclassified\_f\_Comamonadaceae,* and the genus of the Hgcl\_clade with a relatively high relative abundance. Compared with other groups, the species with significant differences in the Shapingba District were *Acinetobacter* [37], *Vogesella*, and *Pseudomonas*, all of which had a higher relative abundance, whereas *Cavicella* had a lower relative abundance. The species with significant differences in the Hechuan District were *Aquabacterium*, *Novosphingobium*, and *Cavicella* with a higher relative abundance, whereas *norank\_Vicinamibacteraceae* had a lower relative abundance. The species with significant differences in the Beibei District had a higher relative abundance of *norank\_Vicinamibacteraceae* and a lower relative abundance of *Vogesella.* The species with significant differences in the Jiangbei District were *Acinetobacter*, *Vogesella,* and *Aquabaterium*; all of which had relatively low relative abundances.

**Figure 7.** Kruskal–Wallis H test analysis results. (The left side of the vertical axis represents the species name and the right side is the *p*-value. The abscissa represents the relative abundance. H represents the Hechuan District, S represents the Shapingba District, J represents the Jiangbei District, and B represents the Beibei District).

#### 3.4.2. LefSe Analysis