**3.3 Comparison of slow earthquakes between simulation and observational results**

Fig. 9e shows that a large aseismic slip event occurs locally between LA and the SA belt (indicated by the ellipse) approximately one year before the megathrust earthquake. The spatiotemporal region enclosed by the ellipse in Fig. 9b shows that the large aseismic slip event triggers slow earthquakes in the SA belt with a shorter recurrence interval compared to other areas (|Strike| > 20 km). The migration distance in the SA belt corresponds to the size of the large aseismic slip region (|Strike| < 20 km). This behavior is similar to the longterm SSE (slow slip event) observed at Bungo Channel in 2003, where nearby LFT (low frequency tremor) migration had occurred either at a shorter recurrence interval or nearly continuously for several months (Obara, 2010). Therefore, the activity of the LFTs may be useful to estimate the duration and the location of local aseismic slip events, such as the long-term SSE and the preseismic slip in the deeper part of LA.

As pre- and post-seismic changes, intense LFT activity began to occur almost directly below the 2004 Parkfield earthquake about three weeks before the earthquake and has continued only apart from the hypocenter over for four years (Nadeau & Guilhem, 2009; Shelly, 2009), which means that the distance from the hypocenter of triggering earthquake may also affect the sensitivity of LFT to pre- and post-seismic slip (Shelly, 2009).

Our simulation shows that the moment release rate of slow earthquakes near the locked region of LA becomes higher about one month before the megathrust earthquake as mentioned in the section 3.2, and Figs. 9ab show that a shorter recurrence interval of slow earthquakes in the SA belt occurs several years after the megathrust earthquake. These simulation results are consistent with those results observed in Parkfield for the pre- and post-seismic stages of nearby large earthquakes.

Since SSE, pre- and post-seismic slips are all transients of interplate-slip faster than *V*pl, the observations reported by Obara (2010), Nadeau & Guilhem (2009), and Shelly (2009) would support our suggestion that the preseismic slip of megathrust earthquakes can be practically detected by monitoring slow earthquake migrations, even if the actual *dc* is less than several centimeters.


*3.4. Variance Decomposition*

In order to compare the relative contribution degree of unit standard deviation in *EC* and *NU*, we further adopt the variance decomposition. The variance decomposition results of period 5, period 10, and period 20 are displayed in Table 3 for different groups: (a) eastern, (b) central, (c) northeastern, and (d) western regions. As shown in Table 4, for provinces in the eastern, central, and western regions, the variance decomposition of *EC* is dominated by its own shock, and *NU* has a small effect. However, for provinces in the northeastern region, the variance decomposition of *EC* is dominated by the effect of *NU*, and its own shock has a small effect. For provinces in the eastern region, the variance decomposition of *NU* to its own shock is almost equal to the effect of *EC*. However, for provinces in central, northeastern, and western regions, the variance decomposition of *NU* is dominated by its own shock, and *EC* has a small effect. Combined with the variance decomposition of both variables, the effect of *EC* on *NU* is much larger than the effect of *NU* on *EC* for provinces in the eastern region, yet it is the opposite for provinces in central, northeastern, and western regions. These variance decomposition results also further confirm that energy consumption may bring about greater advantages relative to its disadvantage for provinces in the eastern region, yet it is imperative to cope with the possible cost of environmental pollution related to energy consumption for the provinces in central, northeastern, and western regions.

**Table 4.** Variance decomposition.

So as to visually display the relative contribution degree of unit standard deviation in *EC* and *NU*, we further draw the variance decomposition results of *EC* and *NU* in the form of a coordinate axis. Figure 5 shows the variance decomposition of *EC* to *NU* for different groups: variance decomposition of *EC* to *NU* (Figure 5a) and variance decomposition of *NU* to *EC* (Figure 5b)*,* respectively. As shown in Figure 5a, *NU* has a remarkable influence on *EC* for the provinces in the northeastern region, and the contribution rate is the largest among the four groups. Additionally, *NU* has an increasing contribution rate to unit standard deviation in *EC* for all groups. As explained above, these results imply that it is imperative to cope with the possible cost of environmental pollution related to energy consumption. As shown in Figure 5b, *EC* has decreasing contribution rate to unit standard deviation in *NU* for provinces in the central, northeastern, and western regions, yet *EC* has an increasing contribution rate to unit standard deviation in *NU* for provinces in the eastern region, which became the largest after the tenth forecast period. These results also imply that environmental pollution related to energy consumption may gradually increase as the pace of new-type urbanization accelerates for provinces in the central, northeastern, and western regions, which are consistent with the estimation results of the PVAR model.

**Figure 5.** Variance decomposition of *EC* and *NU*.

#### **4. In-Depth Analysis of Energy-Related Data**

