A method for monitoring spatial variation of soil animals in black soil based on eDNA

Abstract

The application discloses a kind of black soil animal spatial variation monitoring methods based on eDNA, it is related to soil ecological monitoring technical field.The method includes: obtaining the eDNA sequencing data of different black soil sample points, different field treatment under soil sample, constructs OTU abundance table, classification annotation table and sample metadata and carries out data cleaning;Hierarchical screening target soil animal OTU;OTU abundance is standardized to Hellinger conversion;In the same place, calculate processing response value based on treatment group and control group;From response effect intensity, direction consistency, site coverage and leave-one-out stability four dimensions, calculate comprehensive score, and screen cross-site stable indicating OTU;Based on classification information and response direction vector aggregation, construct soil animal spatial variation module, and output result.The application does not need to rely on soil physical and chemical indexes and geographic distance data, can identify indicating group with robust spatial variation characteristics, realize high-throughput monitoring.

Description

Technical Field

[0001] This invention relates to the field of soil ecological monitoring technology, and in particular to a method for monitoring spatial variation of soil animals in black soil based on eDNA. Background Technology

[0002] Soil fauna in black soil is a core component of the soil food web, directly participating in key ecological processes such as organic matter decomposition, microbial community regulation, nutrient mineralization, and soil aggregate formation. Different geographical locations and agricultural management practices (such as tillage methods and fertilization systems) lead to significant spatial heterogeneity in the composition and abundance of soil fauna communities. Therefore, accurate monitoring of the spatial variation of soil fauna communities in black soil is of great significance for revealing the impact of agricultural management practices and geographical location on soil fauna community structure and for establishing a repeatable soil fauna monitoring technology system.

[0003] Current monitoring methods mainly rely on physical separation and morphological identification, which have limitations such as low throughput, long cycle time, high dependence on experts, and difficulty in standardization. Environmental DNA (eDNA) technology provides a non-destructive, high-throughput alternative. Its routine process includes sample collection, DNA extraction, target region amplification, high-throughput sequencing, and downstream diversity analysis.

[0004] Environmental DNA (eDNA) technology enables non-destructive, high-throughput monitoring of soil fauna communities by detecting residual biological genetic material in soil samples. A typical eDNA analysis workflow includes sample collection, DNA extraction, amplification and sequencing, sequence quality control, OTU clustering, taxonomic annotation, and alpha diversity, beta diversity, and differential analysis.

[0005] However, the existing technology has the following defects and shortcomings: (1) Conventional eDNA analysis usually takes the α / β diversity of the whole community as the core result, which often obscures the key indicator species (OTU) that have low abundance but have a consistent response trend to specific treatments at multiple locations; (2) When comparing between groups, the existing methods usually directly perform overall mixed comparison (such as directly merging the treatment groups and control groups at all locations), lacking the analytical logic of first pairing at the same location and then verifying across locations, which makes the analysis results easily affected by the inherent environmental noise between locations, reducing the sensitivity and accuracy of treatment effect detection; (3) In actual sampling, the treatment settings at different locations are often not completely symmetrical, and some locations may lack specific control groups. The existing standard process lacks an effective processing mechanism for such non-equilibrium data, and forced analysis will lead to improper sample pairing; (4) The existing technology cannot output an indicator biological module that has both consistent treatment response direction and cross-location spatial stability, making it difficult to directly provide spatially intuitive biological evidence for the study of carbon cycle in black soil; (5) When there is a lack of accurate latitude and longitude coordinates, the traditional distance decay model and spatial autocorrelation analysis fail.

[0006] Therefore, how to provide a method that can effectively identify indicator groups with robust spatial variation characteristics in black soil animals without relying on soil physicochemical indicators and geographical distance data, based solely on eDNA OTU abundance and classification information, is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0007] To address the shortcomings of the existing technologies, the present invention aims to provide a method for monitoring spatial variation of soil animals in black soil based on eDNA. By constructing an analysis link of intra-site pairing comparison, cross-site directional consistency screening, and spatial response module aggregation, the method solves the technical problems in the existing technologies, such as the ease with which low-abundance indicator groups are masked, cross-site comparisons being affected by environmental noise, difficulty in processing non-equilibrium sampling data, and the lack of spatially stable indicator biological module output.

