The study was conducted in the Jianfengling National Nature Reserve (JFR), which is located in Ledong county in the southwestern region of Hainan Island, China (latitude and longitude ranges: 18°23′N-18°50′N and 108°36′E-109°05′E, altitude range: 0–1,412 m above sea level). The region comprises mostly of an old-growth tropical montane rainforest (163 km²) at the northern edge of tropical Asia (Jiang et al., 2002; Wang et al., 2014). It has a tropical monsoon climate that is characterized with distinct wet (May to November) and dry (December to April of the second year) seasons. The mean annual air temperature is 24.5°C, and varies between 19.4°C (January) to 27.3°C (June), whereas there has a yearly mean relative humidity is 88%. The annual rainfall in the region ranges from 1,000 to 3,600 mm, with an average precipitation of 2,449 mm, of which 80–90% occurs from May to October (Zhou et al., 2013; Liu et al., 2021). This region has an irregular topography with granite as the predominant soil parent material. The most common soil type is the montane lateritic red or yellow earth (Xu et al., 2013).

A 60-ha plot was established in the tropical montane rainforest of JFR (Figure 1). It had 1,500 basic plots (20 m × 20 m) with a width of 600 m (East–West) and a length of 1,000 m (South–North; Figure 2). The plot establishment, included topographic mapping (using electronic whole-station theodolites), corner-post setting, and an initial tree census (adapting the techniques developed by the Center for Tropical Forest Science; Condit, 1998), was started in March 2009 and completed in December 2012. The altitude of the plot ranged from 866.3 m to 1016.7 m (Xu et al., 2015).

The primary vegetation of this 60-ha plot is dominated by Lauraceae, Fagaceae, Rubiaceae and Palmaceae families, along with the presence of Hopea hainanensis from Dipterocarpaceae, which is a symbol of a typical tropical flora in Southeast Asia. Here, the mean canopy is 23.7 m height, and the mean total basal area per ha is 56 m² (Jiang and Lu, 1991).

SR was measured by LI-8100 (LI-COR, Inc., Lincoln, Nebraska United States) from 20^th October to 23^th November in 2015. In the measurement mode, the LI-8100 system monitors the changes in CO₂ concentration inside the chamber by using an infrared gas analyzer (IRGA). The soil chamber had an internal volume of 4,824 cm³ with a circular contact area to soil of 317.8 cm². Polyvinyl chloride collars were placed in the center of each basic plot (Figure 2). A collar, 19.6 cm in internal diameter and 5 cm high, was buried in the soil at a depth of 2–3 cm at each point to minimize soil disturbance at least 1 week before the measurement. Plant was not present inside the collars, but small litter and branches were left in the collar. The measurements were replicated three times at each sampling point and the data was collected from 09:00 to 14:00 (Beijing standard time), which represented SR in that day (Sha et al., 2005). This time-frame was regarded as the most suitable period of the day for such types of studies because soil moisture, temperature and SOC were relatively stationary and did not influence the measurements (Yao et al., 2011).

No rainfall occurred during 9-14 h on all days from 20^th October to 23^th November in 2015 through the experimental plot. At the same time, soil temperature was measured at a depth of 5 cm using a hand-held digital thermometer (SP-E-17 thermometer, Jinzhengmao Instruments China Inc., Beijing), while the chamber temperature at a depth of 5 cm was measured by a conventional thermometer.

The temperature sensitivity of SR was usually described by Q₁₀, which was determined by the following formula:

where b was the soil temperature response coefficient and e was the natural logarithm.

A 60-ha (600 m × 1,000 m) plot was systematically divided into grids of 40 × 40 m, which generated 416 nodes that were considered as base points. Together with each base point, two additional sampling points were located at random combinations of 2 and 5 m, 2 and 15 m, and 5 and 15 m, along a random compass bearing away from the associated base point. A total of 1,248 topsoil (0–10 cm) samples were collected. Soil samples were air dried, sieved through a 2 mm mesh to remove coarse fragments. SOC was measured with the help of H₂SO₄-K₂Cr₂O₇ oxidation method (Schumacher, 2002). SN was measured by the Micro-Kjeldahl method (Bremner, 1960), and the available nitrogen was estimated with the help of a micro-diffusion method (Mulvaney and Khan, 2001). SP were measured by using an inductively coupled plasma atomic-emission spectrometer (Thermo Jarrell Ash Co., Franklin, United States) and HNO₃-HClO₄ soil solution (McDowell and Sharpley, 2001). The level of available phosphorus was estimated colorimetrically based on measurements from a 0.03 mol l⁻¹ NH₄F and 0.025 mol l⁻¹ HCl soil solution (Allen et al., 1974).

