A method for constructing a soybean drought and high temperature combined resistance prediction model

Abstract

The application is a construction method of a soybean drought and high temperature combined stress resistance prediction model, and relates to the technical field of molecular stress resistance identification and molecular breeding, comprising: converting a normalized expression matrix into a recoverable proportion, dividing the recoverable proportion of a same batch reference variety HN65 to obtain a relative recoverable proportion, and outputting a relative recoverable proportion matrix. In the application, the individual's own baseline is taken as a reference to construct the recoverable proportion, and the degree of regression of the marker member to the individual's own baseline in the recovery period is associated with the deviation amplitude of the marker member in the stress period. Whether the marker member is up-regulated or down-regulated in the stress period, the recovery ability can be measured on the same scale without modeling according to the response direction.

Description

Technical Field

[0001] This invention relates to the field of plant stress resistance identification and molecular breeding technology, and in particular to a method for constructing a prediction model for combined drought and high temperature stress resistance in soybeans. Background Technology

[0002] Drought and high temperature are the main abiotic stresses restricting soybean yield and quality. In major soybean producing areas, drought and hot, dry winds often occur successively or simultaneously around the flowering period, and the yield loss caused by combined drought and high temperature stress is significantly higher than that caused by a single stress. Regarding the identification of soybean stress resistance, existing technologies have developed multiple routes: phenotypic identification based on field or pot water control experiments; evaluation based on physiological indicators such as germination rate, survival rate, yield, and antioxidant enzyme activity; and genotyping based on molecular markers and omics detection, providing a technical foundation for screening drought-resistant and heat-tolerant germplasm. However, combined drought and high temperature stress is not a simple superposition of two single stresses; their damage is synergistically amplified during the flowering period, and the recovery ability after stress relief is more directly related to the final grain yield and quality. How to distinguish the two sources of stress under combined stress, characterize the recovery dynamics after stress relief at the molecular level, and transform the identification results into stress resistance grades that can be used for early breeding selection remains a technical challenge for soybean stress resistance identification.

[0003] For example, in the prior art, patent application CN201510537521.5 discloses a method for identifying drought resistance of soybeans during germination using PEG-6000. It simulates drought stress by using a suitable concentration of PEG-6000 solution, compares the differences in morphological and physiological indicators of different soybean varieties during germination, and uses the membership function method to comprehensively evaluate drought resistance. The identification process is not limited by season, but its identification object is a single drought stress during germination, and the evaluation is based on static phenotypic indicators during the stress period. Its applicability to combined drought and high temperature stress during flowering and the recovery dynamics after stress relief needs to be expanded.

[0004] Existing technologies have been used extensively to assess the drought and heat resistance of soybeans and to evaluate the effects of combined crop stress, but the following problems remain unresolved:

[0005] Identification is mostly carried out for single stresses. Under combined drought and high temperature stress, it is difficult to distinguish the stress source attributes of the identification markers and there is a lack of molecular-level identification methods for joint stress resistance.

[0006] Existing indicators mostly characterize the degree of damage caused by stress or static genotypic differences, but they are difficult to characterize the recovery dynamics of regression to baseline after stress is relieved. There is a lack of homologous pairing between molecular regulatory information and physiological recovery outcome. Summary of the Invention

[0007] The purpose of this invention is to provide a method for constructing a combined drought and high temperature resistance prediction model for soybeans in order to solve the above-mentioned problems.

[0008] To achieve the above objectives, the present invention adopts the following technical solution:

[0009] A method for constructing a prediction model for the combined drought and high temperature resistance of soybeans includes:

[0010] Step 1, Design of a phased compound stress experiment: A combination of drought and high temperature stress was applied to soybean materials containing the reference variety HN65. Leaf samples were collected at the baseline sampling point T0 before the stress treatment, the sampling point Ts at the time of stress relief, the auxiliary sampling point T3 in the early recovery period, and the final judgment sampling point T7 in the recovery period. The samples were divided into two homologous pairs.

[0011] Step 2, Biomarker Combination Screening and Determination: Based on leaf sample expression data, recovery period-specific miRNA-mRNA biomarker pairs that are negatively correlated with target mRNA expression levels are screened to form biomarker combinations;

[0012] Step 3, dual-platform quantitative expression detection: The expression levels of marker members of the marker combination in the leaf sample were determined by quantitative PCR or microarray, and the normalized expression matrix was obtained.

[0013] Step 4, Construction of recoverable proportion function: The recoverable proportion is converted from the normalized expression matrix, divided by the recoverable proportion of the reference variety HN65 in the same batch to obtain the relative recoverable proportion, and the relative recoverable proportion matrix is ​​output.

[0014] Step 5, trajectory parameterization joint determination: normalized rebound rate and rebound amplitude parameters are calculated from the relative recoverable ratio matrix, and the comprehensive rebound rate and comprehensive rebound amplitude are obtained by weighted aggregation and the stress resistance level is divided. The boundary zone material is verified by the physiological composite recovery index of another homologous pair. Grade I or II materials are delineated as HN65 type stress-resistant germplasm.

[0015] Preferably, the stress treatment adopts a four-group design: the control group CK has normal water supply and normal temperature management throughout the process; the drought single stress group D is subjected to drought stress only; the high temperature single stress group H is subjected to high temperature stress only; the drought and high temperature combined stress group DH is subjected to both drought and high temperature stress simultaneously, and the drought and high temperature combined stress group DH first starts the drought treatment, and then adds the high temperature treatment after the relative soil moisture content drops to the set threshold.

[0016] Based on the differential expression results between the drought-single stress group D and the high-temperature-single stress group H, the stress source attribute was labeled for each marker:

[0017] Individuals who showed differential expression only in drought-single stress group D and not in high-temperature-single stress group H were identified as drought-responsive; individuals who showed differential expression only in high-temperature-single stress group H and not in drought-single stress group D were identified as high-temperature-responsive; and individuals who showed differential expression in both drought-single stress group D and high-temperature-single stress group H were identified as combined-responsive.

[0018] Preferably, the stress relief rehydration method is to irrigate once to 75% to 80% of the field water holding capacity, and the cooling method is to restore the temperature to 26 degrees Celsius during the day and 22 degrees Celsius at night;

[0019] When rehydration and cooling are not performed synchronously, a double-anchor-release rule is adopted:

[0020] The transition period starts with the first termination event and the recovery period ends with the subsequent termination events. The number of days in the recovery period is calculated from the zero point of the recovery period.

[0021] The time axis for calculating the recoverable proportion of drought-responsive biomarkers is anchored at the rehydration time, the time axis for calculating the recoverable proportion of high-temperature-responsive biomarkers is anchored at the cooling time, and the time axis for calculating the recoverable proportion of combined-responsive biomarkers is anchored at subsequent de-emergence events.

[0022] Preferably, the screening of the biomarker combination includes recovery period specificity determination and attribute labeling:

[0023] Candidate miRNAs were differentially expressed in comparisons of T3 relative to Ts or T7 relative to Ts to ensure that the biomarkers carry dynamic information of the recovery period rather than just information of the stress period;

[0024] Candidate miRNAs were differentially expressed between the reference variety HN65 and the susceptible control variety at T3 or T7 time points to ensure that the markers could distinguish between resistant and susceptible genotypes.

[0025] Candidate miRNAs were differentially expressed in comparison of Ts to T0 to confirm that they were drought or high-temperature response molecules rather than constitutively differentially expressed molecules.

[0026] Preferably, the biomarker combination consists of 8 to 15 miRNA-mRNA biomarker pairs, and the composition of the three types of biomarker pairs within the combination meets the requirement of class balance:

[0027] The proportion of drought-responsive biomarkers should be no less than 20%, the proportion of high-temperature-responsive biomarkers should be no less than 20%, and the proportion of combined-responsive biomarkers should be no less than 30%.

[0028] Preferably, the normalization adopts a dual-track normalization rule:

[0029] All markers in the miRNA quantification branch are normalized only with the miRNA internal reference, and all markers in the mRNA quantification branch are normalized only with the mRNA internal reference. The normalized expression levels of the two branches are calculated independently, and no absolute comparisons are made across branches. The correlation between miRNA members and mRNA members is always compared using the dimensionless quantity of recoverable proportion from step four.

[0030] The miRNA internal control was selected from miR156a and miR171a, with the geometric mean of their circulation thresholds used as the normalization benchmark. The mRNA internal control was selected from a combination of actin gene and elongation factor gene, with the geometric mean of their circulation thresholds used as the normalization benchmark.

[0031] Reactions with a cycle threshold greater than 35 are judged as undetectable, and the corresponding normalized expression level is marked as below the detection limit, denoted as DL marker. The DL marker is passed to step four along with the data.

[0032] Preferably, in step four, the formula for calculating the recoverable proportion is:

[0033] ;

[0034] in, Represents the recovery period The recoverable proportion of each phase, The value can be 3 or 7; Represents the recovery period Normalized expression levels for each time phase; Normalized expression level at time point Ts, representing the moment when the coercion is lifted; Normalized expression level at time T0 before stress treatment;

[0035] The calculation of the recoverable ratio sets the segmentation rules for the denominator protection:

[0036] The technical standard deviation of expression levels was estimated using the variation data from all technical repetitions in step three, and twice the technical standard deviation was taken as the protection threshold. ;

[0037] when and If the time exceeds the detection limit, switch to the logarithmic field equivalent form for calculation;

[0038] when When carrying the DL mark, the absolute change is calculated using the detection limit conversion value as the denominator. The recoverable proportion and its relative recoverable proportion obtained by branch calculation are marked as relative judgment data. They only participate in the relative comparison and ranking between materials in the same batch, and do not participate in the absolute magnitude judgment based on the value of 1 as the full recovery benchmark. When they participate in the weighted aggregation in step five, the weight of the corresponding marker member is halved.

[0039] Preferably, the weights of each marker member are determined by germplasm panel validation data:

[0040] Based on the absolute value of the correlation coefficient between the relative recoverable proportion of each marker member in the germplasm panel and the confirmed joint stress resistance phenotype in the field, the values ​​were weighted according to the normalized proportion.

[0041] Weights are subject to a balance constraint at the biomarker category level: the sum of the weights of the three biomarker members—drought-responsive, high-temperature-responsive, and combined-responsive—accounts for no less than 20%, 20%, and 30% of the total weights, respectively.

[0042] Preferably, in step five, the normalized rebound rate is obtained by phase parameterization of the relative recoverable ratio of the two recovery phases T3 and T7, which characterizes the process speed of the recovery motion.

[0043] The rebound amplitude parameter is taken as the relative recoverable proportion of the final judgment phase T7, which describes the final state of the recovery motion; the weighted aggregation is to weight and aggregate the normalized rebound rate and rebound amplitude parameter of all marker members according to the weights specified in step four, and the resulting comprehensive rebound rate and comprehensive rebound amplitude constitute the molecular recovery dual-parameter coordinates, with each material corresponding to a coordinate point in the two-dimensional parameter space.

[0044] Preferably, in step five, the boundary zone is defined as the area where the relative distance between the coordinate point formed by the material's overall springback rate and overall springback amplitude and the nearest grade threshold does not exceed 5%.

[0045] Materials falling into the boundary zone will undergo physiological channel review; materials not falling into the boundary zone will be directly classified according to the initial assessment level.

[0046] Physiological pathway verification is as follows:

[0047] Unfrozen leaf samples from the same source pairing and preservation step 1 were taken, and superoxide dismutase activity, peroxidase activity, catalase activity and malondialdehyde content were simultaneously measured at four time points: T0, Ts, T3 and T7. The physiological recovery ratio of each physiological index was calculated according to the recoverable ratio formula in step 4, and then aggregated into a physiological composite recovery index according to the physiological weights specified by the germplasm panel.

[0048] The rules for review and adjudication are as follows:

[0049] If the physiological composite recovery index of the material falling into the boundary zone reaches the physiological threshold corresponding to the next higher level, the level will be upgraded by one level.

[0050] If the physiological composite recovery index is lower than the lower limit of the physiological threshold corresponding to this level, the level will be downgraded by one level; otherwise, the initial level will remain unchanged.

[0051] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0052] 1. This invention constructs a recoverable proportion by using an individual's own baseline as a reference, and correlates the degree to which biomarker members regress to their own baseline during the recovery period with the magnitude of their deviation during the stress period. Regardless of whether biomarker members are upregulated or downregulated during the stress period, recovery capacity can be measured on the same scale without having to model separately according to the response direction.

[0053] 2. This invention uses a reference variety as a calibrator to compare the recoverable proportion of the material to be tested with that of the reference variety in the same batch. The differences in batch environment are offset during the comparison process, and the resulting relative recoverable proportion can be used across materials, batches, and even laboratories, providing a unified measurement basis for data exchange between different breeding units. Attached Figure Description

[0054] Further details, features, and advantages of this application are disclosed in the following description of exemplary embodiments in conjunction with the accompanying drawings, in which:

[0055] Figure 1 This is a flowchart of the method of the present invention;

[0056] Figure 2 This is a schematic diagram of the composite stress staged test design and sampling time axis of the present invention. Detailed Implementation

[0057] Several embodiments of this application will now be described in more detail with reference to the accompanying drawings to enable those skilled in the art to implement this application. This application may be embodied in many different forms and for various purposes and should not be limited to the embodiments set forth herein. These embodiments are provided to make this application thorough and complete, and to fully convey the scope of this application to those skilled in the art. The embodiments described do not limit this application.

[0058] Unless otherwise defined, all terms used herein (including technical and scientific terms) shall have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains. It will be further understood that terms such as those defined in commonly used dictionaries shall be interpreted as having a meaning consistent with their meaning in the relevant field and / or the context of this specification, and shall not be interpreted in an idealized or overly formal sense unless expressly defined herein.

[0059] Example 1

[0060] Its specific implementation method is combined with the appendix Figure 1 and Figure 2 Please provide a detailed explanation.

[0061] In this embodiment, it includes:

[0062] Step 1: Design of staged compound stress experiments:

[0063] The purpose of designing a multi-stress staging test is to provide a raw data foundation for the prediction model with a clear time axis, distinguishable stress sources, and homologous pairing of molecular and physiological data. This step directly determines the calibration quality of the subsequent recoverable proportional function and rebound trajectory parameters.

[0064] The experimental materials were set up in three levels:

[0065] The first level is the reference variety HN65. Heinong 65, which has passed drought resistance identification, was selected as the reference variety and is referred to as reference variety HN65. It has been confirmed to have outstanding performance in drought resistance identification during the bud and seedling stages. When drought stress is aggravated, the total root length, root surface area, root volume and root-shoot ratio show an increasing trend. It belongs to the germplasm with both tolerance and recovery mechanisms.

[0066] The second level consists of stress-sensitive control varieties, which are known stress-sensitive varieties that have experienced significant yield reduction under combined drought and high temperature stress and slow recovery after rehydration.

[0067] The third level consists of the germplasm population to be tested, including intermediate breeding materials, strains, or natural populations, with a minimum of 30 accessions. This is used for model parameter calibration and determination of grade thresholds. All materials must have plump seeds with a germination rate of no less than 95%, and must be carefully selected from the same batch before sowing to eliminate the interference of seed quality differences on baseline expression.

[0068] Cultivation conditions combined pot cultivation and artificial climate chambers. The cultivation substrate was prepared with peat moss, vermiculite, and perlite in a 3:1:1 volume ratio, with each pot containing the same amount of soil. Before sowing, the soil was uniformly irrigated to 75%-80% of field capacity. During the seedling stage, water and fertilizer were kept consistent, with slow-release compound fertilizer applied as a single basal application to avoid nutritional differences between individuals caused by topdressing. Each treatment of each material was planted with at least 15 pots, with 2 seedlings per pot, arranged in a completely randomized block design. The pots were rotated every 3 days to eliminate position effects.

