Variability in Regional Ecological Vulnerability: A Case Study of Sichuan Province, China

Rapid urbanization and natural hazards are posing threats to local ecological processes and ecosystem services worldwide. Using land use, socioeconomic, and natural hazards data, we conducted an assessment of the ecological vulnerability of prefectures in Sichuan Province for the years 2005, 2010, and 2015 to capture variations in its capacity to modulate in response to disturbances and to explore potential factors driving these variations. We selected five landscape metrics and two topological indicators for the proposed ecological vulnerability index (EVI), and constructed the EVI using a principal component analysis-based entropy method. A series of correlation analyses were subsequently performed to identify the factors driving variations in ecological vulnerability. The results show that: (1) for each of the study years, prefectures with high ecological vulnerability were located mainly in southern and eastern Sichuan, whereas prefectures in central and western Sichuan were of relatively low ecological vulnerability; (2) Sichuan’s ecological vulnerability increased significantly (p = 0.011) during 2005–2010; (3) anthropogenic activities were the main factors driving variations in ecological vulnerability. These findings provide a scientific basis for implementing ecological protection and restoration in Sichuan as well as guidelines for achieving integrated disaster risk reduction.


Introduction
The concept of vulnerability has been applied widely in fields such as disaster studies (Yang et al. 2015;Li et al. 2016;Armaş et al. 2017), sustainable development (Turner et al. 2003), urban growth (Hong et al. 2016), and climate science (Cinner et al. 2013). Ecological vulnerability is generally conceptualized as the potential of an ecosystem to modulate its response to external interference and stressors at a specific spatial scale (Williams and Kapustka 2000;Qiu et al. 2015;Beroya-Eitner 2016). Research into ecological vulnerability has become an important aspect of study on sustainable development (Lin and Ho 2003). Spatial assessments of ecological vulnerability importantly enable the identification of areas that are vulnerable to disturbance, providing a basis for the formulation of measures aimed at controlling the degradation of ecological environments and promoting regional eco-environmental development (Ying et al. 2007;Qiao et al. 2013;Hou et al. 2015;Hong et al. 2016).
Ecological vulnerability assessment can be conducted at different hierarchical scales (De Lange et al. 2010), including population level (De Lange et al. 2009), community/habitat level (Ocaña et al. 2019), ecosystem level (Ippolito et al. 2010), and landscape level (Ding et al. 2018). The landscape level is considered a suitable scale for regional ecological vulnerability study (Qiu et al. 2007;Jiang et al. 2008), therefore we adopted the landscape-level ecological vulnerability assessment in this study to reveal regional ecological vulnerability distribution. A comprehensive indicator framework has been widely applied to evaluate ecological vulnerability at landscape level (Song et al. 2010;Zhang et al. 2014;Kumar et al. 2015), which is composed by two major parts.
The first part is establishing an ecological vulnerability indicator (EVI) system to reflect regional ecosystem vulnerability. An EVI system often entails the landscape metrics, which quantitatively describe the spatial characteristics of landscape pattern (Mörtberg et al. 2007;Qiao 2007). Qiu et al. (2007) incorporated landscape metrics of landscape isolation, fractal dimension, and fragmentation in the EVI system, to reflect the vulnerability of the ecoenvironment; Ortega et al. (2012) used metrics related to landscape composition and configuration to assess the landscape vulnerability to wildfire in Spain; Zang et al. (2017) applied landscape metrics of the number of patches, patch density, mean patch size, and class area to identify the impact of landscape pattern changes on regional ecological vulnerability. However, a shortcoming of landscape metrics is that they primarily focus on the landscape structure, while omitting its functions (Zang et al. 2017). Graph-based topological indicators may be more effective in bridging the gap between ecological structures and their functions (Kupfer 2012), offering a more systems-oriented approach for quantifying landscape processes and functions. However, to date, topological indicators have rarely been incorporated into EVI systems (Kupfer 2012). In this work, we combine landscape pattern information and topological indicators in an EVI system. The second part is to determine the weights of indicators. Common methods applied to determine indicator weights include analytic hierarchy process (Song et al. 2010;Sahoo et al. 2019), artificial neural network (Oppio and Corsi 2017), and principle component analysis (PCA) (Uddin et al. 2019). We used a PCA-based entropy method to acquire the relatively objective weights for the EVI system. Sichuan Province, with abundant natural resources, is not only an important ecological barrier in the upper reaches of the Yangtze River, but also one of the hotspots of ecological conservation research in the world (Fang 2019). Sichuan Province is the most important industrial center in western China, and is currently undergoing a rapid urbanization. At the end of 2010, the urbanization rate of the province was 40.18%, increasing to 47.69% in 2015 (National Bureau of Statistics of China 2016). The urbanization speed of Sichuan Province was higher than the national average during the period of the 12th five-year plan (Chen 2017). Meanwhile, Sichuan is a natural hazardprone region. On 12 May 2008, the Ms 8.0 Wenchuan Earthquake struck the region, leading to a loss of 35,679 ha of forestland that comprised 7.66% of the total forest area in the province (Guo 2012;Wang, Fu et al. 2012). Rapid urbanization and frequent natural hazards and disasters threaten the local ecological environment, and ecological conservation has become an essential task.
An assessment of ecological vulnerability advances understanding of the process of change of ecological systems facing various external disturbances. Our aim was to combine landscape metrics and topological indicators derived from graph theory to establish an EVI system and apply it to evaluate the regional ecological vulnerability of Sichuan Province of 3 years: 2005, 2010, and 2015. Importantly, we sought to analyze the factors driving spatiotemporal variations in ecological vulnerability to promote the protection and restoration of local ecologies.

