The study was conducted in the Sanjiang Plain, located in northeastern Heilongjiang Province, China. This region represents one of China’s largest concentrated distribution areas of freshwater wetlands, with a total area of 10.89×10⁴ km² (Wang et al., 2013). Our research focused specifically on the Honghe National Nature Reserve (133°34′38″E~133°46′29″E, 47°42′18″N~47°52′07″N) and its surrounding areas, covering approximately 21,836 hectares of pristine freshwater swamp wetland (Figure 1). The typical soil conditions in this area are mainly marsh soil and meadow soil. These soils generally have a high organic matter content and water-holding capacity, but their physical and chemical properties are significantly affected by hydrological changes and human activities (Li et al., 2018).

The range of human disturbance in the study area covers from the strictly protected core wetlands to the edge areas with dense human activities, forming a complete gradient of disturbance intensity (Li et al., 2019). Specific types of disturbances include light to moderate grazing, medium to high-intensity agricultural drainage and reclamation, as well as habitat edge effects caused by road construction. This clear gradient provides an ideal natural experimental field for systematically examining the response mechanism of plant functional traits to differentiated anthropogenic stress (Liu, 2012).

The main plant species in the area are herbaceous plants adapted to the wetland environment. The dominant species include Sedge and false sedge of the Cyperaceae family, as well as reed and purple grass of the Poaceae family (Lu and Zhang, 2014). This study systematically selected 14 dominant herbaceous plants as the analysis objects, which belong to 8 families and 12 genera, representing different ecological adaptations from hygrophytes to mesophytes, and can well reflect the functional and structural characteristics of the wetland plant community in this area.

In order to accurately test whether there is spatial autocorrelation in the sample plot, we use the classical indicator of Moran’s index. The value range of Moran’s index is between -1 and 1, which can intuitively reflect the nature and degree of spatial autocorrelation (Luke McCormack et al., 2012). Constructing an appropriate spatial weight matrix is a crucial step in the calculation process. We use distance criteria to determine the spatial weight matrix based on the actual geographical location information of the sample site (Ma et al., 2024). Set a reasonable distance threshold. When the distance between two plots is less than this threshold, their corresponding elements in the spatial weight matrix are denoted as 1, indicating that these two plots are spatially correlated; Otherwise, it is recorded as 0.

By calculating the Moran’s index based on plant functional trait data and relevant soil factor data from 24 plots, and conducting significance tests, we found significant spatial autocorrelation (Mou et al., 2018). This result implies that the plant functional traits or soil factors between adjacent plots are not independent of each other, but rather have a certain degree of similarity or correlation.

Given the significant spatial autocorrelation, it must be fully considered in the subsequent construction of statistical models. We choose to incorporate spatial coordinates into the statistical model and use a spatial lag model for analysis. This model introduces a spatial lag term for the dependent variable, fully considering the influence of adjacent plots on the observed values of the current plot. For example, when analyzing the relationship between plant leaf and root traits and soil factors under different artificial disturbance gradients, spatial lag models can more accurately analyze the actual effects of each factor.

In this study, to ensure the scientific validity of the analysis results, we conducted comprehensive and detailed normality and homogeneity of variance tests on all measurement data (Nosetto et al., 2024).

In terms of normality testing, we used the Shapiro Wilk Test. This testing method has high testing efficiency for small sample data and can accurately determine whether the data conforms to the assumption of normal distribution. We strictly follow the testing steps and input the data into professional statistical software to obtain the Shapiro Wilk test statistic and corresponding P-value for each dataset. If the P-value is greater than the pre-set significance level (usually 0.05), the data is considered to follow a normal distribution; On the contrary, it is determined that the data does not follow a normal distribution.

We used Levene’s test for homogeneity of variance. This testing method has relatively loose requirements for the distribution of data and can effectively evaluate the homogeneity of variance between different groups of data. By calculating the Levin test statistic and corresponding P-value, when the P-value is greater than the significance level, it indicates that the variance of each group of data is homogeneous; If the P-value is less than the significance level, it indicates uneven variance (O’Neill et al., 2013).

We did not directly discard data that did not meet the assumptions of normal distribution or homogeneity of variance after testing. Instead, we made appropriate transformation attempts, such as logarithmic transformation, square root transformation, etc., to make them meet the prerequisites for subsequent analysis as much as possible (Qu et al., 2018). If the conversion still cannot meet the requirements, we will use the Kruskal Wallis test in non parametric testing methods for subsequent analysis to ensure that scientific and reasonable research conclusions can be drawn under various data characteristics.