[0008] To achieve the above objectives, the present invention employs the following technical solution: the method includes: Environmental DNA sequencing data of soil samples from different black soil sampling sites and under different field treatments were obtained. An OTU abundance table, a classification annotation table, and sample metadata were constructed. The OTU abundance table and the sample metadata were then cleaned. Based on the aforementioned classification annotation table, a stratified screening strategy was used to screen target soil animal OTUs from the OTU abundance table; The abundance of the selected target soil animal OTUs was normalized using Hellinger transformation to obtain normalized abundance. For samples where both treatment and control groups exist at the same location, the treatment response values ​​of each of the target soil animal OTUs at each location are calculated based on the standardized abundance. The treatment response values ​​of each target soil animal OTU are summarized across all effective locations. A comprehensive spatial stability response score is calculated from four dimensions: response effect intensity, response direction consistency, effective point coverage, and leave-one-out cross-validation stability. Based on the comprehensive spatial stability response score, cross-site stability indicator OTUs are selected. Based on the classification information and cross-location response direction vector of the cross-location stability indicator OTU, a soil animal spatial variation module is constructed by module aggregation. Output the results of soil animal spatial variation monitoring; The treatment response value is used to characterize the relative abundance change of the treatment group relative to the control group within the same location, and the effective point is the location where the effective treatment response value can be calculated. The treatment group is the sample group that received the field management measures to be evaluated, and the control group is the sample group that received routine field management measures within the same location.

[0009] In one approach, the data cleaning of the OTU abundance table and the sample metadata includes: It is mandatory that the sample number in the sample metadata exactly matches the sample number in the OTU abundance table; For samples with only a single sample or unclear treatment labels at the same location, the sample is removed from the calculation pool of the treatment response value, and only the location effect analysis used to describe the species composition of the location is retained.

[0010] In one approach, the step of using a stratified screening strategy to screen target soil animal OTUs from the OTU abundance table includes: OTUs that retain the metazoan class of Metazoa are preserved. Within the metazoa, target taxa closely related to the function of the black soil ecosystem were identified, including at least one of the following: Enchytraeida, Lithobiomorpha, Rhabditida, Psocodea, Mesostigmata, and Philodinida. Remove OTUs that are not annotated to the effective class or order level, as well as OTUs that clearly belong to plants, fungi, or protists.

[0011] In one approach, the normalization of the abundance of the screened target soil animal OTUs based on Hellinger transformation includes: The relative abundance is calculated based on the total sequencing abundance of each sample, and a Laplace smoothing correction factor, which is a positive number less than 1, is introduced when calculating the relative abundance. The relative abundance is transformed by square root to obtain the basic value of the Hellinger distance matrix, which is used as the normalized abundance.

[0012] In one scheme, the Laplace smoothing correction factor is 0.65.

[0013] In one scheme, calculating the treatment response value of each target soil animal OTU in each location based on the standardized abundance includes: The normalized abundance mean of each target soil animal OTU in the treatment group and the control group at the same location was calculated respectively. The logarithmic response ratio of the normalized abundance mean of the treatment group to the normalized abundance mean of the control group was used as the treatment response value of the target soil animal OTU at that location. A treatment response value greater than zero indicates that the target soil animal OTU is enriched in the treatment group, which is a positive response; a treatment response value less than zero indicates that the target soil animal OTU is enriched in the control group, which is a negative response.

[0014] In one scheme, the intensity of the response effect is the median of the absolute values ​​of the processed response values ​​at each of the effective points; The response direction consistency is the ratio obtained by summing the treatment response direction signs of the target soil animal OTU across all its effective points, taking the absolute value, and then dividing it by the number of effective locations of the OTU. The value range of the response direction consistency is 0 to 1. When the response direction consistency is greater than 0.8, the corresponding target soil animal OTU is defined as having strong cross-location consistent responsiveness. The effective point coverage rate is the ratio of the number of locations from which effective processing response values ​​can be calculated to the total number of locations; The stability of the leave-one-out cross-validation method is obtained by removing one location at a time and recalculating the overall response direction based on the processed response values ​​of the remaining locations. Each removal of a location is recorded as a leave-one-out test. If the overall response direction recalculated in this leave-one-out test is consistent with the direction calculated based on all valid points, it is recorded as a stable test. The stability of the leave-one-out cross-validation method is the ratio of the number of stable tests to the total number of leave-one-out tests.

