Functional redundancy is an important tool for justifying and prioritizing species protection in forest ecosystem, but it is a scale-dependent. If functional redundancy really exists, functional trait composition tends to have higher predictive ability of community assembly than species composition. Thus, comparing the differences in the predictive ability of community assembly between species and functional trait compositions across spatial scale represents a useful tool to quantify how functional redundancy varies across spatial scales. Here, we used variation partitioning in combination with distance-based Moran’s eigenvector maps to compare the differences in the predictive ability of community assembly between species composition and functional trait composition across spatial scales (20, 30, 40, 50, and 100 m) in a 20-ha subtropical forest plot. We found that functional trait composition possessed higher predictive ability of niche-based abiotic filtering process than species composition within 40 m. At 50 and 100 m scales, both species and functional trait compositions had approximately equal predictive ability of dispersal limitation processes. Thus, functional redundancy can only exist within 40 m scale but not 50 and 100 m scales. As a result, priority species loss protection should be performed at 50 and 100 m scales.
Introduction
Understanding the relative contributions of niche- and neutral-based processes to the diversity and assembly of species-rich plant communities remains a fundamental goal in ecology (Chave, 2013; Chesson, 2000; Ricklefs, 2004; Zhang et al., 2018). Niche-based processes predict strong correspondence between the distance decay in species similarity and the spatial scale of environmental variation or neighborhood interactions (Ackerly, 2004; Gilbert & Lechowicz, 2004; Nekola & White, 1999), whereas neutral theories predict a distance decay in species similarity that is purely determined by the spatial signature of dispersal processes (Hubbell, 2001). Therefore, quantifying patterns of species similarities across spatial scales has been one of the key tools in assessing the relative important of niche- and neutral-based processes in community assembly (Siefert, Ravenscroft, Weiser, & Swenson, 2013).
Recent developments in variation-partitioning methods combined with distance-based Moran’s eigenvector maps provide a useful set of tools to assess the relative contributions of niche-based and neutral-based processes on variation in species similarity and has been widely employed to study the processes that structure plant communities (Legendre et al., 2009; Punchi-Manage et al., 2014). However, this approach has limitations when the metric of species similarity is based purely on differences in species relative abundance. Just as Leinster and Cobbold (2012) pointed out, using species relative abundance to reflect species similarity is based on a crude model in which, as long as species have different species relative abundance, they are assumed to have nothing in common, contrary to what every biologist knows. Ecologically, a species is a collection of individuals with phenotypic and behavioral traits that determine intraspecific variation between different individuals (Cadotte, Carscadden, & Mirotchnick, 2011; Leinster & Cobbold, 2012; McGill, Enquist, Weiher, & Westoby, 2006). This view of species as an assemblage of traits make community ecologists increasingly realize that community assembly processes may be more reflected in the “traits of individuals (i.e., trait composition)” than in the “species composition” due to functional redundancy (Cadotte et al., 2011; Mouillot, Graham, Villéger, Mason, & Bellwood, 2013). That is because, functional redundancy is based on the observation that different species have the same functional roles in communities and ecosystems so that alterations in species composition cannot affect community processes and even ecosystem functioning (Lawton & Brown, 1993; Rosenfeld, 2002). In contrast, functional trait composition can directly reflect species’ functional roles in community and as a result can have higher predictive ability of community assembly processes, when functional redundancy really exists (Cadotte et al., 2011; McGill et al., 2006). However, to date, data support is still ambiguous.
It should be noted that theory based on the classical Lotka–Volterra competition model shows that stable coexistence does not allow functional redundancy, as stable coexistence requires functional complementarity among species (Loreau, 2000). Nevertheless, spatial and temporal abiotic variability may allow functional redundancy to exist at small spatial and temporal scales but not at the larger scales due to the abiotic environment for sustaining stable coexistence (Loreau, 2004). Based on this scenario, functional trait composition should have higher predictive ability of community assembly processes at the small scales wherein functional redundancy can be allowed. In contrast, as functional redundancy cannot be supported at the large scales, species composition and functional traits composition may have equal predictive ability of community assembly processes at large scales. However, to date, data supports remain ambiguous.
Niche-based deterministic processes are implicitly based on plant functional attributes that are direct measures of plant physiological tolerance to the abiotic environment and of competitive ability that together determine plant fitness (Grime, 2006; Westoby & Wright, 2006). Therefore, niche and neutral processes should also cause distinctive spatial structure in functional traits distributions (Siefert, 2012). This spatial structure in functional traits can also in turn be partitioned using a combination of constrained ordinations and distance-based Moran’s eigenvector maps. Therefore, comparing the results derived from variation-partitioning analyses based on species and traits can improve our understanding of the contribution of various processes to community assembly relative to analyses based on species or traits alone.
In this study, in a 20-ha individual tree mapped plot in a species-rich subtropical forest in Southern China, we constructed an extensive database of 20 key functional trait and abundance measurements on 112 tree species (accounting for 95% of all individuals ≥1 cm diameter at breast height), and several topographic and soil variables that represent potentially important environmental attributes, and employed variance partitioning to evaluate the, respectively, predictive power of functional and species composition in niche-based processes (habitat filtering) and neutral-based processes (e.g., dispersal limitation) at multiple spatial scales (20, 30, 40, 50, and 100 m scales). Our recent work has found that abiotic filtering and dispersal limitation dominate community assembly of this 20-ha mega plot at relatively small (20–40 m) and large (50–100 m) scales, respectively (Zhang et al., 2018). This indicated that strong functional redundancy should exist at small (20–40 m) scales due to the strong abiotic filtering-induced trait convergence at the community level. In contrast, as dispersal limitation cannot generate either trait convergence or divergence at the community level (de Bello et al., 2012; Zhang et al., 2018), functional redundancy may not be supported at the large (50–100 m) scales. As a result, we hypothesized that functional trait composition possesses higher predictive ability of abiotic filtering than species composition at the relatively small scales (20, 30, and 40 m) due to the strong abiotic filtering-induced functional redundancy. In contrast, at larger scales (50 and 100 m) wherein functional redundancy cannot be supported by the strong dispersal limitation, both species and functional trait compositions have comparable predictive ability of community assembly.
Methods
Study Site
The study site is located in the Dinghu Mountain (DHM) Reserve (112°30′39″–112°33′41″E, 23°09′21″–23°11′30″N) in Guangdong Province, Southern China. Owing to Tibetan plateau uplift and south subtropical monsoon climate, DHM has the unique subtropical monsoon broad-leave forest with high species diversity. The Reserve comprises low mountains and hilly landscapes (total area 1,155 ha), with altitudinal range from 14 to 1,000 m, covered by subtropical forests. The site has a south subtropical monsoon climate with a mean annual temperature of 20.9°C, and mean monthly temperature of 12.6°C and 28.0°C in January and July, respectively. Average annual precipitation is 1,929 mm, with most of the precipitation occurring between April and September (Li et al., 2009). A 20-ha (400 × 500 m) permanent forest plot was established in the core area of DHM in 2005. The plot features rough terrain with a steep hillside in the southeast corner. Topography varies with ridge and valley in the plot, and the elevation ranges from 240 to 470 m (Figure 1). All free-standing individual stems with diameter at breast height ≥1 cm have been identified, labeled, and mapped.
Environmental Heterogeneity Sampling
To quantify environmental heterogeneity, we measured both topography and soil properties at different scales. The topography of the DHM plot was quantified by measuring elevation at the four corners of each cell of a 20-m grid. Elevation values at the 5-m cell size was interpolated by ordinary kriging from 20-m data, while the values for larger cell sizes (i.e., 20 , 30 , 40, 50, and 100 m) were based on averages of the 5-m cells. For each cell size, we calculated the mean elevation, slope, convexity, and aspect of each grid cell (Harms, Condit, Hubbell, & Foster, 2001). To quantify soil properties at different scales, we selected 208 quadrats (20 × 20 m2) in the 20-ha plot (distributed regularly at 30-m intervals) and collected 3 0 to 20 cm samples within each quadrat which were analyzed individually to measure the following soil variables: soil bulk density, soil carbon, total nitrogen, total phosphorus, total potassium, available nitrogen, available phosphorus, available potassium, soil water content, and pH. We then used geostatistical methods (ordinary kriging) following John et al. (2007) to obtain estimates of environmental variables at each spatial scale (20, 30, 40, 50, and 100 m).
Functional Trait Collection
In this study, we measured 20 plant functional traits including leaf area (LA; cm2), leaf lamina thickness (Thk; cm), leaf dry matter content (LDMC; mg kg−1), petiole length (Pl; m), petiole dry matter content (Pdm; mg kg−1), petiole density (Pd; g cm−2), wood density (WD; g m−3), specific leaf area (SLA; cm2 g−1), nitrogen content per leaf mass (Nmass; mg kg−1), phosphorus content per leaf mass (Pmass; mg kg−1), leaf chlorophyll concentration (Chl; g m−2), maximum net CO2 assimilation rate (Aarea; µmol s−1), photosynthetic nitrogen use efficiency (µmol mol−1 s−1), photosynthetic phosphorus use efficiency (mmol mol−1 s−1), instantaneous water use efficiency (WUEi; µmol mol−1), leaf turgor loss point (ψtlp; Mpa), sapwood-specific conductivity (ks; kg m−1 s−1 MPa−1), leaf-specific conductivity (kl; kg m−1 s−1 MPa−1), stomatal conductance per unit area (gsa; mmol m−2 s−1), and stomatal conductance per unit mass (gsmass; mmol g−1 s−1) for a total of 112 tree species in the 20-ha plot, sampling 3to 5 individual adult trees of each species (Table 1). All functional traits were determined as described in previous studies (Zhang et al., 2018). The performance of each of the 20 functional root traits was shown in Table 2. The procedures for trait measurements were given in detail as below.
Table 1.
Mean Value of the Measured 20 Functional Traits.
Table 2.
Functional Traits and Their Performance in Ecological Strategies.
Chl and leaf and petiole morphological traits (LA, Pl, Pd, Thk, SLA, LDMC, Pdm)
Chl was evaluated as the average of three points on each leaf by a portable chlorophyll meter (SPAD 502, Plus Chlorophyll Meter; Konica Minolta, Ramsey, MI, USA) based on a significant positive relationship with total chlorophyll. LA, Pl, and Pd were determined using a scanner (CanoScan LiDE 700 F) and analyzed with an image processing software (ImageJ, version 1.43 u; National Institute of Mental Health, Bethesda, MD, USA). Thk was measured twice on each side of the main vein at the widest part of each leaf (to avoid major veins) using a micrometer. Leaves were then dried at 60°C for 72 h and weighed to determine leaf and petiole dry weight. Individual leaf size and petiole density were calculated from the leaf scans using ImageJ (Rasband, 1997); SLA was calculated as leaf size per unit of dry leaf mass (g), and LDMC was calculated as fresh leaf mass (g) per unit of dry leaf mass (g). Pdm was expressed as the ratio of petiole dry mass to petiole fresh mass.
Leaf nitrogen and phosphorus content per unit mass
For LA, 20 fully expanded leaves from the top of three to five mature individuals for each species were measured with an LA meter (Li-3000A; Li-Cor, Lincoln, NE, USA). Leaves were oven-dried at 70°C for 48 h to determine dry mass. SLA (cm2 g−1) was calculated as LA per dry mass. The ovendried leaves were then ground to fine powder, leaf nitrogen content per unit mass (Nmass; mg kg−1) was determined by Kjeldahl analysis, and leaf phosphorus content per unit mass (Pmass; mg kg−1) was determined using atomic absorption spectrophotometry.
Leaf gas exchange rate
Measurements of maximum net CO2 assimilation rate (Aarea) and stomatal conductance per unit area (gsa) were processed between 9:00 and 11:00 a.m. on sunny days with a Li-6400 portable photosynthesis system (Li-6400; Li-Cor). Based on preliminary trials, photosynthetic photon flux density was set at 1,500 µmol m−2 s−1 to ensure that light-saturated photosynthetic rates were measured for all species. Ambient CO2 and air temperature were maintained at 390 µmol mol−1 and 28°C, respectively. Before data were recorded, leaves were exposed to the above conditions for about 5 min to allow photosynthetic parameters to stabilize. Three to five individuals were selected for each measurement and five to six sun leaves were selected from each individual. Then, maximum CO2 assimilation rate per unit mass (Amass, µmol g−1 s−1) was calculated as SLA × Aarea/10. Stomatal conductance per unit mass (gsmass, µmol g−1 s−1) was calculated as SLA × gsa/10. WUEi was calculated as Amass/gs. Photosynthetic nitrogen use efficiency was calculated as Amass/Nmass. Photosynthetic phosphorus use efficiency was calculated as Amass/Pmass.
WD, branch and leaf hydraulic conductivity, and sapwood density
A total of 10 healthy and leaf-bearing branches (6–8 mm in diameter) from the top of three to five mature individuals for each species were cut off in early morning, sealed in black plastic bags with moist towels, and transported to the laboratory immediately. Before measurement, all of the branch samples were recut under water, and the cut ends were trimmed with a razor blade. The branch segments used in experiments were about 20 to 25 cm long. To remove air embolisms, branch segments were perfused with a filtered (diameter: 0.2 µm) 20 mmol KCl solution at a pressure of 0.1 MPa for 20 min. Each segment was then connected to a hydraulic conductivity-measurement apparatus following the method in Sperry, Donnelly, and Tyree (1988). An elevated water reservoir supplied the same perfusion solution to the segment, with a head pressure of about 6 KPa. Water flow through the segment was allowed to equilibrate for about 10 min, after which the mass of water flux though the segment over time (in seconds) was measured. Maximum hydraulic conductivity of the segment (kh) was calculated as kh = FL/ΔP, where F is the flow rate (kg s−1), ΔP is the pressure gradient (MPa) through the segment, and L is the length of the segment (m). Sapwood-specific hydraulic conductivity (ks, kg m−1 s−1 MPa−1) is equivalent to kh divided by the mean value of sapwood cross-sectional area of both ends of the branch segment. Leaf-specific hydraulic conductivity (kl, kg m−1 s−1 MPa−1) is calculated as kh/LA.
Sapwood density (WD) was determined from the same branch segments that were used for hydraulic conductivity measurements. The volume of fresh sapwood (with bark and pith removed) was determined by the water displacement method (Poorter et al., 2010), and its dry mass was subsequently determined after oven-drying at 70°C for 72 h. Then, WD (g cm−3) was calculated as the ratio of dry mass to fresh volume.
Leaf pressure–volume relationships
Leaf-bearing branches from three to five individuals of each species were harvested and transferred to the laboratory where the basal ends of the branches were immersed in distilled water and recut (5 cm removed). The branch samples were rehydrated until leaf water potential was greater than −0.05 MPa. Leaves were first weighed to obtain the initial fresh mass and then immediately placed in a pressure chamber to determine the initial water potential. Leaf mass and water potential were measured periodically during slow desiccation in the laboratory. Finally, leaves were oven-dried for 72 h at 70°C to determine their dry mass. Leaf water potential at turgor loss point (Ψtlp) was determined with a pressure–volume relationship analysis program developed by Schulte and Hinckley (1985).
Statistical Methods
Construction of functional trait composition across spatial scales
Functional trait composition can be described as community weighted mean (CWM) value. As one of the good indicators for functional trait–environment relationships, CWM describes how species respond to the environment and tends to show high sensitivity to environmental changes (Díaz & Cabido, 2001; Díaz, Cabido, & Casanoves, 1998; Vandewalle et al., 2010). Moreover, CWM can directly capture functional similarity among species and thus can be the direct measure of functional redundancy (Cadotte et al., 2011). Hence, we computed the CWM values of the 20 traits for each spatial scale. Here, species trait values were weighted by their quadrat abundance, and the weighted values were then summed over all species in the quadrat for each scale. These values for each quadrat and each trait were tabulated with quadrats in rows and CWMs of trait values in columns.
Evaluating the relative contribution of environmental, spatial factors in determining both species and functional trait compositions at multiple spatial scales
We firstly used Moran’s eigenvector mapping Legendre and Legendre (2012) to quantify spatial structure in functional trait composition at each spatial scale. Moran’s eigenvector mapping was based on the principal coordinates of neighbor matrix (PCNM) axes (Borcard & Legendre, 2002; Borcard, Legendre, Avois-Jacquet, & Tuomisto, 2004), which could also be used to describe the spatial structure or correlation of both species and functional trait compositions (Liu, Swenson, Zhang, & Ma, 2013). All spatial variables were represented by PCNM eigenfunctions calculated by principal coordinate analysis. Then, we used the R function poly to quantify the polynomial terms for each variable measured such that both linear and nonlinear relationships between abiotic variables and both species and functional trait compositions were analyzed. Based on an initial investigation of trait–environment relationships (Zhang et al., 2015), we used the first- and second-order terms for soil bulk density and pH, and the first through third-order terms for total nitrogen, available nitrogen, total phosphorus, available phosphorus, total potassium, available potassium, soil water content, soil carbon, mean elevation, slope, and convexity. Aspect (degrees from north) was decomposed into east–west and north–south orientation using sin (aspect) and cos (aspect). We also used the method developed by Blanchet, Legendre, and Borcard (2008; “packfor” library in R; R Core Team, 2013) to forward-select which of the polynomials of all measured abiotic variables and spatial variables were significant predictors of both species and functional trait compositions across spatial scales. Finally, we used variance-partitioning method (Legendre & Legendre, 2012) to allocate both species composition and functional trait composition variations as arising from the four complementary contributions (components): (a) “purely abiotic” (proportions that only can be explained by abiotic factors), (b) “spatially structured abiotic” (spatial structure in both species and functional trait compositions induced by abiotic variables), (c) “purely spatial” (proportions that only can be explained dispersal limitation), and (d) “undetermined” (Zhang et al., 2018). All variance partitioning was done across spatial scales using the function “varpart” in R.
Results
Our variance-partitioning results demonstrated that functional trait composition can directly reveal the dominating role of abiotic filtering in determining community assembly at relatively small scales (20, 30, and 40 m), whereas species composition cannot achieve this. That is because abiotic filtering represented by the sum of purely abiotic variables and spatially structured abiotic variables explained large proportions (61%–66%) of the variations in functional trait composition (Figure 3). In contrast, less than 50% of the variations in species composition were captured by abiotic filtering (Figure 2). At relatively large scales (50 and 100 m), both species composition and functional trait composition can uncover the dominating role of dispersal limitation in community assembly. That is because dispersal limitation represented by purely spatial variables (12 to 340 PCNM eigenfunctions) explained a large proportion of the variance in both species composition (51%–55%) and functional trait composition (61%–68%; Figures 2 and 3).
Discussion
Topography and soil variables all show different spatial structure; therefore, the importance of abiotic filtering increases with spatial scale (Borcard et al., 2004; John et al., 2007; Legendre et al., 2009; Portmann, Solomon, & Hegerl, 2009). However, these findings are all based on the effects of abiotic filtering on species composition. As functional redundancy can be supported at the relatively small scales (Loreau, 2004), compared with species composition, functional trait composition may better reflect the importance of abiotic filtering in community assembly at relatively small scales. Indeed, our variance-partitioning results showed that abiotic variables explain large proportions (61%–66%) of functional trait composition at relatively small scales (20, 30, and 40 m scales), indicating the dominating roles of abiotic filtering in community assembly of subtropical forest community at relatively small scales (20, 30, and 40 m scales). However less than 50% proportions of species composition are explained by abiotic variables at relatively small scales (20, 30, and 40 m scales) and thus cannot directly reveal the dominating roles of habitat filtering. Hence, assembly rules based on functional trait composition is indeed better than species composition in revealing the dominating roles of niche-based habitat filtering in community assembly at relatively small scales (20, 30, and 40 m scales).
Species composition contains 112 species abundance variables, whereas functional trait composition only includes 20 CWMs. Hence, based on the statistic assumption, species composition should show higher explaining power of community assembly compared with functional trait composition. However, the variance-partitioning results showed that functional trait composition but not species composition can reveal the dominating role of abiotic filtering at relatively small scales (20, 30, and 40 m scales). One possible reason is that strong functional redundancy should exist at relatively small scales, thereby making alteration in species composition cannot affect community assembly processes (Loreau, 2004). As CWM can directly reflect traits for dominating species (Cadotte, 2017), trait convergence in dominating species resulting from the strong abiotic filtering at relatively small scales (20–40 m) predominated community assembly (Zhang et al., 2018), and thus make many species becomes redundant. As a result, our first hypothesis that “functional trait composition possesses higher predictive ability of abiotic filtering than species composition at the relatively small scales (20, 30, and 40 m) due to abiotic filtering-induced strong functional redundancy” was supported.
At 50 and 100 m scales, the effects of abiotic variables on both species composition and functional trait composition tapered off to 5% and null. The dominant and majority effect on both species and functional trait compositions is simply due to dispersal limitation (51%–55% and 61%–68% at 50 and 100 m scales, respectively). One possible reason is that dissimilarity at large cell sizes (100 × 100 m) is quite small, suggesting that forests are relatively homogenous at such scales (De Cáceres et al., 2012; Sreekar, Katabuchi, Nakamura, Corlett, & Slik, 2018). This may be why only dispersal limitation effected tree community assembly at large spatial scales. Indeed, our previous work have found that variances of our measured abiotic variables were high at scales from 20 to 40 m, while at 50 and 100 m scales, the variances of these three abiotic variables were rather low (Zhang et al., 2018). This indicated that forests are relatively heterogeneous at scales from 20 to 40 m but homogeneous at 50 and 100 m scales. Thus, it is not surprised to see only dispersal limitation dominate community assembly at 50 and 100 m scale. As dispersal limitation can allow for neither trait convergence nor trait divergence, our results confirm that functional redundancy cannot be supported within 40 m scale (Loreau, 2004). However, we are for the first time to reveal that except for the abiotic environments for maintaining coexistence (Loreau, 2000), dispersal limitation may also be one possible mechanism for preventing functional redundancy at relatively large scales. In addition, our second hypothesis that, at 50 and 100 m scales wherein functional redundancy cannot be supported by the strong dispersal limitation, both species and functional trait compositions have comparable predictive ability of community assembly, was also supported.
It is important to note that our unexplained proportions could be attributed to unmeasured spatially structured abiotic variables and neutral processes (Legendre et al., 2009). Although we have measured a number of important abiotic factors including topographic and soil variables, some contributions of unmeasured abiotic variables (e.g., cation exchange capacity, micronutrient availability) and neutral processes to the observed spatial structures of functional and species composition cannot be ruled out and merit further investigation (Siefert et al., 2013). Moreover, although we collect 20 functional traits that have direct influence on plant performance and fitness, by no means is our data set an exhaustive set of all important functional traits. Several traits that are important for plant reproduction, dispersal, and recruitment are not included in this study. For instance, seed mass is a key trait that may be involved in competition and dispersal ability, thus representing an important axis of plant life history differentiation (Coomes & Grubb, 2003; Turnbull, Paul-Victor, Schmid, & Purves, 2008; Zhang, Gilbert, Zhang, & Zhou, 2013). Arguably therefore, our assessments on the relative importance of dispersal effects at larger scales (50 and 100 m scales) need to be further tested with more dispersal related trait data once they become available. In addition, the explaining power of abiotic and spatial variables for species compositions is inconsistent with Legendre et al. (2009) which shows that variation of species composition explained by abiotic and spatial variables increases and decreases with cell size, respectively. This may indicate a possibility of site-specific patterns, which also merit future investigation in other forest plots.
Taken together, our results for the first time to provide the empirical proof for the scale-dependent intensity of functional redundancy. More important, our results further revealed that the differences in the predictive abilities of community assembly between species composition and functional traits composition varied with the spatial-dependent intensity of functional redundancy. Niche-based assembly rules based on functional trait composition is indeed better than species composition at relatively small scales (20–40 m), wherein strong functional redundancy can be supported. However, species composition and functional trait composition have the equal predictive ability of community assembly processes at large scales, wherein strong dispersal limitation cannot allow for functional redundancy.
Implications for Conservation
The acceleration of biodiversity loss has impaired ecosystem functioning including energy transformation and matter cycling, which in turn lead to newly generated biodiversity loss (Rosenfeld, 2002). As a result, maintaining the integrity of ecosystem function is the best way to minimize species loss (Walker, 1992; Wellnitz & Poff, 2001). Based on this scenario, functional redundancy is an important tool for justifying and prioritizing species protection (Rosenfeld, 2002). That is because strong functional redundancy tends to make some species have similar contributions to ecosystem functioning, and species loss therefore has little impact on ecosystem functioning (Lawton & Brown, 1993). Thus, priority conservation effort should be put when there is little or no redundancy (Walker, 1992). Our results have found that functional redundancy cannot exist at 50 and 100 m scales for our tropical forest megaplot, as a result, priority species loss protection should be performed at 50 and 100 m scales. As dispersal limitation is the key factor for not allowing for functional redundancy at 50 and 100 m scales, active seedling of lost species may be a good way to recover species loss at large scales, if the species loss really exists.
Acknowledgments
The authors would also like to thank Christine Verhille at the University of British Columbia for assistance with English language and grammatical editing of the manuscript.
Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was funded by the National Natural Science Foundation of China (31770469), a start-up fund from Hainan University (KYQD (ZR) 1876), the Chinese Academy of Sciences through its Hundred Talent Program and Knowledge Innovation Project (KSCX2-EW-J-28), Strategic Priority Research Program of the Chinese Academy of Sciences (XDB31030000) and the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry of China.