Finally, the pollution-related calculations are included for an in-depth analysis. We collect supplementary data, such as carbon emissions per real CNY 10,000, energy used per real CNY 10,000 (energy consumption intensity), the ratio of output value in the tertiary industry to real GDP, and the average per capita real GDP. All the relevant data are from each year's China Statistical Yearbooks, China Energy Statistical Yearbooks, each provincial statistical yearbook, and so on. As original data for the carbon emissions can be obtained from 2003 to 2019, the sample data span from 2003 to 2019. Figure 6 shows these calculations for different groups: carbon emissions (Figure 6a), energy consumption intensity (Figure 6b), industrial structure (Figure 6c), and per capita real GDP (Figure 6d), respectively.

**Figure 6.** Some energy-related data.

As shown in Figure 6a, the largest carbon emitter appears in the western region, followed by the northeastern and central regions, and the least carbon emitter happens to be in the eastern region. As the pace of the new-type urbanization accelerates, provinces in the eastern region may make great strides to reduce the pollution of carbon emissions related to energy consumption. These results undoubtedly denote that there is an EKC relation. As shown in Figure 6b, the western region has the highest energy consumption intensity, followed by the northeastern and central regions, and finally the eastern region. As energy consumption intensity represents the efficiency of energy use, the characteristics of this variable in the four different groups denote that the energy is the most efficiently used in the eastern regions, followed by central and northeastern regions, and finally the western region. As the pace of new-type urbanization accelerates, provinces in the eastern region may successfully change their energy mix, and a decrease in carbon emissions is also can be anticipated. For example, pollution-free electrical energy can be massively put to use. To sum up, these two variables both confirm the causal relationship that new-type urbanization leads to energy consumption negatively for provinces in the eastern region, and new-type urbanization leads to energy consumption positively for provinces in the central, northeastern, and western regions over time, which seems to be in agreement with the estimation results of the PVAR model.

Figure 6c,d further provides the economic reality for the EKC relation and the causal relationship between the two. As shown in Figure 6c, the ratio of output value in the tertiary industry to real GDP in the eastern region is the highest among the four groups, followed by the northeastern region, the western region, and the central region. As shown in Figure 6d, per capita real GDP in the four groups is the same as Figure 6c. As the economy develops in the eastern region, there may be a trend toward the tertiary industry, which produces low-pollution products. Therefore, new-type urbanization leads to negative energy consumption. Provinces in the central, northeastern, and western regions are

generally eager to raise the per capita real GDP level, so standards are lax concerning relevant environmental regulations. Even some high-pollution industries are encouraged in these regions, and spontaneously their ratio of output value in the tertiary industry to real GDP is behind the eastern region. Hence, the economic reality is consistent with the EKC relation and the causal relationship between the two.


**4. Discussion**

As a chemical for pretreatment to remove external DNA on the lorica, we selected commercial bleach, which can be used on samples through the proper combination of bleach concentration and exposure duration time (seconds). Commercial bleach has been used mainly to prevent or eliminate contamination in DNA analyses (Table 1), but at the same time, it affects the body tissue of zooplankton, which can lead to the disintegration of rotifer loricas, and thus the release of rotifers' gut contents [33,35]. In this study, it was found that each rotifer species showed different response times when treated with commercial bleach at the same concentration, and the duration time for lorica survival differed by its characteristics (Table 3); the shortest time was for *Asplanchna* sp., having the softest lorica and the longest time for *Trichocerca* sp. and *Brachionus* sp. having lorica that are not easily damaged [49]. In the case of the genus *Keratella*, although its lorica has been suggested as a hard cover, which can be protective against mechanical interference by daphnids and predation by invertebrate predators [36], the loss of inner contents occurred through the mouth parts and not through lorica disintegration following the commercial bleach treatment. Therefore, based on the response time observed for *Asplanchna* sp. showing the shortest response time for disintegration after treatment with 2.5% diluted commercial bleach, we set the pretreatment time to 210 s, as this is an appropriate standard pretreatment method universally applicable to all rotifer taxa.

After applying this pretreatment method to raw water samples from the reservoir, it was confirmed that DNA fragments of rotifer food sources detected in non-treated samples were completely removed; Chlorophyceae, Diatomea, Cyanobacteria, and Ciliophora, except for bacteria (Figure 3, N, T1, T2). Further sequencing analyses indicated that the bacteria detected were mainly *Bacillus* sp. which is known to be tolerant and survive various removal treatments such as disinfection [50] (Table 4). *Bacillus* sp., gram-positive bacterium, has commonly been found in soil and other environments. It has been reported that *Bacillus* plays important roles in the lysis of bloom-forming blue-green alga and the control of their biomass in aquatic ecosystems [51,52]. Therefore, when we applied the pretreatment determined from this study to DNA analysis of gut content of rotifers, we cannot distinguish the source origins of bacteria detected in rotifer species, whether they came from contamination, water samples, or rotifer gut content, like the results of *Asplanchna* (Figure. 4D, A1~3). Since bacteria is one possible main food source for rotifers [53], a suitable pretreatment method for eliminating extracellular bacterial DNA should additionally be developed.