[0069] The stress treatment period was selected during the flowering stage, namely the soybean R1 to R2 stage.

[0070] There are three reasons for choosing this period:

[0071] Firstly, the flowering period is the growth window where the combined stress of drought and high temperature causes the most severe yield loss. Pollen viability is sensitive to high temperature and pod formation is sensitive to water, and the combined damage is amplified synergistically during this period.

[0072] Secondly, during the flowering period, leaf metabolism is vigorous, the miRNA-mRNA regulatory network is active, and biomarker expression signals are strong.

[0073] Third, the recovery ability after rehydration during this period has the highest correlation with the final grain quality, ensuring that molecular indicators are linked to breeding outcomes.

[0074] The stress treatment adopted a four-group design: the control group CK had normal water supply and normal temperature management throughout the process; the drought single stress group D was subjected to drought only; the high temperature single stress group H was subjected to high temperature only; and the drought and high temperature combined stress group DH was subjected to both stresses simultaneously.

[0075] The significance of setting up groups D and H is to label the stress source attribute for each marker: those that are differentially expressed only in group D and not in group H are identified as drought-responsive; those that are differentially expressed only in group H and not in group D are identified as high-temperature-responsive; and those that are differentially expressed in both groups D and H are identified as joint-responsive. This attribute labeling provides a basis for the classification of marker combinations in step two, and also provides data for the calibration of marker class balance constraints in step four.

[0076] The drought stress was implemented using a combination of water cut-off and soil moisture monitoring. Irrigation was stopped from the start of the stress, and the relative soil moisture content was monitored every 12 hours using the dry weighing method or time domain reflectance probe. When the relative soil moisture content dropped to 35% and the deviation did not exceed 2 percentage points, the stress was determined to have reached the identification intensity and was maintained at this intensity. The maintenance method was to replenish water daily according to the evaporation limit to stabilize the moisture content within the set range, and the maintenance period was 7 days.

[0077] The implementation of high-temperature stress is as follows:

[0078] The temperature in the artificial climate chamber is raised to 38 degrees Celsius from 10:00 AM to 4:00 PM daily, and maintained at 30 degrees Celsius at night. The light cycle consists of 14 hours of light and 10 hours of darkness, with the relative humidity controlled between 55% and 65%, for a total of 7 days.

[0079] The operation sequence of the combined stress group DH was as follows: first, drought treatment was initiated, and after the relative soil moisture content dropped to the set threshold, high temperature treatment was added to simulate the real situation where drought arrived in the field before hot and dry winds; the stress maintenance period of the DH group was 7 days from the date of the high temperature addition.

[0080] The rules for coercion release operations and their time anchoring are as follows:

[0081] After the maintenance period, the re-watering method is to irrigate once to 75% to 80% of the field water holding capacity, and the cooling method is to restore the temperature to 26 degrees Celsius during the day and 22 degrees Celsius at night.

[0082] Considering that rehydration and cooling may not be synchronized in actual experiments, a dual release anchor point rule is established: the starting point of the transition interval is defined by the first release event, and the zero point of the recovery period is defined by the subsequent release event, with the recovery period days calculated from the zero point. For drought-responsive biomarkers, the time axis for calculating the recoverable proportion is corrected with the rehydration time as the anchor point; for high-temperature-responsive biomarkers, the time axis for calculating the recoverable proportion is corrected with the cooling time as the anchor point; and for combined-responsive biomarkers, the subsequent release event is used as the anchor point. This rule ensures that data from different batches of experiments with asynchronous release can be normalized to the same time axis.

[0083] There are four sampling time points: T0 is the baseline sampling point before stress treatment, which is sampled 1 day before the stress is initiated; Ts is the sampling point at the moment of stress relief, which is sampled within 2 hours after the rehydration or cooling operation is completed; T3 is the auxiliary sampling point in the early recovery period, which is sampled on the 3rd day from the zero point of the recovery period; T7 is the sampling point for the final judgment of the recovery period, which is sampled on the 7th day from the zero point of the recovery period.

[0084] All sampling was uniformly completed between 9:00 AM and 10:00 AM to eliminate the interference of circadian rhythms on miRNA and mRNA expression. The sampled tissue was the middle leaflet of the fully expanded uppermost trifoliate compound leaf. Each replicate consisted of at least 5 leaflets from the same leaf position. Immediately after sampling, the samples were flash-frozen in liquid nitrogen and then transferred to -80 degrees Celsius for storage.

[0085] Each sample was divided into two parts: one part was used for RNA extraction and expression quantification, and the other part was used directly for the determination of physiological indicators such as antioxidant enzyme activity without freezing. The two samples were homologously paired to ensure that the molecular and physiological data in step five came from the same biological population, thus preventing sample mismatch in molecular-physiological pairing. Biological replicates were set to at least three independent replicates for each material, each treatment, and each time point.

[0086] Environmental parameters were recorded throughout the process.

[0087] From 7 days before the stress was initiated to the 7th day of the recovery period, the relative soil moisture content, canopy temperature, air temperature and humidity were recorded daily to form a batch environmental record.

[0088] This file serves two purposes: first, it triggers the redoing of batch data when environmental parameters of a certain batch deviate from the set value beyond the allowable range; second, it investigates the cause of the batch's anomaly by combining the environmental parameter file when the recoverable proportion of a certain batch of reference variety HN65 deviates from its historical file range.

[0089] Step 2, screening and determining the combination of markers:

[0090] The purpose of biomarker combination screening is to identify a set of miRNA-mRNA biomarker pairs that are specific to the recovery period, have clear mechanisms, and are balanced in categories from the whole genome. This combination is the source of molecular information for the predictive model, and its screening criteria directly determine the mechanism specificity and cross-material applicability of the model.

[0091] In this scheme, a marker member refers to a single miRNA or mRNA molecule in the combination, and a marker pair refers to a combination of a miRNA member and an mRNA member with a negative regulatory pairing relationship.

[0092] The first stage of the screening was sequencing data collection. Leaf samples were collected from the two extreme materials, the reference variety HN65 and the stress-sensitive control variety, at four time points (T0, Ts, T3, and T7) obtained in step one. Each time point was repeated three times, and small RNA sequencing and transcriptome sequencing were performed simultaneously.

[0093] Small RNA sequencing libraries were ligated with adapters, reverse transcribed, and amplified before being sequenced, with read lengths set to cover small RNA fragments of 18 to 30 nucleotides. Transcriptome sequencing employed a conventional mRNA enrichment library construction strategy with paired-end sequencing. Both types of libraries were constructed using the same batch of homologous RNA samples to ensure strict temporal correspondence between miRNA and mRNA expression data.

[0094] The second stage of screening is data preprocessing and molecular identification. After adapter sequence removal and low-quality read filtering, the raw small RNA data retains sequences of 18 to 30 nucleotides in length, which are then aligned to the soybean reference genome.

[0095] Known miRNAs were identified by alignment with soybean entries in miRBase, while novel miRNAs were identified by precursor hairpin structure prediction. Prediction required both a minimum folding free energy threshold and support from an asterisk sequence. Transcriptome data were quality-controlled and aligned to the same reference genome for quantification at the gene level. miRNA and mRNA expression levels were normalized using transcript counts per million reads to eliminate sequencing depth differences.