Study Area and Data
We first introduce the study area of Sichuan Province, then describe the data involved in the present study.

Study Area
Sichuan Province, located at 97°21 0 E to 108°31 0 E and 26°03 0 N to 34°19 0 N, is situated in the upper reaches of the Yangtze River in southwestern China. The province covers an area of over 485,000 km 2 , and comprises two geographically distinct areas (Fig. 1a). The western part of Sichuan evidences a rugged terrain and a complex geological structure comprising numerous mountain ranges. By contrast, the eastern part of the province is mostly located within the Sichuan Basin, of which the westernmost part is the Chengdu Plain, and the easternmost part features a terrain comprising folds and valleys. The dramatic differences in the terrains lead to highly variable climate (Fig. 1b) and an abundance of natural resources.

Data
We used various kinds of data to conduct the assessment of Sichuan's ecological vulnerability and to identify factors driving variations in ecological vulnerability. Table 1 provides details on the datasets. Socioeconomic data include the annual statistics of gross domestic product (GDP), the number of industrial enterprises above the designated size (IEADS), highway mileage, and the urbanization rate of each prefecture. Natural hazards data include frequency of geological hazards and the magnitude of earthquakes. All raster data were uniformly resampled to a resolution of 250 m 9 250 m for this study.

Methodology
We first identified ecological patches and then established an ecological network. We subsequently developed the EVI system, incorporating the attributes of landscape pattern and the topological attributes of ecological patches. Principal component analysis (PCA) and the entropy method were used to obtain the EVI.

Identification of Ecological Patches
An ecological patch is a relatively homogeneous area that differs from its surrounding areas (Pinto and Keitt 2009). By providing habitats and necessary resources for species, patches contribute significantly to the maintenance of biodiversity (Forman 1995). We applied the following three criteria (Kong et al. 2010;Yu et al. 2012;Lee et al. 2014) to identify ecological patches.
First, we defined ecological patches as comprising two land use/land cover types: (1) forestlands, including thick forests, open woodland, shrubbery, and immature plantations; and (2) grasslands with low, moderate, and high grass cover, respectively. Compared with other land cover types (for example, cropland and construction areas),  forestland and grassland are subjected to lower degrees of anthropogenic disturbance, and provide multiple ecological services (Myers et al. 2000;Ellison et al. 2012). Second, we set the minimum area of an ecological patch at 20 km 2 based on the spatial scale of the study area and with reference to criteria adopted in previous studies (Kong et al. 2010;Yu et al. 2012;Lee et al. 2014).
Third, no two patches shared a common point or a common edge. We deemed patches within an ecological network should be spatially isolated from each other. It is assumed that if two patches shared common points or edges, they would be able to spontaneously exchange materials, energy, and information and could therefore be merged.
Applying the above criteria, we extracted ecological patches, calculated their landscape metrics, and established an ecological network for a topological analysis.