The field investigation was conducted in August 2023, and a total of 24 sample plots were set up along the human-induced disturbance gradient in the Honghe National Nature Reserve. The intensity of interference is classified into three levels - mild, moderate and severe - based on a comprehensive assessment of three quantifiable dimensions: road distance, grazing evidence and farmland proximity (Shi, 2021). Specifically, the severely disturbed plots are adjacent to roads or the edges of farmlands, with a distance from the roads not exceeding 100 meters, and are accompanied by obvious evidence of grazing activities, such as the amount of livestock manure per unit area exceeding 5 per 100 square meters, or significant reduction in vegetation height and obvious trampling marks on the soil surface. Moderately disturbed plots are 100 to 500 meters away from roads and are indirectly affected by surrounding farmlands or show signs of mild grazing. The amount of feces in these plots is generally 1 to 5 per 100 square meters. The mild disturbance sample is located at the core of the buffer zone, far from roads and farmlands. The distance from the road is more than 500 meters, and there are almost no traces of grazing activities. The amount of feces does not exceed 1 per 100 square meters or is completely absent (Song et al., 2023).

Set up a 10-meter by 10-meter large plot at the center of each plot, and arrange three 1-meter by 1-meter small plots along the diagonal for detailed vegetation investigation (Tang et al., 2024). Collect complete plants of dominant herbaceous plants such as purple-flowered wild grass, Carex, and reed. Meanwhile, soil samples were collected using a soil drill, sealed for storage and transported back to the laboratory for analysis. The coordinates of the sample plots were recorded using handheld GPS, and the soil temperature was measured in situ using a ground temperature meter.

The research area of this study has been scientifically planned and reasonably divided into 24 disturbance gradient plots. The division of these 24 plots is not arbitrary, but takes into account various environmental factors and differences in disturbance levels, aiming to comprehensively cover the ecological conditions under different levels of disturbance and provide a rich and representative sample basis for subsequent research (Thapa et al., 2022).

Within each carefully divided plot, we further set up three sampling points. The selection of these sampling points follows the principle of random and uniform distribution, striving to accurately reflect the overall ecological characteristics of the plot and avoid data distortion caused by the concentration or deviation of sampling points.

For each sampling point, we divided it into four quadrants and conducted measurements of plant functional characteristics for each quadrant. This detailed division method enables us to gain a deeper understanding of the various functional characteristics of plants from multiple perspectives and levels, and obtain more precise and accurate data (Wang et al., 2009).

Through this hierarchical sampling design, we ultimately calculate the total number of samples. The specific calculation method is: multiply 24 plots by 3 sampling points of each plot, and then multiply by 4 quadrants of each sampling point, that is, 24 × 3 × 4 = 288 samples.

Specific Leaf Area (SLA), as a key indicator for measuring the functional characteristics of plant leaves, reflects the area size corresponding to the unit dry weight of leaves (Wang et al., 2019). It can comprehensively reflect the characteristics of plants in resource acquisition, utilization strategies, and growth adaptation. The calculation method strictly follows the standard protocol in the Functional Ecology Handbook, and the specific steps are as follows:

Firstly, collect the leaves. Randomly select several healthy, mature, and representative leaves within each selected sample quadrant. To ensure the accuracy and representativeness of the data, the sampling should cover leaves from different growth positions and lighting conditions.After the collection is completed, immediately use a precise area measurement tool, such as a leaf area meter, to measure the total area of the collected leaves (in square centimeters). The measurement process should ensure that the blades are placed flat to avoid measurement errors caused by folding, curling, and other factors.Subsequently, place the measured area of the blades into an oven and dry them to a constant weight at a suitable temperature (usually 70-80 °C). The drying time depends on the thickness and moisture content of the leaves, and usually takes several hours to tens of hours. After drying, use a precision balance to weigh the dry weight of the blades in grams.Finally, according to the definition of specific leaf area, the specific leaf area is calculated by dividing the total leaf area by the dry weight of the leaf. The formula is: SLA (square centimeters/gram)=total leaf area (square centimeters) ÷ dry weight of the leaf (grams).

In each quadrat, 3–5 intact, healthy individuals of dominant species were collected, ensuring the preservation of the root-shoot connection. Samples were immediately refrigerated (<5 °C) and transported to the laboratory for processing. Leaves and roots were scanned and analyzed using ImageJ and WinRHIZO Pro 2012b software, respectively (Wang et al., 2024).

The following functional traits were measured: Leaf Area (LA), saturated leaf fresh weight, leaf dry weight (after 48 h at 65 °C), Specific Leaf Area (SLA, leaf area per unit leaf dry mass), and Leaf Dry Matter Content (LDMC, leaf dry mass per unit saturated fresh mass). Leaf total carbon (LCC), nitrogen (LNC), phosphorus (LPC), and potassium (LKC) contents were determined using a total organic carbon analyzer and an elemental analyzer after ball milling.