[0015] In one scheme, the comprehensive score of the spatial stability response is the product of the normalized response effect strength, the response direction consistency, the effective point coverage, and the leave-one-out cross-validation stability. The step of selecting cross-site stability indicator OTUs based on the comprehensive spatial stability response score includes: selecting the top 10% of target soil animal OTUs in the comprehensive spatial stability response score as cross-site stability indicator OTUs; and marking target soil animal OTUs that have calculated effective treatment response values ​​at only a single location but have extremely high absolute values ​​of treatment response values ​​as location-specific response OTUs, and not including them in the module aggregation.

[0016] In one scheme, the module aggregation based on the classification information of the cross-location stability indicator OTU and the cross-location response direction vector includes: The processing response direction of each of the cross-location stability indication OTUs in each location is transformed into a response feature vector; When two OTUs belong to at least the same family or order in taxonomy, and the Jaccard similarity coefficient of their response feature vectors is greater than 0.75, the two OTUs are aggregated into the same soil animal spatial variation module. After successful aggregation, the module abundance is calculated, which is the sum of the relative abundances of each OTU within the module; The aggregated modules are recursively executed to calculate the in-location processing response value and evaluate the cross-location stability response. Based on the response mode, the output results are divided into four categories: location response type, processing response type, location dependent type, and stability indication type.

[0017] Beneficial effects of this invention: (1) Only the OTU table and basic metadata generated by eDNA sequencing are needed. There is no need to measure expensive soil physicochemical indicators or perform GPS positioning of sample points. The implementation cost is low and the scope of application is wide.

[0018] (2) Through rigorous sample matching and target group screening, interference from non-target organisms and abnormally labeled samples was eliminated, significantly improving the accuracy of downstream analysis.

[0019] (3) It effectively separates the treatment effect and the location environment effect, and can extract the true impact of the treatment measures on the biological community more purely compared with the traditional direct mixing comparison method.

[0020] (4) Four evaluation indicators were introduced: response effect intensity, response direction consistency, effective point coverage, and leave-one-out cross-validation stability. In particular, the original weighted comprehensive scoring model can screen out low abundance indicator OTUs that are spatially robust and have clear ecological indication significance.

[0021] (5) Based on the aggregation strategy of dual constraints of classification status and response direction vector, scattered OTUs are transformed into functional response modules with ecological significance, making the monitoring results more macroscopic and more interpretable; the detection sensitivity of treatment effect is higher under the same amount of data, and the location pairing analysis can still be completed under non-equilibrium sampling, making the monitoring conclusions more reproducible.

[0022] (6) A treatment plan for missing paired treatment locations was clearly proposed (for location effect analysis only), avoiding statistical errors caused by forcibly constructing false controls.

[0023] (7) By introducing leave-one-out cross-validation, this method naturally has the ability to eliminate outlier interference when evaluating indicator OTUs, which improves the reliability and robustness of monitoring conclusions. Attached Figure Description

[0024] Figure 1 A flowchart of the method for monitoring spatial variation of soil animals in black soil based on eDNA provided in this embodiment of the invention; Figure 2 This is a flowchart of data acquisition and preprocessing provided in an embodiment of the present invention; Figure 3 A flowchart for stratified screening of target soil animal OTUs provided in this embodiment of the invention; Figure 4 A flowchart for calculating the local processing response value and evaluating the cross-location stability response provided in this embodiment of the invention; Figure 5A flowchart illustrating the aggregation of the soil animal spatial variation module provided in an embodiment of the present invention. Detailed Implementation

[0025] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0026] The terms and abbreviations involved in this invention are explained as follows: Environmental DNA (eDNA): Genetic material released into the environment by organisms; in this invention, it refers to the genetic material of soil animals remaining in black soil samples; OTU (Operational Taxonomic Unit): A taxonomic unit obtained by clustering based on sequence similarity in microbial and environmental DNA research; in this invention, OTU is used as the basic analytical unit for soil animal groups; Hellinger transformation: A standardized method for square root transformation of relative abundance; α-diversity: A diversity index characterizing species richness and evenness within a single sample; β-diversity: A diversity index characterizing differences in community composition among different samples; Treatment group (g) and control group (h): The sample group receiving the target field treatment and the control group used as a reference within the same location.