Statistical outliers were eliminated with the help of Grubbs outlier test (the numbers of valid data points were 1,454). Normality of the datasets was assessed with the one-sample Kolmogorov–Smirnov (K-S) test before data analyzes, the datasets would be logarithmically converted if datasets was not a normal distribution. All the analyzes were done in R 3.3.0. Principal component analysis (PCA) was used to reduce the number of parameters when several parameters reflected the same underlying process.

The spatial mean SR of the entire 60-ha plot was 4.327 ± 0.0411 (standard error) μmol m⁻² s⁻¹, and the coefficient of variation (CV) of SR was 36.2%. These findings were confirmed by semivariogram analysis (Figure 3), and the best-fit semivariogram models and important parameters were provided in Table 1. The range of SR autocorrelation was 112.7 m, and the NSR was 68.1% indicated that there was a moderate spatial dependency of SR (Table 1). The contour maps, obtained by ordinary kriging for SR (Figure 4), showed a patchy distribution pattern.

On the basis of a multiple regression tree, 1,454 sampling plots were divided into six types of low altitude (altitude < 914.4, type 1), high altitude(altitude > 991.5, type 2), low convexity and low slope (altitude; 963.5–991.5, convexity < 2.5, slope < 23.6, type 3), low convexity and high slope (altitude; 963.5–991.5, convexity < 2.5, slope ≥ 23.6, type 4), low convexity (altitude; 914.4–963.5, convexity < 2.5, type 5), and high convexity (altitude; 914.4–991.5, convexity ≥ 2.5, type 6; Figure 5). Stepwise linear regression analysis results showed that the SP was a major influencer for types 1, 2, and 5, whereas the SOC and soil available nitrogen were major influencer for types 3 and 4, respectively (Table 2).

The linear mix-effect model analysis showed that SR was significantly positive correlated to SP, and negatively to the slope in all subplots (Figure 6; Table 3). The low-nitrogen and high-nitrogen plots were distinguished by k-means clustering method; the soil average nitrogen was 1.930 ± 0.00379 g kg⁻¹ (n = 831), 2.245 ± 0.00538 g kg⁻¹ (n = 623), respectively. In low-nitrogen plots, SR was significantly positive related to SP, and negatively related to the slope (Figure 7; Table 4). In high-nitrogen subplots, SR was significantly positively related to SOC (Figure 8; Table 5). Also, there was a significant index relationship between SR and soil temperature (Q₁₀ of 1.92; Figure 9).

The cumulative proportion of altitude, convexity, slope, and the four soil properties, explained by PC1 and PC2, was 57.1% of the total variance. This showed the distribution of topography and soil properties on the PC1-PC2 axis, indicated a high SR occurs in the habitats of high altitude, high SOC and high SP (Table 6; Figure 10).

Quantifying the spatial heterogeneity of SR was a challenging matter, which we accomplished with the help of estimating CV. Here, we found that the CV of tropical montane rainforests in Jianfengling was 36.2%. It was consisted with previous study results of CV for SR of tropical forests in Malaysia (33%) (Katayama et al., 2009) and in Indonesia (36%) (Ishizuka et al., 2005). These were reported over 1,600 m² (10 × 10 m) and 567 m² (8 × 10 m) grids, respectively. In another Malaysian forest, Kosugi et al. (2007) observed that the CV of SR increased with subplot size, from 29 to 39%, for plot size ranging from 25 to 2,500 m², respectively. Further, the CV for SR was 42.7% for 50 lattice positions within a 100 m × 200 m plot (Adachi et al., 2005). Moreover, the spatial variability of SR in perhumid Amazonian and Malaysian rainforests were characterized by CV ranging from 24 to 45% (Schwendenmann et al., 2003; Sotta et al., 2004), and from 26 to 62% (Kosugi et al., 2007; Katayama et al., 2009), respectively. The CV of SR for rainforests in Xishuangbanna was 42% in the rainy seasons and 38% in the dry seasons (Song et al., 2013). Overall, these results showed that the CV of SR in tropical forests generally varied from 20% to 50%.

In this study, SR of 60-ha plot had a moderate spatial dependence. Some scholars found SR had a strong spatial dependence in a Southeast Asian tropical rainforest (Kosugi et al., 2007). Investigation of a 20-ha plot, in a seasonal tropical rainforest in Xishuangbanna, also revealed a strong spatial heterogeneity (Song et al., 2013). In this study, the spatial heterogeneity of SR, at small scales (below 20 m) was caused by random factors, accounting for 68.1%, while it was caused by spatially autocorrelated factors at relatively large scales (20 to 112.7 m). There was little information about the spatial structure of SR in tropical forests, and additional studies were needed.