[0096] The third stage of screening is the recovery period specificity determination, which consists of four criteria. Candidate miRNAs must simultaneously meet the first two criteria and complete the attribute labeling according to the last two criteria. Criterion 1 is recovery differential expression: the candidate miRNA must meet the following fold change in T3 relative to Ts or T7 relative to Ts. And the false detection rate is less than 0.05, of which The ratio of normalized expression levels represents the ratio of the two comparison objects. In criterion one, the two comparison objects are samples at the recovery period time point and samples at the time of stress relief, respectively. In criterion two, the two comparison objects are samples of the reference variety HN65 and samples of the susceptible control variety, respectively. This criterion ensures that the biomarker carries dynamic information of the recovery period rather than only information of the stress period.

[0097] Criterion 2 is genotypic heterogeneity: the differential expression of candidate miRNAs between the reference variety HN65 and the susceptible control variety at time points T3 or T7 also meets the criteria. Furthermore, with a false detection rate of less than 0.05, this criterion ensures that the biomarker has the ability to distinguish between resistant and susceptible genotypes.

[0098] Criterion 3 is stress-induced labeling: Candidate miRNAs are differentially expressed in comparison of Ts relative to T0, confirming them as drought or high-temperature response molecules rather than constitutively differentially expressed molecules. Criterion 4 is stress source labeling: Based on the differential expression results of groups D and H in step 1, candidate miRNAs are labeled as drought-responsive, high-temperature-responsive, or combined-responsive.

[0099] The fourth stage of screening is target gene prediction and pairing confirmation. For candidate miRNAs that have passed the third stage, the sequence complementarity scoring rule is used to predict their target mRNAs. The scoring rule takes into account the number of mismatches in complementary regions, the number of GU pairs, and the accessibility of target sites. It is supplemented by degradome sequencing data to verify the miRNA-mediated cleavage sites, and only miRNA-mRNA pairs supported by cleavage evidence are retained.

[0100] Subsequently, negative correlation verification was performed: the Pearson correlation coefficient between miRNA expression and target mRNA expression at four time points (T0, Ts, T3, and T7) for each pair was calculated, and only pairs with a correlation coefficient not exceeding -0.6 were retained. This threshold ensures that the negative regulatory relationship of the pairings truly exists in the recovery dynamics rather than being sequencing noise.

[0101] The formula for calculating the negative correlation coefficient is:

[0102] ;

[0103] in, Representing the The negative correlation coefficient of a pair of markers; This indicates that the miRNA in the marker pair is in the first stage. Normalized expression levels at each time point; This represents the average expression level of miRNA at all time points. This indicates that the marker targets the mRNA in the first phase. Normalized expression levels at each time point; The average expression level of the target mRNA at all time points; This represents the total number of time points involved in the calculation.

[0104] The fifth stage of screening involves prioritizing functional annotations. For confirmed miRNA-mRNA pairs, priority is given based on the functional annotations of the target genes.

[0105] Priority will be given to target gene annotations involving heat shock proteins, dehydration response element binding transcription factors, heat shock transcription factors, aquaporins, reactive oxygen species scavenging systems, and osmotic regulatory substance synthesis pathways.

[0106] Target genes annotated as having housekeeping functions and not related to stress responses are downgraded. As a preferred implementation, the candidate pool may include drought or high-temperature response miRNAs reported in previous studies, such as gma-miR166a and its target mRNA, miR398a and its target mRNA, miR156 family and its target SPL transcription factor, etc. However, the above-mentioned reported molecules are only used as candidate sources, and whether they can enter the final combination still needs to be tested by all the criteria in this step, and literature reports are not used as a direct basis for selection.

[0107] The sixth stage of screening was the panel validation and finalization of the combination. A panel containing 20 to 30 soybean germplasms was constructed. The field-confirmed combined stress resistance phenotypes of each germplasm had been pre-confirmed through two-year repeated identification of field yield loss rate and rehydration recovery rate. All germplasms in the panel were subjected to combined stress treatment according to step one, and samples T0, Ts, T3, and T7 were collected. The expression levels of miRNA and mRNA in the candidate pairs were detected by quantitative PCR. The correlation coefficient between the relative recoverable proportion of each pair (calculated according to the method in step four) and the physiological recovery indicators of each germplasm (i.e., the recovery range of antioxidant enzyme activity and malondialdehyde content during the recovery period as measured in step five) was calculated. Only pairs with an absolute correlation coefficient of not less than 0.5 were retained.

[0108] Simultaneously, the correlation coefficient between the relative recoverable proportion of each pair and the confirmed combined stress resistance phenotype in the field is calculated, and the obtained data is used for biomarker weighting in step four. The final number of combinations is 8 to 15 pairs, and the composition of the three types of biomarkers within the combination meets the balance requirement: the proportion of drought-responsive pairs is not less than 20%, the proportion of high-temperature-responsive pairs is not less than 20%, and the proportion of combined-response pairs is not less than 30%. This ensures that the information of the two types of stress in the combined stress has an independent carrier in the combination, and avoids the combination being dominated by a single stress signal and losing its "combined" identification ability.

[0109] The seventh stage of screening was stability testing. For all biomarker members of the finalized combination, intra-batch and inter-batch repeatability tests were performed: intra-batch testing required that the coefficient of variation of three technically repeated quantifications of the same RNA sample not exceed 5%; inter-batch testing required that the correlation coefficient of expression levels of the same material in two independent experimental batches not be less than 0.9.

[0110] The marker pair to which the marker member that fails the stability test belongs is removed from the portfolio, and candidate pairs are added according to the priority ranking of the fifth stage until all marker pairs in the portfolio pass the test.

[0111] The internal reference system was also determined concurrently with the combination:

[0112] miRNAs that are stably expressed across all time points and genotypes are screened from small RNA sequencing data as miRNA internal control candidates, and stably expressed mRNAs are screened from transcriptome data as mRNA internal control candidates, providing internal control sources for the dual-track normalized quantification in step three.

[0113] At this point, the molecular composition, stress source attributes, target gene mechanism annotation, and quantitative stability of the biomarker combination have all been determined, forming a list of biomarkers that can be directly used to establish detection methods, and we proceed to step three.

[0114] Step 3, Dual-platform expression quantification:

[0115] The purpose of dual-platform expression quantification is to transform the miRNA-mRNA biomarker combination determined in step two into time-aligned, cross-platform comparable, and cross-batch mergeable expression data, providing standardized input for the recoverable proportion calculation in step four.

[0116] This step sets up two parallel detection branches: quantitative PCR and a custom chip. The two branches can be used independently or converted to each other through a conversion factor to adapt to the equipment conditions of different breeding units.

[0117] The first step is RNA extraction and quality control. Take the leaf samples preserved in step one and freeze them in liquid nitrogen. Grind them into a fine powder in a liquid nitrogen environment. Extract them using a total RNA extraction reagent with a small RNA retention process. The key point is that the small RNA component must not be discarded during the extraction process, because conventional column extraction will retain molecules with less than 200 nucleotides. An extraction system that is clearly labeled to retain small RNA must be selected.

[0118] The extracted RNA was digested with DNase to remove residual genomic DNA. Quality control criteria were as follows: an A260 to A280 ratio between 1.8 and 2.1 as determined by UV spectrophotometry indicated acceptable protein contamination; and an RNA integrity score of at least 7 as determined by capillary electrophoresis indicated undegraded RNA. Samples failing any criterion were re-extracted; those still failing triggered the resampling procedure in step one. Each sample was divided into two portions after extraction, one for miRNA quantification and the other for mRNA quantification.

[0119] The second step is the miRNA quantification branch. Since miRNA molecules are only about 22 nucleotides long, conventional reverse transcription primers cannot bind directly. Therefore, stem-loop reverse transcription primers are used: the 5' end of the primer forms a stem-loop backbone, and the 3' end contains a sequence that is complementary to the end of the target miRNA by 6 to 8 nucleotides. After reverse transcription, an extended complementary strand is formed, and then quantitative PCR amplification is completed using a miRNA-specific upstream primer and a stem-loop universal downstream primer.