[0027] Example 1: A method for monitoring spatial variation of soil animals in black soil based on eDNA like Figure 1 As shown, this embodiment provides a method for monitoring spatial variation of soil animals in black soil based on eDNA. It is achieved by constructing a systematic analysis process that integrates data cleaning, target screening, effect calculation, stability assessment and module aggregation, specifically including steps S101 to S107.

[0028] Step S101: Soil animal eDNA acquisition and data preprocessing. Soil samples were collected from different black soil sampling sites and under different field treatments, and eDNA was extracted and sequenced. After obtaining the raw data, three core datasets were constructed: an OTU abundance table, denoted as X=[x ij ], where x ij The sequence abundance of the i-th OTU in the j-th sample is indicated; the classification annotation table contains the OTU number and classification information of kingdom, phylum, class, order, family, genus and species; the sample metadata includes at least the sample number, site and treatment group.

[0029] The specific procedure for extracting eDNA and sequencing is as follows: Total DNA is extracted from soil samples using a soil DNA extraction kit; the hypervariable region of metazoan 18S rRNA gene V4 is used as an amplification marker, and the primer pair is TAReuk454FWD1 (5′-CCAGCASCYGCGGTAATTCC-3′) and TAReukREV3 (5′-ACTTTCGTTCTTGATYRA-3′); after purification and library construction, the amplified products are sequenced using the Illumina paired-end sequencing platform with PE250 sequencing; the raw sequences are subjected to quality control and adapter removal using FastP software, and OTUs are generated by clustering at 97% similarity using the UPARSE algorithm; the representative sequences of the OTUs are annotated with the SILVA database using RDPclassifier (confidence threshold 0.8) to complete the classification of kingdom, phylum, class, order, family, genus and species, and only OTUs with the kingdom annotation of Metazoa are entered into the target group screening in step S102 to ensure that the results of the stratified screening step are reproducible.

[0030] like Figure 2 As shown, the data cleaning criteria for this step are: it is mandatory that the sample metadata completely matches the sample number in the OTU abundance table; for samples with only a single sample or unclear treatment labels at the same location, they are directly removed from the subsequent treatment effect calculation pool, and only the location effect analysis used to describe the species composition of the location is retained.

[0031] Samples excluded from the treatment response value calculation pool were used only for site effect analysis: the OTU composition of these samples was summarized by site, the α diversity descriptive statistics (OTU richness, Shannon diversity index) of each site were calculated, and principal coordinate analysis (PCoA) was performed based on Bray-Curtis distance to output the species composition background description of each site; this part of the results was not used for any comparison between treatment groups and control groups.

[0032] Step S102: Screening of target soil animal OTUs and data dimensionality reduction. For example... Figure 3 As shown, based on the classification annotation table, a stratified screening strategy was adopted: (1) retaining OTUs at the Metazoa (metazoa) level; (2) within Metazoa, further identifying target groups closely related to the function of the black soil ecosystem, including but not limited to Enchytraeida, Lithobiomorpha, Rhabditida, Psocodea, Mesostigmata, and Philodinida; (3) removing OTUs that were not annotated to the effective class / order level or that clearly belonged to plants, fungi, or protists, to ensure that subsequent analyses were based on high-quality soil animal data.

[0033] Step S103: OTU Abundance Standardization Based on Hellinger Transformation. To avoid statistical interference from sequencing depth differences and a large number of zero values, the selected target OTUs are standardized. First, the relative abundance is calculated based on the total sequencing abundance of each sample, where the total sequencing abundance of the j-th sample is... Calculate according to formula (1): (1) In equation (1), the summation iterates through all K target soil animal OTUs in sample j, i.e., N_j only sums the counts of target OTUs after screening in step S102, which serves as the denominator of the relative abundance in equation (2) to ensure that the sample diameter is consistent under the same sequencing depth.

[0034] To reduce counting noise in low-frequency OTUs and ensure the computability of zero-value samples in subsequent logarithmic operations, a Laplace smoothing correction factor is introduced. relative abundance Calculate according to formula (2): (2) In equation (2), K is the total number of target OTUs in the sample. Set the relative abundance to a positive number less than 1, preferably 0.65, to avoid over-smoothing. Perform a square root transformation on the relative abundance to obtain the base values ​​of the Hellinger distance matrix. As shown in equation (3): (3) The innovation of this standardization procedure lies in the fact that, compared to directly using relative abundance, Hellinger transformation can give low-frequency OTUs and high-frequency OTUs reasonable weights in the Euclidean distance space, avoiding the bias of high-abundance species dominating the analysis results.