Ecological Vulnerability Index System
We selected representative indicators to construct the index system, which is the foundation of vulnerability analysis. The EVI system comprised five landscape metrics (total core area, landscape shape index, fragmentation of number, fragmentation of shape, and aggregation index) and two topological indicators (betweenness centrality and node vulnerability), respectively, depicting the landscape pattern and topological properties of patches.

Landscape Metrics
The five landscape metrics (Table 2) were calculated for each prefecture using the Fragstats 4 program (McGarigal et al. 2012). The positive indicators used in the calculation of the EVI were fragmentation of number, fragmentation of shape, and the landscape shape index, while the negative indicators were the total core area and aggregation index.

Topological Indicators
Topological indicators derived from graph theory can be used to measure the connectivity of patches. An ecological network can facilitate the topological analysis when it is analyzed as a graph, wherein the nodes and edges respectively represent ecological patches and corridors. The ecological network of Sichuan Province was established by two steps (Yu et al. 2012). The first step entailed the establishment of the cost surface, considering the following five factors: road density, construction areas, water bodies, vegetation coverage, and NPP. The second step entailed the establishment of potential corridors, indicating species' potential ability to disperse between patches. The least-cost path algorithm was applied to search for all potential dispersal corridors between each pair of patches. We subsequently applied Prim's algorithm to extract the minimum spanning tree representing the final ecological network.
We calculated two topological indicators based on the ecological network: betweenness centrality (BC) and node vulnerability (NV), and applied them in the EVI system. They are negative indicators for the calculation of EVI.
Betweenness centrality indicates the number of times that one node lies on the shortest paths of two other nodes. The node with the highest value of betweenness exerts the greatest control over the biological flow within the network (Estrada et al. 2009). Therefore, patches with higher BC values are more resilient in the face of disturbances.
where g st denotes the number of the shortest paths between nodes s and t and n i st denotes the number of the shortest paths between nodes s and t passing through node i.
Node vulnerability indicates the loss of network efficiency after the removal of a single node from the network (Grubesic et al. 2008;Chen and Hero 2013;Hu et al. 2016). A higher NV value of a node corresponds to a greater loss in efficiency with its removal.
where E denotes the efficiency of the entire network, E i denotes the efficiency of the entire network after the removal of node i, and d ij is the cost distance between nodes i and j.