In the results based on DNA analysis, we used each group-specific primer for detection of targeted groups to confirm their presence/absence. As far as we know, there is no information about HNF-specific primer [54], so we applied instead a universal primer set for eukaryotes (Euk-A and Euk-B) which has been used to detect HNF (Table 2). We, therefore, carried out an additional experiment to define the applicability of the HNF primers set to rotifers. As a result, this primer set amplified all possible rotifer species from our study site except for *Asplanchna* sp. (Figure A1). These results provide a proper explanation for why the primer set did not work for all our samples (Figure 4F). The usage

of the primer sets that act specifically for each targeted biological community can help in improving detection accuracy for a targeted species group. However, there remain some limitations in verifying the e ffectiveness of a determined pretreatment on biological communities where specific primers have not yet been developed, such as HNF. In spite of these limitations, the results of the applied pretreatment method to *Asplanchna* sp. showed that specific food sources were detected in the gut content. Choanoflagellate, HNF species, has the habitat selection characteristic of being attached to phytoplankton species, and consequently it is expected that rotifers can eat Choanoflagellate indirectly in the process of eating phytoplankton, or select it as their food source directly [55]. In the case of *Meira nashicola*, which is a kind of yeast-like fungi species, although whether or not *M. nashicola* exists in aquatic ecosystem needs further research, it is considered a valid result of Asplanchna gut contents because parasitic fungus on phytoplankton, such as cyanobacteria, are known to feed on rotifers as alternative food sources [56]. So, when the limitations related to the detection of bacteria and HNF will be resolved, rotifer gut contents can be analyzed by pretreating with alcohol and commercial bleach as we recommend in the present study. Our study used traditional primer sets information; however, Adl et al. (2019) [54] recently revised the classification and nomenclature of Eurkaryotes and recommend some primer sets (rbcL, 18S V4) for Diatomea and Ciliophora. Therefore, we should apply these primer sets according to this new system for further study.

The main goal of the present study was to develop a pretreatment process that eliminated extracellular DNA fragments adhering to the Rotifera lorica and employing DNA barcoding, in order to accurately identify rotifer gut contents, thereby providing a better understanding of rotifer feeding behavior. We devised an experimental design for rotifer gut content analysis on the basis of DNA technology (DNA barcoding) while hypothesizing that feeding behavior (food selectivity) of rotifers with species-specific masticatory apparatus, e.g., the trophi, is dependent on the trophi characteristics. In this process of developing an experimental design, a pretreatment process for removing extracellular DNA as well as the cells attached to the rotifer lorica is essential in isolating accurately the DNA of the food sources remaining within each rotifer gut. Therefore, we selected appropriate chemicals for pretreatment and tried to establish the proper treatment bleach concentration (%) and duration time (seconds) by observing the response time for di fferent types of lorica firstly fixed by 60% diluted ethanol and secondly treated with 2.5% diluted commercial bleach for 210 s. The final pretreatment process was tested on a water sample and a rotifer species (*Asplanchna* sp.) to verify its e ffectiveness. We conclude that the pretreatment process for rotifer worked e ffectively in removing extracellular DNA while enabling identification of selected food source taxa of rotifers using DNA barcoding. In this study, single PCR products from group-specific primers and the general eukaryotic primers for HNF were sequenced by the cloning and Sanger method. In forthcoming studies, the taxonomic diversity of the gut content may be analyzed using next-generation sequencing (NGS) while applying improved methods for the decontamination and selection of primers in the controlled experimental environments. The DNA analysis process of rotifer gut contents, especially the pretreatment process, can allow various approaches for DNA analyses for microinvertebrates whose feeding behavior is not su fficiently understood.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-3417/10/3/1064/s1. DNA sequences: Supplementary\_Raw Sequences.FASTA (text format)

**Author Contributions:** conceptualization, H.-J.O., H.-G.J., and K.-H.C.; methodology, H.J., J.-S.G., K.-H.C. and H.J.; formal analysis, H.J., G.-J.J., S.-J.H. and H.J.; investigation, H.-J.O. and K.-H.C.; writing—original draft preparation, H.-J.O. and K.-H.C.; writing—review and editing, P.H.K., K.-H.C. and H.J.; visualization, H.-J.O., H.-G.J., and K.-H.C.; supervision, K.-H.C. and H.J.; project administration, I.-S.K.; funding acquisition, I.-S.K. and H.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Research Foundation of Korea, grant number NRF-2018R1A6A1A03024314.

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


**Author details**

Viktoras Rudzianskas\* and Rita Kupcinskaite‐Noreikiene

\*Address all correspondence to<EMAIL_ADDRESS>
Institute of Oncology, Medical Academy, Lithuanian University of Health Sciences, Kaunas, Lithuania