[0035] Step S104: Calculation of treatment response value within the site. The core logic of this step is to compare within the same site and avoid direct comparison across sites. For the case where both the treatment group (g) and the control group (h) exist within the same site, the local treatment effect of each OTU is calculated: the normalized abundance mean p̄i,s,g and p̄i,s,h of the i-th OTU in the treatment group and control group at that site are calculated respectively, and the treatment response value is calculated using the log response ratio according to formula (4). : (4) In equation (4), the logarithmic base is uniformly taken as 2 throughout the text, that is, the logarithmic response ratio with base 2. In the examples, the thresholds such as 0.8 and 0.75 are all calculated with base 2; the smoothing factor in equation (4) The Laplace smoothing factor ε defined in equation (2) is the same factor and has the same value. Its continued use in equation (4) is only to prevent the mean of the group after Hellinger transformation from still having a zero value, which would result in the logarithm being undefined and does not constitute a second smoothing. The OTU is enriched in the processing group (positive response). The OTU was enriched in the control group (negative response). The advantage of this data structure is that even if there are huge differences in baseline communities between locations, the logarithmic ratio can effectively filter out background noise between locations and extract the relative abundance changes caused purely by the treatment.

[0036] Step S105: Calculation and comprehensive scoring of cross-location stability response indicators. For example... Figure 4 As shown, for each OTU, the set of processing response values ​​calculated across all valid points (N) is summarized. We assess its spatial stability from four dimensions and construct a comprehensive scoring model.

[0037] (1) Response effect intensity E: Calculated according to equation (5), using the median instead of the mean, to reduce the interference of extreme response values ​​at individual locations on the overall intensity assessment: (5) (2) Consistency of response direction C: Calculated according to formula (6): (6) The C value ranges from 0 to 1. An OTU is defined as having strong cross-location consistent responsiveness if and only if C is greater than 0.8. This threshold setting is a key feature of the invention, ensuring that the selected indicator OTUs exhibit stable bioindication function under different environmental gradients.

[0038] In equation (6), the summation iterates through Ni effective points of the OTU, that is, only the locations where an effective processing response value can be calculated are counted, and the denominator is Ni; missing locations are not included in the numerator or denominator. Response direction consistency C and effective point coverage V=Ni / N measure different attributes: the former measures the degree of consistency of response direction within effective points, and the latter measures the spatial coverage of effective points. The two are independent of each other.

[0039] (3) Effective point coverage V: Calculated according to formula (7), where Ni is the number of locations for which OTU i can calculate an effective response value, and N is the total number of locations participating in cross-location aggregation. This indicator is used to distinguish between general-purpose response OTUs and location-specific OTUs: (7) (4) Leave-one-out cross-validation stability L: One location is removed one by one, and the overall response direction is recalculated based on the response values ​​of the remaining N-1 locations; each removal of a location is recorded as a leave-one-out test. If the overall response direction recalculated in this leave-one-out test is consistent with the direction calculated based on all N locations, it is recorded as a stable test; L is the ratio of the number of stable tests to the total number of leave-one-out tests, calculated according to formula (8). This index is a key indicator for evaluating the spatial robustness of OTUs and effectively controls the risk of false positives caused by a single abnormal sample point: (8) (5) Spatial stability response comprehensive score Q: The present invention constructs a weighted scoring model as shown in equation (9): (9) In equation (9), The maximum value among all OTU response effect strengths is used to... Normalization was performed. The score assigns higher weights to directional consistency and spatial robustness because, for ecological indicator species, the repeatability of response direction is more ecologically indicative than the magnitude of the response. The selection criteria were as follows: the top 10% of OTUs in terms of Qi were selected as cross-site stable indicator OTUs; for OTUs with Ni=1 but extremely high |ri,s|, they were separately marked as site-specific response OTUs and not included in subsequent module aggregations.