Principal Component Analysis
We focused on the EVI at the prefectural level to capture regional differences in ecological vulnerability and, further, to identify their driving factors. Each indicator in the EVI system was calculated based on ecological patches and then assigned to the corresponding prefecture. For the calculation of landscape metrics, the ecological patches were first divided between the prefectures that they covered. Next, the ecological patches within each prefecture were separately entered into the FRAGSTATS software to obtain the landscape level metrics. For topological indicators, we first established the ecological network and subsequently calculated the two topological indicators for each patch. Last, we took the average value of the topological attributes for all of the patches within a prefecture as the topological attribute value for that prefecture.
Inevitably, there were cases where one patch simultaneously belonged to more than one prefecture. In such cases, the value of the topological indicator for this patch was applied in calculations for each prefecture to which it belonged. Principal component analysis is an eigenvector-based method of reducing data dimensionality that was first proposed by Karl Pearson in 1901. It enables a multidimensional system to be simplified and condensed into a smaller set of orthogonal vectors known as principal components (Eq. 4).
where Y denotes the principal components, T denotes the coefficient matrix, X represents the raw data, and p is the number of principal components. Each principal component explains a certain proportion of the total variance of the raw dataset, with components that demonstrate greater variance containing more information than those with less variance Aksha et al. 2019). The number of retained principal components was determined according to the Kaiser criterion (eigenvalues [ 1.00) (Yang et al. 2015).

The Entropy Method
Information entropy, proposed by Shannon in 1948, is a measure of the degree of disorder of a system. Building on this concept, the entropy method offers a relatively objective approach for determining weights. When there is a significant difference in the observed values for the same indicator, the corresponding entropy will be small, indicating that this indicator provides useful information.
Consequently, greater weight should be assigned to this kind of indicator; and the reverse situation also applies. Usually, the potential correlation existing between indicators is not considered when applying the entropy method, which is a drawback of this method. To overcome this drawback, we established a model for evaluating vulnerability by combining PCA and the entropy method. In this model, the principal components obtained in the PCA are inputs for the entropy method. Equations 5 to 10 illustrate the calculation process of entropy method (Du et al. 2014;Zhang et al. 2014;Wang et al. 2018).
Assuming that there are m principal components and n prefectures forming an original matrix Y ij mÂn ð Þ , we first normalized the principal components using Eq. 5.
The specific gravity value Z ij for each element Y ij is calculated using Eq. 6: Equation 7 is used to calculate the entropy of each principal component: where e i denotes the entropy of the i th principal component. The variation degree of each principal component, d i ; is measured using Eq. 8: where d denotes the variation degree of the i th principal component. Further, the weight of each principal component w i is calculated using the following formula: As a final step, we calculated the comprehensive index for assessing the ecological vulnerability of the j th prefecture as follows: where EVI j is the EVI for the j th prefecture. Prefectures with higher EVI values demonstrated higher levels of ecological vulnerability and greater vulnerability in the face of external disturbances than those with lower EVI values.

Correlation Analyses of Variations in Ecological Vulnerability and Disturbance Indicators
We used the Spearman's rank correlation analysis to test the potential correlation between EVI changes and indicators of disturbance. The Spearman's rank correlation method hinges on the ranks of two groups of data and is not sensitive to outliers.

Selection of Disturbance Indicators
The prevailing consensus is that external disturbances caused by anthropogenic socioeconomic activities and natural hazards are likely to influence ecological vulnerability (Aswani and Lauer 2014;Xie et al. 2019). To further explore the key factors driving changes in ecological vulnerability, we selected four typical socioeconomic indicators and two natural hazard indicators.

Socioeconomic Indicators
1. The number of industrial enterprises above the designated size (IEADS): The IEADS is defined as stateowned and non-state-owned industrial companies with annual revenues above RMB 20 million yuan. In recent years, the industrial sector has played an important role in the rapid development of the Chinese economy, with IEADS emerging as the primary consumers of energy and the main providers of the municipal gross industrial product (Wang, Zeng et al. 2012). 2. Highway mileage (HM): Highway mileage reflects the developmental level of regional land transportation (Cao et al. 2016). These data include public motorable roads connecting prefectures, towns, and villages; the mileage of roads passing through urban streets, the lengths of highway bridges, the lengths of tunnels, and the widths of ferry crossings.  (2008) to assign scores to earthquakes with various magnitude. Earthquakes above Ms 7.0 usually affect a wide area. Prefectures located within severely-stricken, hard-stricken, and less-stricken areas were assigned a score of 0.6, 0.3, and 0.1, respectively. Prefectures experienced earthquakes with magnitudes of Ms 5.0-5.9 and 6.0-6.9 were assigned a score of 0.1 and 0.3, respectively. The scores of a prefecture for the periods 2005-2010, and 2010-2015 were summed to assess the impact of earthquakes during the corresponding period.

Results and Discussion
In this section, we presented the spatial distribution (Sect. 4.1) and temporal variations (Sect. 4.2) of ecological vulnerability in Sichuan over the study years, and explored the possible factors that related to the ecological vulnerability variations (Sect. 4.3). Then the key findings were discussed.

Spatial Distribution of Ecological Vulnerability
Ecological patches showed similar patterns for all three study years. In terms of their spatial distribution, the patches were relatively concentrated in the northwestern part of the province, whereas they were scattered and sparse in the southern and eastern parts for each of these years (see  Figure 2b shows the spatial distribution of EVI. The values of EVI are dimensionless, which have been normalized into 0 to 1. Suining Prefecture was omitted because no patches were identified that fitted our criteria. Prefectures with high EVI values were mainly located in southern and eastern Sichuan, whereas prefectures with relatively low EVI values were mainly located in central and western Sichuan. For all 3 years, Guangyuan, Yibin, and Liangshan evidenced the highest degrees of ecological vulnerability. Patches in these prefectures were highly fragmented, relatively scattered, and small. By contrast, prefectures with widely distributed vegetation, such as Aba, Ganzi, and Ya'an, demonstrated relatively lower degrees of vulnerability.

Temporal Variations of Ecological Vulnerability
To understand differences in EVI values for the three study years, we performed the non-parametric Kruskal-Wallis one-way analysis of variance (Kruskal-Wallis one-way ANOVA) along with multiple comparison procedures (MCPs) using the MATLAB R2018b (student edition) software.
The result of Kruskal-Wallis one-way ANOVA revealed significant differences in EVI values for the 3 years (v 2 ¼ 9; p ¼ 0:011). The values obtained for the overall ecological vulnerability in 2010 and 2015 were significantly higher than the value for 2005 (Fig. 3).
The MCPs provided more detailed information about the differences between the medians of EVI of the 3 years. The EVI value in 2010 (median = 0.500) was significantly higher than that in 2005 (median = 0.409), since p value = 0.0143 \ 0.05 = a. Thus, a significant increase in the EVI value corresponded to a significant increase in the degree of ecological degradation in Sichuan Province during this period. We also observed that the EVI variance in 2010 was greater than that in 2005, indicating that regional ecological vulnerability was more unevenly distributed and that intraregional differences in ecological vulnerability have amplified. During 2005-2010, the top three prefectures evidencing the highest EVI percentage changes were Deyang (58.1%), Aba (57.7%), and Panzhihua (48.1%) (Fig. 2a, underlined in blue; Table 3). Deyang and Aba are located in areas that were struck by the Wenchuan Earthquake (Fig. 2b). Considering previous studies (Guo 2012;Wang, Fu et al. 2012) have concluded that there was severe ecological damage after the Wenchuan Earthquake, we inferred that natural hazards may be a possible factor driving EVI change, and therefore involved natural hazards indicators in the correlation analysis.
The MCPs revealed that the difference between the EVI values in 2015 (median = 0.483) and in 2010 was not significant, although there were local variations in EVI values. The difference between the EVI values in 2015 and 2005 was marginally significant (p value = 0.051). Thirteen out of 20 prefectures showed decreased EVI values between 2010 and 2015, and there was a small overall decrease comparing the median of 2015 with the median of 2010 (Fig. 3). However, the statistical results showed that the difference was insignificant. These findings imply that the ecological system has recovered to some extent, but it was still worse than the situation in 2005, indicating an urgent need for improved ecological conservation planning.

Correlation Analysis Results
We performed the Spearman's rank correlation analysis to test the potential correlation between EVI and indicators of disturbance. The results are shown in Fig. 4.
Compared with natural hazard indicators, changes in socioeconomic indicators were evidently more correlated to the EVI change (Fig. 4), indicating that human activities may be the main factors driving variations in ecological vulnerability. The result is coherent with a few previous studies (Li and Fan 2014;Liou et al. 2017), which concluded that anthropogenic activities indeed intensified the eco-environmental vulnerability.
Of the four socioeconomic indicators, highway mileage and the number of IEADS, reflecting the intensity of human construction activities, are significantly correlated with the EVI (Fig. 4a, b). The rank correlation of the percentage change in the EVI and the percentage change in the number of IEADS was highly significant (r = 0.59, p \ 0.001) (Fig. 4a). The likely reason is that productive activities within large-sized enterprises exert intensive and sustained impacts on the ecological environment in the long run. A significant positive rank correlation (r = 0.61, p \ 0.001) between the percentage change in the EVI and the percentage change in highway mileage was also observed (Fig. 4b). The construction of highways potentially results in continuous disruption of the original habitats and migration corridors of various species, thus shifting the landscape pattern and further reducing the efficiency of the exchange of materials and energy. A study of Hawbaker et al. (2006) showed a consistent result that The severely-stricken areas and the hard-stricken areas are short for severely-stricken areas of the Wenchuan Earthquake and hard-stricken areas of the Wenchuan Earthquake, which are provided by the Ministry of Civil Affairs of China roads remove habitat, alter adjacent areas, subdivide wildlife populations, as well as interrupt and redirect ecological flows. Li et al. (2010) also concluded that road construction influenced the local ecosystem by deteriorating landscape fragmentation. This finding calls for an indepth ecological impact assessment when planning construction projects.
Gross domestic product and urbanization rate, which provide more comprehensive macroscopic descriptions of socioeconomic development, did not show significant linear correlations with the EVI (Fig. 4c, d). This might be due to the fact that GDP and urbanization rate are related to multiple factors, so their relationships with EVI cannot be revealed by a linear correlation analysis.
None of the natural hazard indicators evidenced rank correlation with EVI (Fig. 4e, f), which suggests that ecosystems are resilient to natural disturbances on a long time scale. This result was echoed by Xiong et al. (2018): by 2015, the average vegetation coverage in Wenchuan Earthquake stricken area had recovered to pre-earthquake level. Besides, we noticed that the influence of earthquakes on the EVI was even less significant than that of geological hazards. A possible explanation might be that earthquakes occurred less often than geological hazards, and the selfrepairing capacities of ecosystems enabled rapid recovery from disturbances of low frequency.
Although there was a catastrophic event,-the Wenchuan Earthquake-during the period 2005-2010, a more significant factor contributing to increased EVI value might be the rapid increase in human construction activities. The construction activities contributed to a decline in the area of ecological patches and to severe fragmentation , thereby increasing ecological vulnerability. During the period 2010-2015, construction activities were reduced, and the overall ecological vulnerability decreased. In 2010-2015, although the majority of prefectural EVI values declined, EVI of Yibin, Lesha, and Ziyang Prefectures continued to increase. Simultaneously, highway mileage of these three cities increased by 43%, 35%, and 25%, respectively, ranking the first, second, and sixth among all cities in Sichuan. Since the change rate of highway mileage has a significant positive correlation with the EVI change rate, we infer that rapidly increasing highway mileage contributed to the continued increase of EVI in Yibin, Lesha, and Ziyang.
The result of EVI distribution provided a general description of ecological vulnerability situation in Sichuan Province. However, when it comes to assessing the vulnerability to a specific disturbance-for example, vulnerability to wildfire-the current EVI system may not be suitable. The EVI system needs to incorporate more eventspecific indicators-for example, the percentage of total landscape area covered by fine-grained forest-agriculture  mixtures and the total length of agricultural-forest edges per unit area (Ortega et al. 2012) to reflect the susceptibility of the ecological system to a specific disturbance. To further explore the potential relationship between EVI change and natural hazard factors, long-term studies on ecological vulnerability using higher resolution data on the intensity, frequency, and scope of various natural hazards are required. Additionally, this study did not involve much analysis between ecological vulnerability and land use/land cover changes or climate changes, which can become a crucial topic for future work.

Conclusion
In this study, we developed an integrated index system for assessing ecological vulnerability that combines landscape metrics with the topological attributes of ecological patches. The inclusion of topological attributes in a system for evaluating ecological vulnerability enables the quantification of both the landscape pattern and process, thereby contributing to a more holistic understanding of ecological vulnerability. A PCA-based entropy method was applied in the calculation of EVI to provide reasonable weight assignment and avoid correlations between indicators. We applied this index system in Sichuan Province to assess the ecological vulnerability of each prefecture for the years 2005, 2010, and 2015. The Spearman's rank correlation analysis was subsequently performed to determine the likely factors driving variations in ecological vulnerability. We conclude that: 1. Overall, the eastern and southern parts of Sichuan were more ecologically vulnerable than central and western Sichuan.

Ecological vulnerability in Sichuan significantly
increased between 2005 and 2010 and stayed high during 2010 to 2015, indicating an urgent need for improved ecological conservation planning. 3. Indicators related to intense human construction activities contributed more to variations in EVI than natural hazard indicators and macroscopic socioeconomic indicators. Among the socioeconomic indicators, the ones that concretely represent intense human activities (for example, highway mileage and the number of enterprises) were correlated more Fig. 4 Results of the Spearman's rank correlation test. PCEVI denotes the percentage change of the ecological vulnerability index; PCNIEADS denotes the percentage change of the number of industrial enterprises above the designated size; PCHM denotes the percentage change of highway mileage; PCGDP denotes the percentage change of gross domestic product; PCUR denotes the percentage change of urbanization rate; FGH denotes the frequency of geological hazards; and EII denotes the earthquake impact index. r correlation coefficient. p significance of the correlation significantly with variations in ecological vulnerability than were macroscopic indicators (for example, urbanization rate and GDP). These findings indicate that in recent decades, intensive anthropogenic construction activities have been a key factor influencing the shift in landscape patterns and affecting ecological processes. 4. The proposed assessment framework provides a new method for ecological vulnerability analysis. The framework can be extended easily to multi-scale studies, and the indicator system can be adjusted to local contexts. adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.