[0120] Primers for each miRNA marker member must be specifically validated: Amplification is performed separately using the target miRNA and its closely related family members (differences of 1 to 2 nucleotides) as templates, requiring a cycle threshold difference of at least 5 cycles between the target and closely related members to prevent cross-amplification. The quantitative PCR reaction system and procedure are performed according to the dye method, with each reaction set up in 3 technical replicates. The acceptance range for amplification efficiency is 90% to 110%, and a standard curve is constructed using 5-fold serial dilutions of the template for determination.

[0121] The melting curve must have a single peak; reactions with extraneous peaks are considered invalid, and primer dimers are investigated. The internal control for the miRNA branch uses the stable miRNA control combination screened in step two, stage seven. As a preferred implementation, miR156a and miR171a can be used as two internal controls, with the geometric mean of their cycle thresholds used as the normalization benchmark. The geometric mean is less sensitive to abnormal fluctuations in a single internal control compared to the arithmetic mean.

[0122] The third step is the mRNA quantification pathway. Homologous RNA samples are taken, and reverse transcription is initiated by mixing oligo-dT primers and random hexamer primers to obtain complementary DNA covering the full length. Specific primers for each mRNA marker member across exon linkers are designed. The cross-linker design can eliminate false positive amplification caused by residual genomic DNA.

[0123] The amplification system, number of technical replicates, amplification efficiency acceptance range, and melting curve criteria are all consistent with those of the miRNA branch. The internal control for the mRNA branch uses the stable mRNA internal control gene screened in step two. As a preferred implementation, a combination of actin and elongation factor genes can be used, with the geometric mean also used as the normalization benchmark.

[0124] The fourth step involves a dual-track normalization rule. miRNAs and mRNAs belong to different molecular types, with different reverse transcriptase systems, primer binding methods, and amplification kinetics. If the same set of internal controls is used for normalization, the difference in system efficiency between the two pathways will be introduced into the expression levels, making the expression levels of miRNAs and mRNAs incomparable on an absolute scale. This disrupts the same-scale calculation of the recoverable ratio between miRNAs and mRNAs in step four. Therefore, a dual-track normalization rule is established: all markers in the miRNA pathway are normalized only using the miRNA internal control, and all markers in the mRNA pathway are normalized only using the mRNA internal control. The normalized expression levels of the two pathways are calculated independently, without cross-branch absolute comparisons. The correlation comparison between miRNAs and mRNAs is always performed using the dimensionless recoverable ratio from step four. Since the recoverable ratio is a ratio structure within the same pathway, the internal control factor is naturally canceled out in the numerator and denominator, thus eliminating the influence of dual-track differences in principle.

[0125] The formula for calculating normalized expression level is:

[0126] ;

[0127] in, Normalized expression levels of representative marker members; The value represents the difference between the cyclic threshold of the target marker member in the test sample and the cyclic threshold of the geometric mean of the same internal reference, minus the corresponding difference of the calibration sample. This represents the quantitative cycle threshold. The rule for defining the detection limit is as follows: a reaction with a cycle threshold greater than 35 is considered undetectable, and the corresponding normalized expression level is marked as below the detection limit, denoted as the DL marker. This marker is passed along with the data to step four to trigger the corresponding branch of denominator protection.

[0128] The fifth step is the chip branch, which is a customized expression chip for the final biomarker combination. miRNA probes and mRNA probes are immobilized on the same chip: the miRNA probes are designed according to the full-length complementary sequence and the hybridization temperature window is optimized, and the mRNA probes are designed by selecting specific segments in the target gene transcript. Each biomarker member is equipped with no less than 3 repeat probes.

[0129] The samples were fluorescently labeled, hybridized, washed, and scanned to obtain the raw signals. Chip data normalization was performed in two stages: the first stage was intra-chip normalization, which used the signal of the exogenous dopant reference as a benchmark to correct the labeling and hybridization efficiency; the second stage was inter-chip normalization, which used quantile normalization to align the signal distributions of different chips. The chip expression level of each marker member was taken as the median of its repeat probe signal.

[0130] The sixth step is to establish the cross-platform conversion factor. To make the data from the qPCR and microarray platforms interchangeable and merging, a cross-platform conversion sample set is set up: at least 10 leaf RNA samples are selected, and the expression levels of their biomarkers cover the complete dynamic range from near the detection limit to high expression. The same sample is quantified on both the qPCR platform and the microarray platform.

[0131] For each biomarker member, a linear fit was performed with qPCR-normalized expression level as the independent variable and microarray-normalized expression level as the dependent variable. The fitted equation is as follows:

[0132] ;

[0133] in, Representing the Normalized expression levels of individual marker members on the chip platform; Representing the Cross-platform conversion slope of individual marker members; Representing the Normalized expression levels of each biomarker member using qPCR platform; Representing the The cross-platform conversion intercept of each biomarker member. The coefficient of determination for each biomarker member should not be lower than 0.98 to be acceptable. If it is lower than this value, check the chip probe or qPCR primer of that biomarker member, replace it and refit. and Once determined, the conversion parameters are archived as fixed conversion parameters for that biomarker member. When the chip batch, fluorescent labeling system, or qPCR reagent batch is changed, it must be recalibrated using a cross-platform converted sample set. Data obtained from any platform can be converted into equivalent data for another platform using conversion factors, enabling cross-platform merging modeling and result verification.

[0134] The seventh step is the data output rules. The data unit output for each material is a normalized expression matrix: the rows correspond to the marker members, the columns correspond to the four time points T0, Ts, T3, and T7, and the matrix elements are the normalized expression values.

[0135] The accompanying output information includes the DL label for each element, the technical repeatability coefficient of variation, and the detection platform to which it belongs. Elements with a technical repeatability coefficient of variation exceeding 5% are marked as low-confidence data and are retested before calculation in step four. This expression matrix is ​​the input for step four.

[0136] Step 4, Construction of the recoverable proportional function:

[0137] The purpose of constructing the recoverable proportion function is to convert the expression data output in step three into a recovery metric that is comparable across materials and batches, with the individual as the reference. This metric is the mathematical operand of the trajectory parameterization in step five, and its numerical stability and comparability directly determine whether the entire prediction model can be established.

[0138] For any biomarker member of any material, the expression levels at four time points are extracted from the normalized expression matrix in step three, and denoted as follows: Normalized expression level at time T0 before stress treatment; Normalized expression level at time point Ts, representing the moment when the coercion is lifted; Normalized expression level at time point T3, representing the auxiliary sampling point in the early recovery period; Normalized expression level at time T7, the final sampling point during the recovery period.

[0139] The four marker pairs are defined and assigned values ​​separately for miRNA and mRNA members, and the two members of each marker pair are calculated independently.

[0140] The recoverable proportion characterizes "the proportion of a marker member's regression to its baseline during the recovery period relative to its deviation amplitude during the stress period," calculated as follows:

[0141] ;

[0142] in, Represents the recovery period The recoverable proportion of each phase, The value can be 3 or 7; Represents the recovery period The normalized expression level of each phase, i.e. or ; Normalized expression level representing the moment the coercion is lifted; This represents the normalized expression level at baseline before stress treatment. The formula is directionally consistent: for marker members upregulated during stress, A negative value indicates that the expression level drops during the recovery period. If both are negative, dividing them gives a positive value;

[0143] For members with markers of downgraded stress period, When the expression level is positive, it rises during the recovery period. Both are positive, and the division is also positive. There is no need to model separately according to the response direction, which is a structural advantage of recoverable proportions compared to conventional expression ratio indicators.

[0144] An equal value of 1 indicates that the marker member is in the [number]th [year]. Each phase has returned to its baseline, meaning that the molecular-level recovery is complete; A value between 0 and 1 indicates partial recovery; the higher the value, the greater the degree of recovery. A value greater than 1 indicates over-rebound, meaning the recovery movement crosses the baseline and continues in the opposite direction, suggesting that the material has compensatory regulation in this marker member. It should be retained as a preferred record rather than truncated, as the over-rebound amplitude has identification value for some stress-resistant germplasm.