[0040] Step S106: Construction of the soil animal spatial variation module based on response vectors. (e.g.) Figure 5 As shown, it specifically includes: (1) Constructing response feature vectors: converting the response direction of each OTU in each location into a vector. As shown in equation (10): (10) In equation (10), the Jaccard similarity coefficient is calculated using binary encoding: the component with a response direction of +1 in the response feature vector σi is encoded as 1, and the component with a response direction of -1 is encoded as 0. The Jaccard similarity coefficient between σi and σj is the ratio of the number of positions of the components with a value of 1 after encoding to the number of positions of at least one component with a value of 1. When all the component directions of the two vectors are exactly the same, the Jaccard similarity coefficient is 1.

[0041] Module aggregation conditions: For two OTUs to be aggregated into the same spatial variation module, they must simultaneously satisfy: (a) they belong to at least the same family or order in taxonomy; (b) their response feature vectors The Jaccard similarity coefficient is greater than 0.75. (3) Module abundance calculation: After successful aggregation, the abundance of the k-th module in the j-th sample. Calculated according to equation (11), that is, the sum of the relative abundance of each OTU within the module, where Mk is the set of OTUs contained in the k-th module: (11) (4) Recursive analysis: Repeat steps S104 and S105 for the aggregated modules to obtain the module-level processing response characteristics. Finally, the output results are divided into four categories according to the response pattern: location-responsive, processing-responsive, location-dependent, and stable-indicative.

[0042] The four types of criteria are based on two dimensions: response direction consistency C (threshold 0.8) and effective point coverage V (threshold 0.5): (1) Stable indicator type: C > 0.8 and V ≥ 0.5, the module responds in the same direction at most effective points; (2) Location-dependent type: C ≤ 0.8 and V ≥ 0.5, the module can detect the response at most locations but the response direction changes with location; (3) Processing response type: C > 0.8 and V < 0.5, the response is only seen at a few locations but the direction is consistent, mainly driven by processing measures; (4) Location response type: C ≤ 0.8 and V < 0.5, the response is concentrated at a specific location. In Example 2, "stable indicator OTU" corresponds to the stable indicator type, and "location-specific response OTU" and "location-dependent response OTU" correspond to the location response type and location-dependent type, respectively.

[0043] Step S107: Output monitoring results. Output the comprehensive analysis results from each level, including: group composition, location effect analysis results, etc. Descriptive statistics on diversity, response score ranking table at the OTU / module level, spatial response pattern heatmap, and classification result report.

[0044] Example 2: Application Case of Black Soil Sample Area In this embodiment, the input OTU abundance table contains 1866 OTUs and 89 black soil samples, with complete matching of sample metadata. Based on the screening in step S102, 197 target soil animal OTUs were successfully identified; based on the smooth Hellinger transformation in step S103, the problem of uneven sequencing depth was effectively solved.

[0045] Three typical sites (Sites A, B, and C) were selected in the black soil belt. Each site had two treatments: MCM (conservation tillage management with straw mulch and no-till) and CRS (conventional rotary tillage management) as the treatment group (g). Each treatment at each site was replicated 15 times, resulting in 90 samples. One sample was removed in step S101 due to unclear treatment labeling, leaving 89 samples for analysis. Following step S104, the response values ​​r for 197 OTUs were calculated at each site. Further, following step S105, the stability index of each OTU was calculated. For example, a certain rotifer OTU showed a consistent negative response (C=1.0) at all three sites, and its stability using leave-one-out cross-validation was L=1.0. Its final overall score Qi ranked within the top 10%, and it was labeled as a stable indicator (negative response) OTU.

[0046] Based on step S106, the algorithm identified five nematodes OTUs with highly similar response direction vectors, whose pairwise Jaccard similarity coefficients were the lowest at 0.82, satisfying the aggregation condition and forming a nematode spatial response module. The abundance of this module showed a consistent decreasing trend across the treatment groups. The final results are visually presented in the form of a heatmap, clearly indicating the differentiated response patterns of different soil animal groups between locations and treatments.

[0047] In another, larger-scale monitoring exercise, based on the treatment response results of each OTU across 31 pairable locations, the response effect strength, response direction consistency, effective point coverage, and leave-one-out cross-validation stability were calculated, and a spatially stable response score was further calculated. OTUs were then sorted from highest to lowest spatially stable response score, and the following were selected: stable indicator OTUs that repeatedly occurred at multiple locations with consistent response directions; location-specific response OTUs that occurred only at specific locations; and location-dependent response OTUs whose treatment effect direction varied with location.