Root traits included root length, volume, surface area, root dry weight (after 48 h at 65 °C), Specific Root Length (SRL, root length per unit root dry mass), Root Tissue Density (RTD, root dry mass per unit root volume), and Specific Root Area (SRA, root surface area per unit root dry mass). Root total carbon (RCC), nitrogen (RNC), phosphorus (RPC), and potassium (RKC) contents were analyzed using the same methods as for leaves (Wu et al., 2010).

Soil moisture (SW, soil moisture,%), soil temperature (ST, soil temperature, °C), soil electrical conductivity (EC, soil conductivity, μS/cm), and soil salinity (SS, soil salinity, mg/L) were measured by using soil multi-parameter rapid detector. After natural air-drying, the soil samples were dried in an oven at 105 °C for 24 hours, ground and homogenized by a ball mill, and then passed through a 100-mesh sieve to determine the nutrient content. The total organic carbon content (SC, soil carbon concentration, mg/g) in soil samples was measured by Vario TOC, Elementar, Germany. The contents of total nitrogen (SN, Soil nitrogen concentration, mg/g) and total phosphorus (SP, Soil phosphorus concentration, mg/g) in soil samples were determined by an element analyzer (Vario EL III, Elementar, Germany).

All data were analyzed and charted by Origin 2021, SPSS 25.0, Canoco5, and R language 4.1. 3. Each trait variable was averaged at the quadrat level, and the variation degree of morphological traits and stoichiometric traits was calculated by using the coefficient of variation (CV). Before the analysis, the logarithmic transformation based on 10 was carried out to make the data satisfy the normalized normal distribution. Pearson correlation analysis was used to analyze the correlation between functional traits and environmental factors in different parts of plants. One-way ANOVA was used to test the difference in functional traits of leaves and roots under different human disturbance environments. Redundancy analysis (RDA) was used to explore the relationship between plant leaf and root traits and soil factors under different human disturbance types, to analyze the response mechanism of plants to different human disturbance environments.

RDA analysis identified varying environmental drivers of plant traits across disturbance intensities (Figure 7). In lightly disturbed wetlands (Figure 7a), leaf traits were primarily influenced by EC and SS (19.9% each), showing strong negative correlations with LCC and LNPR. Root traits (Figure 7b) were mainly affected by SW (12.4%) and SCNR (11.1%), with SW positively correlated with RNC and negatively with RCNR.

In moderately disturbed wetlands (Figure 8a), SP emerged as the dominant factor for leaf traits (22.5%), exhibiting strong positive correlation with LPC and negative correlations with LCPR and LCC. For root traits (Figure 8b), SC was the primary driver (22.9%), showing positive correlation with RTD and negative correlations with SRA and RPC.

Under severe disturbance (Figure 9a), ST (36.6%) and SP (22.6%) were the main factors affecting leaf traits. ST positively correlated with LNC and negatively with LCNR, while SP positively correlated with LKC and negatively with LDMC. For root traits (Figure 9b), SN (24.1%) and SP (23.2%) were dominant, with SN negatively correlated with SRA and SRL, and SP positively correlated with RPC and RKC but negatively with RCPR.

Our study demonstrates systematic variations in plant functional traits along the human disturbance gradient in the Sanjiang Plain wetlands. The observed patterns reflect distinct plant adaptive strategies to changing environmental conditions.

The correlation analyses reveal how soil factors shape trait expression (Figure 6). The positive relationships between SW and acquisitive root traits (SRL, SRA) demonstrate how water availability drives root economic strategies. Similarly, the negative correlation between soil nutrients and SRL supports the resource conservation theory—in nutrient-rich soils, plants reduce investment in nutrient-acquisition structures.

The RDA results further elucidate the changing importance of environmental drivers across disturbance levels (Figures 7-9). In mildly disturbed wetlands, where soil factors explain less trait variation (58.65-65.09%), plant traits appear influenced by biotic interactions and historical contingencies. However, as disturbance intensifies, soil factors become dominant drivers, explaining 73.82-79.08% of trait variation in severely disturbed wetlands.

Our findings reveal a continuum of plant strategies across the disturbance gradient. In mildly disturbed wetlands, plants exhibit conservative traits with high tissue density and carbon content, representing a “slow-return” investment strategy. Moderately disturbed environments select for balanced resource allocation with efficient nutrient retention. Under severe disturbance, plants shift toward acquisitive traits with high SLA and SRL, adopting a “fast-return” strategy to capitalize on dynamic conditions.

These patterns align with the leaf and root economic spectra frameworks but extend them by demonstrating how multiple traits coordinate across organs in response to complex environmental gradients. The Sanjiang Plain wetlands thus provide a model system for understanding how human disturbance reshapes plant adaptive strategies through environmental filtering.