[0145] A value less than 0 indicates that the marker member failed to regress during the recovery period and instead deviated further from the baseline, suggesting a failure of the recovery mechanism. Each segment of this value range has a clear physiological meaning, making the model parameters interpretable.

[0146] When the denominator When the absolute value is close to zero, The minute measurement errors are amplified by the ratio structure, necessitating segmented protection. The protection threshold is determined by estimating the technical standard deviation of the expression level using the variation data from all technical repetitions in step three, and denoted as... To protect the threshold, take twice the technical standard deviation.

[0147] The segmentation rules have three branches:

[0148] Branch 1: When When that happens, the main formula of the second stage is used for calculation.

[0149] Branch 2: When and If the value exceeds the detection limit, switch to the logarithmic field equivalent form for calculation. The calculation formula is as follows:

[0150] ;

[0151] in, Represents the recovery period The recoverable proportion of each phase; , , The meaning is consistent with the main formula; This represents a calibration constant to prevent the logarithmic independent variable from being zero, and is taken as half of the minimum normalized expression level of all positive samples of this biomarker member. The logarithmic domain form converts multiplicative noise into additive noise, compressing the variance at the high expression end, ensuring that the recoverable proportion of low-amplitude response biomarker members maintains the same scale as the main formula and can be directly incorporated into subsequent weighting. Branch 3: When When the expression of the DL marker is carried, i.e., below the detection limit during the stress period, the absolute change is calculated using the detection limit as the denominator. The calculation formula is as follows:

[0152] ;

[0153] in, Represents the recovery period The recoverable proportion of each phase; Represents the recovery period Normalized expression levels for each time phase; This represents the normalized expression level corresponding to the detection limit. The recoverable proportion and its relative recoverable proportion calculated by branch three are marked as relative judgment data, which only participate in the relative comparison and ranking among materials in the same batch, and do not participate in the absolute magnitude classification with a value of 1 as the complete recovery benchmark; when participating in the weighted aggregation in step five, the weight of the corresponding marker member is halved. Each The calculation branches are recorded along with the data for review and traceability.

[0154] The three-branch setup resolves the ratio distortion of low-expression, low-response biomarkers in segments, providing a numerical stability guarantee that conventional ratio-based recovery indices do not possess.

[0155] The absolute recoverable proportion is still affected by batch-to-batch environmental differences; therefore, the reference variety HN65 is used as the accompanying calibrator. Each test batch must include complete processing and sampling of the reference variety HN65. The relative recoverable proportion is obtained by dividing the recoverable proportion of the test material by the recoverable proportion of the reference variety HN65 in the same batch. The calculation formula is as follows:

[0156] ;

[0157] in, Representing the The marker member during the recovery period The relative recoverable proportion of each time phase; Representing the material to be tested The marker member during the recovery period The recoverable proportion of each phase; This represents the reference variety HN65 from the same batch. The marker member during the recovery period The recoverable proportion of each phase. A value of 1 indicates that the material's recovery ability in this marker member is comparable to that of the reference variety HN65; a value greater than 1 indicates that it is stronger than the reference variety HN65; and a value less than 1 indicates that it is weaker than the reference variety HN65. When a sample falls within the denominator's protection range, the marker is downgraded within that batch and used only for relative comparisons between materials, without participating in absolute grading. After relative normalization, batch effects are canceled out in the numerator and denominator, allowing data to be combined across years and laboratories.

[0158] The contributions of each biomarker member to the combined stress resistance phenotype differ. The weights are determined from the germplasm panel validation data in Step 2, Stage 6: based on the absolute value of the correlation coefficient between the relative recoverable proportion of each biomarker member in the panel and the field-confirmed combined stress resistance phenotype, weights are assigned according to normalized proportions. The calculation formula is as follows:

[0159] ;

[0160] in, Representing the The weight of each marker member; Representing the The correlation coefficient between the relative recoverable proportion of each biomarker member in the germplasm panel and the field-confirmed combined stress resistance phenotype; This represents the total number of marker members participating in the aggregation, that is, the sum of all miRNA members and mRNA members in the final combination.

[0161] At the same time, the weights are subject to a balance constraint at the marker category level: the sum of the weights of the three marker members of drought response, high temperature response, and joint response shall account for no less than 20%, 20%, and 30% of the total weights, respectively, to prevent information of a certain type of stress from being marginalized in the weighting and to maintain the structural attribute of "joint" identification.

[0162] For each miRNA-mRNA marker pair, examine the sign consistency of the recoverable proportions of the two members:

[0163] Negative regulatory relationships require that the recovery phase miRNA and its target mRNA regress synchronously to their respective baselines, meaning that the recovery ratios of the two should be equal. If opposite signs occur, the biomarker pair is considered a regulatory mismatch in this material, and the weights of both members are halved before participating in subsequent calculations. The mismatch record is then output as auxiliary information indicating abnormal recovery of the regulatory network in this material. This validation filters out false positive signals caused by detection errors or individual regulatory dysregulation at the mechanistic level.

[0164] The data unit output from step four to step five is a relative recoverable proportion matrix: rows correspond to marker members, columns correspond to the two recovery phases T3 and T7, and the matrix elements are... The accompanying information includes the computational branch record for each element, DL tags, and weights. The consistency verification result with the other layer. This matrix is ​​the direct input for trajectory parameterization and joint determination in step five.

[0165] Step 5, Joint determination of trajectory parameters:

[0166] The purpose of trajectory parameterization joint determination is to transform the relative recoverable ratio matrix output in step four into a two-parameter feature that characterizes the recovery dynamics, complete the initial judgment of the stress resistance level, the boundary calibration of the physiological channel and the final delineation of the HN65 germplasm, and is the link in the prediction model outputting the identification conclusion.

[0167] For any marker member of any material, the normalized rebound rate is calculated by using the relative recoverable ratios of the two recovery phases, T3 and T7, as phase parameterization. The calculation formula is as follows:

[0168] ;

[0169] in, Representing the Normalized rebound rate of each marker member; Representing the The relative recoverable proportion of each marker member on day 7 of the recovery period; Representing the The relative recoverable proportion of each marker member on day 3 of the recovery period; This represents the time interval between two recovery phases, which is set to 4 days in this application. The physiological meaning of this parameter is the relative rate at which the marker recovers its own damage per unit time: a positive value indicates that the marker continues to regress to the baseline within the recovery window, and the larger the value, the faster the recovery; a value close to zero indicates that the recovery stalls in the T3 to T7 interval, i.e., plateauing or delayed initiation after early rebound; a negative value indicates that the marker deviates from the baseline in the opposite direction within the recovery window. This parameter must be based on the baseline correction and denominator protection in step four because the rate is a differential structure. The measurement noise of the two phases is amplified by subtraction, and the ratio of the original expression levels without stabilization is not discriminative after differential calculation. However, the relative recoverable ratio after baseline correction, logarithmic domain switching, and relative normalization with reference variety HN65 eliminates batch and internal reference factors in the numerator and denominator, allowing the differential results to be compared across materials.

[0170] The normalized rebound rate obtained in this way has a comparable recovery rate across genotype backgrounds and experimental batches, making recovery rate a quantitatively comparable identification dimension.

[0171] The rebound amplitude parameter is taken as the relative recoverable proportion of the final judgment phase, i.e. ,in, Representing the The rebound amplitude parameter of each marker member; Representing the The relative recoverable proportion of each marker member on day 7 of the recovery period. Describe the final state of the recovery movement and answer the question "How much has been recovered?" Describe the speed of the recovery process and answer the question, "At what relative speed does it recover?"

[0172] Both parameters are indispensable: (Looking only at...) It is impossible to distinguish between fast-rebound and slow-rebound germplasm, as both types may achieve the same amplitude on day 7; only by looking at... It is impossible to confirm whether the recovery has reached an effective level, because for high-rate materials, if the starting point is too low, the final recovery may still be insufficient. Only by combining two parameters can the recovery kinetics be fully characterized.