[0048] In the construction of the spatial variation module, OTUs are aggregated according to their classification information and cross-site response direction vectors. OTUs with the same or similar response patterns are grouped into the same soil animal spatial variation module. For each module, the module abundance, intra-site treatment response value, and spatial stability response score are recalculated, and the module-level spatial variation results are output. These results can be displayed as a site-treatment matrix, OTU heatmap, module heatmap, or response direction matrix.

[0049] Example 3: Auxiliary Verification Results Based on existing eDNA analysis results, current data can identify multiple soil animal groups and show that different groups have different location or treatment response characteristics. Existing results show that: (1) medium and large soil animal communities... Diversity did not show significant differences in location or treatment, and the analysis results are shown in Table 1; (2) Linea order The treatment effect of diversity showed a marginally significant trend (P=0.054), and the analysis results are shown in Table 2; (3) Scolopendra order The location effect and treatment effect of diversity were not significant (P > 0.05), and the results are shown in Table 3.

[0050] In Tables 2 and 3, all F-values ​​have been recalculated using the formula F = MS factor / MS residual, where MS factor = sum of squares of factors / degrees of freedom of factors, and MS residual = sum of squares of residuals / degrees of freedom of residuals. The significance P-value is calculated using an F-distribution with 2 degrees of freedom in the numerator and 7 degrees of freedom in the denominator. The calculation results are as follows: Table 2: Plot effect F = 0.656 (P = 0.548), treatment effect F = 4.569 (P = 0.054), interaction effect F = 0.298 (P = 0.751); Table 3: Plot effect F = 0.625 (P = 0.563), treatment effect F = 0.150 (P = 0.863), interaction effect F = 0.280 (P = 0.764). All data in the tables have been corrected according to the calculation results.

[0051] Table 1 Dominant taxa of medium and large soil animals Results of bivariate analysis ; The above results indicate that even when overall α-diversity is not significant, specific soil animal groups may still exhibit location-based or treatment-based responses. (Cirthales order) The results of the two-factor analysis of diversity are shown in Table 2.

[0052] Table 2. Lineae Results of bivariate analysis ; Stone centipedes The results of the two-factor analysis of diversity are shown in Table 3.

[0053] Table 3 Lithocentipedes Results of bivariate analysis ; Therefore, the spatial stability response screening algorithm described above can be used as an auxiliary analysis method in addition to α-diversity and β-diversity analysis. It should be noted that the above results are only used to demonstrate that this method can identify location and treatment differences in soil fauna communities, and are not intended to demonstrate changes in SOC content, organic carbon stability, or carbon sink efficiency.

[0054] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for monitoring spatial variation of soil animals in black soil based on eDNA, characterized in that: The method includes: Environmental DNA sequencing data of soil samples from different black soil sampling sites and under different field treatments were obtained. An OTU abundance table, a classification annotation table, and sample metadata were constructed. The OTU abundance table and the sample metadata were then cleaned. Based on the aforementioned classification annotation table, a stratified screening strategy was used to screen target soil animal OTUs from the OTU abundance table; The abundance of the selected target soil animal OTUs was normalized using Hellinger transformation to obtain normalized abundance. For samples where both treatment and control groups exist at the same location, the treatment response values ​​of each of the target soil animal OTUs at each location are calculated based on the standardized abundance. The treatment response values ​​of each target soil animal OTU are summarized across all effective locations. A comprehensive spatial stability response score is calculated from four dimensions: response effect intensity, response direction consistency, effective point coverage, and leave-one-out cross-validation stability. Based on the comprehensive spatial stability response score, cross-site stability indicator OTUs are selected. Based on the classification information and cross-location response direction vector of the cross-location stability indicator OTU, a soil animal spatial variation module is constructed by module aggregation. Output the results of soil animal spatial variation monitoring; The treatment response value is used to characterize the relative abundance change of the treatment group relative to the control group within the same location. The effective point is the location where the effective treatment response value can be calculated. The treatment group is the sample group that received the field management measures to be evaluated, and the control group is the sample group that received routine field management measures within the same location.