In our study, SR exhibited spatial autocorrelation at the larger scales, across the subplots, with semivariogram ranges of 112.7 m, which have been caused by the influence of large-scale edaphic or topographic gradients, as previously reported by La Scala et al. (2000). They have measured spatial autocorrelation in CO₂ emissions as a result of a heterotrophic activity (bare soil) at scales from 29 to 58 m (La Scala et al., 2000). Kosugi et al. (2007) determined the ranges between 4.4 and 24.7 m on a 50-by 50-m plot for a tropical rainforest. The average range of the spatial autocorrelation for SR was 40 m in rainforests in Xishuangbanna (Song et al., 2013). Furthermore, Pringle and Lark (2006) detected a range of 174 m by using 156 sampling points on a 1,024-m transect, for which incubated bare soil samples were used to measure SR. Since the scales of spatial autocorrelation of SR appear to be ecosystem-and species-specific, nested sampling designs (i.e., two or more sampling intervals) in situations such as ours were recommended (Fox et al., 2015). The scale of measurement of spatial autocorrelation (less than 20 m) from this study might be used as a guide to stratify future sampling schemes in tropical montane rainforests or similar ecosystems. Since the range values were always greater than the chosen grid size, the sampling design was inferred to be adequate for the study site. In addition, column by one column mechanical point sampling was applied in our experiment in order to eliminate the effect of temporal variables. The style was to first do odd columns, followed by even columns. Random factors for mechanical point sampling were not considered; therefore, 1,500 pairs of parameters about distance and angle were randomly generated from 1,500 sub-sample centers as the center; then 1,500 sites were filtered to measure SR for further research.

The relationships illustrated in this study advances our understanding of the ecological variables that affected spatial patterns of SR within a forest. This was necessary for enhancing the understanding of mechanistic basis of SR models (Chen et al., 2011; Martin et al., 2012). Spatial relationships were found between SR, soil environment, and topography, but models were not identical across all the plots of 60-ha of sampled area. This suggested that SR might have responded to ecological covariates at varying spatial and temporal scales, created a complex pattern of SR that was not easily modeled. Spatial variations of SR in the forests have been ascribed to differences in soil C (Soe and Buchmann, 2005; Saiz et al., 2006). Many measurements suggested that SOC could be one of the main determinants of SR, particularly at large spatial and temporal scales (Giardina et al., 2003; Davidson and Janssens, 2006). Since it involved the conversion of organic to inorganic C, SR was ultimately controlled by the supply of C (Wan et al., 2007). And significantly positive relationships between SR and SOC had been reported (Chen et al., 2010). Further, SR was reported to be strongly related to SOC in a tropical rain forest in French Guiana (Epron et al., 2006). It was important to know the factors and mechanisms that controlled the spatial pattern of SR for scaling up and predicted patterns of soil C emissions. These results suggested that several soil topography and chemical factors, such as altitude, SOC, and SP, controlled the spatial distribution of SR at this site. Physical and chemical conditions of soil were also related directly and indirectly to SR, increased the complexity of SR. The next step in the process of investigation would be to independently confirm the inferences drawn in our study by generating more datasets.

We observed that the effect of convexity on SR was complex: in the lower N range, the relationship between SR and convexity was negative (Figure 7). However, in the higher N range, the relationship was positive (Figure 8). These effects could offset each other when the data was generalized over the whole of the 60-ha plot, and there was no significant relationship observed initially between convexity and SR (Figure 6). Similarly, the effect of available phosphorus on SR was also two-fold: when N was in the lower range, the relationship between SR and available phosphorus was positive (Figure 7). However, in the higher N range, this relationship became negative (Figure 8). These effects were offset when data of the complete 60 ha-plot was not parsed and no significant relationship between available phosphorus and SR was observed (Figure 6). Tropical ecosystems that were typically developed on old and nutrient-impoverished soils were considered to be rather P limited (Vitousek et al., 2010).

Although soil temperature played a leading role in controlling the temporal patterns of SR, this factor rarely explained the spatial variation of SR (Sotta et al., 2004). Here, we found a weak correlation between the soil temperature and SR (R² < 0.3; Figure 9). Furthermore, topographically induced microclimates and variations in soil temperature and soil water content could also cause spatial heterogeneity by affecting soil’s ability to retain C, water, and nutrients (Kang et al., 2000).

Yang et al. (2018) found that the average emission rate of soil CO₂ during daytime was very close to those for the time from 9:00 to 12:00. Therefore, it was crucial to measure soil CO₂ at this time points as these were the best representatives of the daily flux means. In our study, the measurements of the 1,500 sampling points were carried out during 09:00 to 12:00, and were basically the representatives of the average SR in that day. In our study, the spatial heterogeneity was surveyed only in the dry season. Therefore, more work will be needed to conduct in wet season to fully understand the temporal variation of SR in the tropical rainforest.