[0173] The two parameters of all marker members are weighted and aggregated according to the weights specified in step four. The calculation formula is as follows:

[0174] ;

[0175] ;

[0176] in, The overall rebound rate of the representative material; The overall rebound range of representative materials; Representing the The weight of each marker member; Representing the Normalized rebound rate of each marker member; Representing the The rebound amplitude parameter of each marker member; This represents the total number of marker members participating in the aggregation, that is, the sum of all miRNA members and mRNA members in the final combination. and Together they form the molecular recovery two-parameter coordinates, with each material corresponding to a coordinate point in the two-dimensional parameter space.

[0177] The grade threshold is marked on the germplasm panel in step two: based on the field-confirmed combined stress resistance phenotype of the panel materials, all panel materials are ( , The coordinate points are clustered according to the confirmed joint stress resistance level in the field, and the threshold for the five-level division is determined along the rate axis and the amplitude axis respectively.

[0178] The specific rules for determination are as follows:

[0179] Germplasm panel materials were divided into five phenotype groups based on the combined stress resistance phenotypes confirmed in the field. The overall rebound rate of materials in each phenotype group was calculated. Within-group mean and overall rebound magnitude The mean within a group is used as the midpoint between the mean values ​​of two adjacent phenotypic groups, and the threshold for dividing the adjacent levels is taken as the mean value within the group.

[0180] The five levels are: Level I, Not less than 1 and A recovery rate no lower than that of the same batch of the reference variety HN65 indicates that both the recovery endpoint and recovery rate have reached or exceeded those of the reference variety HN65, and the product is judged to be equivalent to or stronger than the reference variety HN65; Level II. Slightly below 1 and Approaching the level of the reference variety HN65, indicating that the recovery ability is close to that of the reference variety HN65; Level III. and All are between the Level II and Level IV thresholds, indicating moderate recovery ability; Level IV, Significantly lower than 1 or Approaching zero indicates a slow recovery; Level V, Below the threshold between Level IV and Level V or A negative value indicates that the recovery mechanism has failed, and the material is classified as sensitive. The initial judgment rule is: if the coordinates of the material fall within a certain level range, it is initially judged as belonging to that level; and For those belonging to different grade ranges, the lower grade is used as the initial assessment grade to reflect the constraint of the weakest link.

[0181] There is a boundary zone near the initial assessment threshold that is not robust. The boundary zone is defined as the area where the relative distance between the coordinate point and the nearest assessment threshold does not exceed 5%. Materials falling into the boundary zone will undergo physiological channel review; materials not falling into the boundary zone will be directly graded according to the initial assessment.

[0182] The data source for the physiological pathway is the unfrozen leaf samples that were preserved from the same source pair in step 1: superoxide dismutase activity, peroxidase activity, catalase activity and malondialdehyde content were measured simultaneously at four time points: T0, Ts, T3 and T7. The physiological recovery ratio of each physiological indicator was calculated according to the recoverable ratio formula in step 4, and then aggregated into a physiological composite recovery index according to the physiological weights specified on the panel.

[0183] Simultaneously, post-harvest grain quality indicators, including protein content and 100-grain weight retention rate, are measured as long-term validation items for recovery outcomes. The formula for calculating the physiological composite recovery index is:

[0184] ;

[0185] in, Represents the physiological composite recovery index; Representing the The weights of the physiological indicators were normalized and calibrated by the correlation coefficient between the recovery rate of the indicator in the germplasm panel and the confirmed combined stress resistance phenotype in the field. Representing the The recovery rate of each physiological indicator during the recovery period is calculated using the same structure as the aforementioned recoverable rate. The total number of physiological indicators involved in the aggregation. Physiological Complex Recovery Index. The physiological thresholds corresponding to each level are determined according to the same midpoint rule of the means of adjacent phenotype groups as the molecular two-parameter level thresholds.

[0186] The rules for review and adjudication are as follows:

[0187] If the material falling into the boundary zone, If the physiological threshold corresponding to the next higher level is reached, the level is raised by one level; if If the grade falls below the lower limit of the corresponding physiological threshold, the grade will be downgraded by one level; otherwise, the initial grade will remain unchanged. Grain quality indicators are not directly involved in the grading for the current year, but serve as the basis for the model to retrospectively verify the threshold in subsequent planting cycles.

[0188] The judgment unit for each material output includes: final stress resistance level, two-parameter coordinates ( , ), members of each marker and Detailed information, including biomarker pair regulatory mismatch records, physiological composite recovery index (if verification is initiated), and all calculated branches and quality control markers, is collected. Materials ultimately classified as Grade I or II are designated as HN65-type stress-resistant germplasm and will be used for breeding or further regional trials. Grade III materials will be retained or disposed of based on breeding objectives. Grade IV and V materials will be eliminated in the direction of stress-resistant breeding. Biomarker pair regulatory mismatch records serve as supplementary information, indicating abnormal recovery sites in the material's regulatory network and providing molecular-level references for parental selection.

[0189] Before the model is put into use, it must be independently validated: Select no less than 20 soybean materials with confirmed combined stress resistance phenotypes in the field and which have not participated in threshold calibration to form a validation group. The testing personnel who do not know the phenotypic information shall perform double-blind testing according to the complete process of steps one to five. The model's judgment level shall be compared with the confirmed combined stress resistance level in the field. The validation is considered successful if the judgment accuracy rate is no less than 85%.

[0190] The concordance rate of Grade I and Grade II materials is calculated separately, requiring a minimum of 90%, as the cost of false negatives in delineating stress-resistant germplasm is higher than that of misjudgments in intermediate grades. Repeatability verification requires that the judgment grade of the same material be consistent across three independent batches, and that the coefficient of variation of the two-parameter coordinates of the reference variety HN65 between batches not exceed 10%. If verification fails, the stability and weighting of the biomarker combination are retrospectively checked, and the combination is updated according to the supplementary rules in step two before recalibration and verification.

[0191] Model identification can be completed in one go during the flowering period, and the cycle from sampling to outputting the judgment conclusion is about 3 weeks, which is significantly shorter than the traditional process that relies on more than two years of field identification. In the breeding program, compound stress treatment and testing can be carried out on individual plants or small groups during the segregating generation, and early elimination can be implemented according to the judgment level. This allows molecular marker-assisted selection to advance from the level of "marker and genotype association" to the level of "restoration kinetic phenotypic prediction", forming a complementary rather than substitutive relationship with conventional phenotypic identification.

[0192] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for constructing a prediction model for the combined drought and high temperature resistance of soybeans, characterized in that, include: Step 1, Design of a phased compound stress experiment: A combination of drought and high temperature stress was applied to soybean materials containing the reference variety HN65. Leaf samples were collected at the baseline sampling point T0 before the stress treatment, the sampling point Ts at the time of stress relief, the auxiliary sampling point T3 in the early recovery period, and the final judgment sampling point T7 in the recovery period. The samples were divided into two homologous pairs. Step 2, Biomarker Combination Screening and Determination: Based on leaf sample expression data, recovery period-specific miRNA-mRNA biomarker pairs that are negatively correlated with target mRNA expression levels are screened to form biomarker combinations; Step 3, dual-platform quantitative expression detection: The expression levels of marker members of the marker combination in the leaf sample were determined by quantitative PCR or microarray, and the normalized expression matrix was obtained. Step 4, Construction of recoverable proportion function: The recoverable proportion is converted from the normalized expression matrix, divided by the recoverable proportion of the reference variety HN65 in the same batch to obtain the relative recoverable proportion, and the relative recoverable proportion matrix is ​​output. Step 5, trajectory parameterization joint determination: normalized rebound rate and rebound amplitude parameters are calculated from the relative recoverable ratio matrix, and the comprehensive rebound rate and comprehensive rebound amplitude are obtained by weighted aggregation and the stress resistance level is divided. The boundary zone material is verified by the physiological composite recovery index of another homologous pair. Grade I or II materials are delineated as HN65 type stress-resistant germplasm.

2. The method for constructing a combined drought and high temperature resistance prediction model for soybeans according to claim 1, characterized in that, The stress treatment adopted a four-group design: the control group CK was given normal water supply and normal temperature management throughout the process; the drought single stress group D was given only drought stress; the high temperature single stress group H was given only high temperature stress; and the drought and high temperature combined stress group DH was given both drought and high temperature stress simultaneously. In the drought and high temperature combined stress group DH, drought treatment was started first, and high temperature treatment was added after the relative soil moisture content dropped to the set threshold. Based on the differential expression results between the drought-single stress group D and the high-temperature-single stress group H, the stress source attribute was labeled for each marker: Individuals who showed differential expression only in drought-single stress group D and not in high-temperature-single stress group H were identified as drought-responsive; individuals who showed differential expression only in high-temperature-single stress group H and not in drought-single stress group D were identified as high-temperature-responsive; and individuals who showed differential expression in both drought-single stress group D and high-temperature-single stress group H were identified as combined-responsive.

3. The method for constructing a combined drought and high temperature resistance prediction model for soybeans according to claim 1, characterized in that, The method for restoring water after stress relief is to irrigate once to 75% to 80% of field capacity, and the method for cooling is to restore the temperature to 26 degrees Celsius during the day and 22 degrees Celsius at night; When rehydration and cooling are not performed synchronously, a double-anchor-release rule is adopted: The transition period starts with the first termination event and the recovery period ends with the subsequent termination events. The number of days in the recovery period is calculated from the zero point of the recovery period. The time axis for calculating the recoverable proportion of drought-responsive biomarkers is anchored at the rehydration time, the time axis for calculating the recoverable proportion of high-temperature-responsive biomarkers is anchored at the cooling time, and the time axis for calculating the recoverable proportion of combined-responsive biomarkers is anchored at subsequent de-emergence events.

4. The method for constructing a combined drought and high temperature resistance prediction model for soybeans according to claim 1, characterized in that, The screening of biomarker combinations includes determination of convalescent specificity and attribute labeling: Candidate miRNAs were differentially expressed in comparisons of T3 relative to Ts or T7 relative to Ts to ensure that the biomarkers carry dynamic information of the recovery period rather than just information of the stress period; Candidate miRNAs were differentially expressed between the reference variety HN65 and the susceptible control variety at T3 or T7 time points to ensure that the markers could distinguish between resistant and susceptible genotypes. Candidate miRNAs were differentially expressed in comparison of Ts to T0 to confirm that they were drought or high-temperature response molecules rather than constitutively differentially expressed molecules.

5. The method for constructing a combined drought and high temperature resistance prediction model for soybeans according to claim 4, characterized in that, The biomarker ensemble consists of 8 to 15 miRNA-mRNA biomarker pairs, and the composition of the three types of biomarker pairs within the ensemble meets the requirement of class balance. The proportion of drought-responsive biomarkers should be no less than 20%, the proportion of high-temperature-responsive biomarkers should be no less than 20%, and the proportion of combined-responsive biomarkers should be no less than 30%.

6. The method for constructing a combined drought and high temperature resistance prediction model for soybeans according to claim 1, characterized in that, Normalization adopts a dual-track normalization rule: All markers in the miRNA quantification branch are normalized only with the miRNA internal reference, and all markers in the mRNA quantification branch are normalized only with the mRNA internal reference. The normalized expression levels of the two branches are calculated independently, and no absolute comparisons are made across branches. The correlation between miRNA members and mRNA members is always compared using the dimensionless quantity of recoverable proportion from step four. The miRNA internal control was selected from miR156a and miR171a, with the geometric mean of their circulation thresholds used as the normalization benchmark. The mRNA internal control was selected from a combination of actin gene and elongation factor gene, with the geometric mean of their circulation thresholds used as the normalization benchmark. Reactions with a cycle threshold greater than 35 are judged as undetectable, and the corresponding normalized expression level is marked as below the detection limit, denoted as DL marker. The DL marker is passed to step four along with the data.

7. The method for constructing a combined drought and high temperature resistance prediction model for soybeans according to claim 1, characterized in that, In step four, the formula for calculating the recoverable proportion is: ； in, Represents the recovery period The recoverable proportion of each phase, The value can be 3 or 7; Represents the recovery period Normalized expression levels for each time phase; Normalized expression level at time point Ts, representing the moment when the coercion is lifted; Normalized expression level at time T0 before stress treatment; The calculation of the recoverable ratio sets the segmentation rules for the denominator protection: The technical standard deviation of expression levels was estimated using the variation data from all technical repetitions in step three, and twice the technical standard deviation was taken as the protection threshold. ; when and If the time exceeds the detection limit, switch to the logarithmic field equivalent form for calculation; when When carrying the DL mark, the absolute change is calculated using the detection limit conversion value as the denominator. The recoverable proportion and its relative recoverable proportion obtained by branch calculation are marked as relative judgment data. They only participate in the relative comparison and ranking between materials in the same batch, and do not participate in the absolute magnitude judgment based on the value of 1 as the full recovery benchmark. When they participate in the weighted aggregation in step five, the weight of the corresponding marker member is halved.

8. The method for constructing a combined drought and high temperature resistance prediction model for soybeans according to claim 7, characterized in that, The weights of each biomarker member were determined using germplasm panel validation data: Based on the absolute value of the correlation coefficient between the relative recoverable proportion of each marker member in the germplasm panel and the confirmed joint stress resistance phenotype in the field, the values ​​were weighted according to the normalized proportion. Weights are subject to a balance constraint at the biomarker category level: the sum of the weights of the three biomarker members—drought-responsive, high-temperature-responsive, and combined-responsive—accounts for no less than 20%, 20%, and 30% of the total weights, respectively.

9. The method for constructing a combined drought and high temperature resistance prediction model for soybeans according to claim 1, characterized in that, In step five, the normalized rebound rate is obtained by phase parameterization of the relative recoverable ratio of the two recovery phases T3 and T7, which characterizes the process speed of the recovery motion. The rebound amplitude parameter is taken as the relative recoverable proportion of the final judgment phase T7, which describes the final state of the recovery motion; the weighted aggregation is to weight and aggregate the normalized rebound rate and rebound amplitude parameter of all marker members according to the weights specified in step four, and the resulting comprehensive rebound rate and comprehensive rebound amplitude constitute the molecular recovery dual-parameter coordinates, with each material corresponding to a coordinate point in the two-dimensional parameter space.

10. The method for constructing a soybean drought and high temperature combined stress resistance prediction model according to claim 9, characterized in that, In step five, the boundary zone is defined as the area where the relative distance between the coordinate point formed by the material's overall springback rate and overall springback amplitude and the nearest grade threshold does not exceed 5%. Materials falling into the boundary zone will undergo physiological channel review; materials not falling into the boundary zone will be directly classified according to the initial assessment level. Physiological pathway verification is as follows: Unfrozen leaf samples from the same source pairing and preservation step 1 were taken, and superoxide dismutase activity, peroxidase activity, catalase activity and malondialdehyde content were simultaneously measured at four time points: T0, Ts, T3 and T7. The physiological recovery ratio of each physiological index was calculated according to the recoverable ratio formula in step 4, and then aggregated into a physiological composite recovery index according to the physiological weights specified by the germplasm panel. The rules for review and adjudication are as follows: If the physiological composite recovery index of the material falling into the boundary zone reaches the physiological threshold corresponding to the next higher level, the level will be upgraded by one level. If the physiological composite recovery index is lower than the lower limit of the physiological threshold corresponding to this level, the level will be downgraded by one level; otherwise, the initial level will be maintained.

Images