2. The method according to claim 1, characterized in that: Data cleaning is performed on the OTU abundance table and the sample metadata, including: It is mandatory that the sample number in the sample metadata exactly matches the sample number in the OTU abundance table; For samples with only a single sample or unclear treatment labels at the same location, the sample is removed from the calculation pool of the treatment response value, and only the location effect analysis used to describe the species composition of the location is retained.

3. The method according to claim 1, characterized in that: The step of using a stratified screening strategy to screen target soil animal OTUs from the OTU abundance table includes: OTUs that retain the metazoan class of Metazoa are preserved. Within the metazoa, target taxa closely related to the function of the black soil ecosystem were identified, including at least one of the following: Enchytraeida, Lithobiomorpha, Rhabditida, Psocodea, Mesostigmata, and Philodinida. Remove OTUs that are not annotated to the effective class or order level, as well as OTUs that clearly belong to plants, fungi, or protists.

4. The method according to claim 1, characterized in that: The standardization process for the abundance of the selected target soil animal OTUs based on Hellinger transformation includes: The relative abundance is calculated based on the total sequencing abundance of each sample, and a Laplace smoothing correction factor, which is a positive number less than 1, is introduced when calculating the relative abundance. The relative abundance is transformed by square root to obtain the basic value of the Hellinger distance matrix, which is used as the normalized abundance.

5. The method according to claim 4, characterized in that: The Laplace smoothing correction factor is set to 0.

65.

6. The method according to claim 1, characterized in that: The calculation of the treatment response value of each target soil animal OTU in each location based on the standardized abundance includes: The normalized abundance mean of each target soil animal OTU in the treatment group and the control group at the same location was calculated respectively. The logarithmic response ratio of the normalized abundance mean of the treatment group to the normalized abundance mean of the control group was used as the treatment response value of the target soil animal OTU at that location. A treatment response value greater than zero indicates that the target soil animal OTU is enriched in the treatment group, which is a positive response; a treatment response value less than zero indicates that the target soil animal OTU is enriched in the control group, which is a negative response.

7. The method according to claim 1, characterized in that: The intensity of the response effect is the median of the absolute values ​​of the processed response values ​​at each effective point; The response direction consistency is the ratio obtained by summing the treatment response direction signs of the target soil animal OTU at all valid points, taking the absolute value, and then dividing it by the number of valid locations of the OTU. The value range of the response direction consistency is 0 to 1. When the response direction consistency is greater than 0.8, the corresponding target soil animal OTU is defined as having strong cross-site consistent responsiveness. The effective point coverage rate is the ratio of the number of locations from which effective processing response values ​​can be calculated to the total number of locations; The stability of the leave-one-out cross-validation method is obtained by removing one location at a time and recalculating the overall response direction based on the processed response values ​​of the remaining locations. Each removal of a location is recorded as a leave-one-out test. If the overall response direction recalculated in this leave-one-out test is consistent with the direction calculated based on all valid points, it is recorded as a stable test. The stability of the leave-one-out cross-validation method is the ratio of the number of stable tests to the total number of leave-one-out tests.

8. The method according to claim 7, characterized in that: The spatial stability response comprehensive score is the product of the normalized response effect strength, the response direction consistency, the effective point coverage, and the leave-one-out cross-validation stability. The step of selecting cross-site stability indicator OTUs based on the comprehensive spatial stability response score includes: selecting the top 10% of target soil animal OTUs in the comprehensive spatial stability response score as the cross-site stability indicator OTUs; For target soil animal OTUs that have a valid treatment response value calculated only at a single location but have an extremely high absolute value of the treatment response value, they are marked as location-specific response OTUs and are not included in the module aggregation.

9. The method according to claim 1, characterized in that: The module aggregation based on the classification information and cross-location response direction vector of the cross-location stability indicator OTU includes: The processing response direction of each of the cross-location stability indication OTUs in each location is transformed into a response feature vector; When two OTUs belong to at least the same family or order in taxonomy, and the Jaccard similarity coefficient of their response feature vectors is greater than 0.75, the two OTUs are aggregated into the same soil animal spatial variation module. After successful aggregation, the module abundance is calculated, which is the sum of the relative abundances of each OTU within the module; The aggregated modules are recursively executed to calculate the in-location processing response value and evaluate the cross-location stability response. Based on the response mode, the output results are divided into four categories: location response type, processing response type, location dependent type, and stability indication type.

Images