Rice nitrogen response regulation network analysis and breeding target identification system and method based on multi-omics data

Abstract

The invention discloses a rice nitrogen response regulation and control network analysis and breeding target identification system and method based on multi-omics data. According to the system, organic combination of regulation and control network construction based on single or multiple varieties of materials, key transcription factor recognition and accurate positioning of regulation and control areas where transcription factors play roles is achieved through an expression-chromatin accessibility correlation research method, and cis-trans effect distinguishing of the regulation and control areas is achieved through a deep learning model. The method comprises the following steps: carrying out nitrogen starvation pretreatment on rice, then carrying out nitrogen resupply, collecting a root sample, and carrying out ATAC-seq and RNA-seq sequencing; an eCAAS method is adopted to construct a regulation and control network, and key transcription factors are identified and accurately positioned; the chromatin accessibility difference of different varieties is predicted through a deep learning model, the cis-action effect and the trans-action effect are distinguished, an upstream transcription factor target is provided for genes dominated by the trans-effect, and haplotype and editable regulatory region targets available for direct breeding are provided for genes dominated by the cis-effect.

Description

Technical Field

[0001] The present invention relates to the technical field of plant molecular biology and bioinformatics, and in particular to a system and method for analyzing rice nitrogen response regulatory networks and identifying breeding targets based on multi-omics data. Background Art

[0002] Nitrogen is a key nutrient limiting rice growth and yield. Global rice production consumes approximately 100 million tons of nitrogen fertilizer annually, but nitrogen use efficiency is generally low, with only 30-40% of applied nitrogen fertilizer being absorbed and utilized by the crop, while the remainder is lost to the environment, causing serious pollution. Improving rice nitrogen use efficiency is crucial for ensuring food security and environmental sustainability.

[0003] The construction of plant nitrogen response regulatory networks aims to identify key regulatory genes that control nitrogen absorption, transport, and assimilation, providing molecular targets for variety improvement. However, existing technologies have the following shortcomings:

[0004] First, traditional gene expression analysis methods can only identify expressed genes and are unable to effectively identify key upstream transcription factors. In complex regulatory networks, it is often a small number of key transcription factors at the top of the regulatory hierarchy that truly determine trait differences. Identifying these transcription factors is crucial for variety improvement. Literature shows that improving a single key transcription factor can often influence the expression of dozens or even hundreds of downstream genes, resulting in significant phenotypic effects.

[0005] Second, existing methods lack the ability to precisely locate regulatory regions. Even if important regulatory relationships are identified, the precise location of these regions cannot be determined, severely limiting the application of precise molecular modification technologies such as gene editing. In current gene editing practices, the lack of precise target information often requires researchers to conduct extensive trial-and-error experiments, which is inefficient and costly.

[0006] Third, there is a lack of systematic methods to distinguish the root causes of differences in gene expression. When differences in target traits are discovered between varieties, it is difficult to determine whether to improve the regulatory region of the target gene (cis-effects) or the upstream transcription factor (trans-effects), resulting in a lack of scientific basis for formulating molecular improvement strategies. This blind approach not only wastes research resources but also reduces the success rate of molecular breeding.

[0007] Fourth, methods for integrating and analyzing multi-omics data are not mature enough. Existing technologies often analyze single-type omics data in isolation, failing to fully exploit the synergistic information from multi-omics data, thus affecting the accuracy and reliability of regulatory network construction.

[0008] Therefore, there is an urgent need to develop a technical solution that can systematically integrate multi-omics time-series data, accurately identify key transcription factors, locate regulatory regions, and formulate personalized improvement strategies. Summary of the Invention

[0009] This paper addresses the complexities of analyzing rice nitrogen response regulatory mechanisms by providing a multi-omics data-based system and method for analyzing rice nitrogen response regulatory networks and identifying breeding targets. This paper innovatively integrates multiple technologies, including time-series multi-omics data analysis, eCAAS regulatory network construction, deep learning prediction, and effect differentiation, to establish a comprehensive technical framework and achieve a systematic solution from regulatory network construction to precise molecular target localization.

[0010] To achieve the above purpose, the technical solution designed by the present invention is as follows:

[0011] The present invention provides a method for constructing a rice nitrogen response regulatory network based on multi-omics time series data, comprising the following steps:

[0012] S1. Multi-omics time-series data collection: Seedlings of a single or multiple rice varieties are placed in a nitrogen-deficient nutrient solution for 7 days of nitrogen starvation pretreatment. The seedlings are then transferred to a normal nutrient solution for culture. Root samples are collected at ≥10 time points within 48 hours. The root samples at each time point are divided into two equal parts for ATAC-seq and RNA-seq sequencing to obtain chromatin accessibility and gene expression data for the root samples at each time point. The collection of the two data at all time points is the multi-omics time-series data;

[0013] S2. Construction of regulatory network and identification of key transcription factors: The expression-chromatin accessibility association study eCAAS method was used to associate the dynamic changes of transcription factor expression with the chromatin accessible regions of the whole genome through a linear mixed model, and the chromatin accessible regions associated with transcription factors were dynamically associated with the expression of the genes closest to them through a linear mixed model to construct a nitrogen-responsive transcriptional regulatory network, that is, a regulatory network including transcription factors, target chromatin accessible regions (regulatory regions) and target genes.

[0014] Furthermore, in step S1, after the ATAC-seq and RNA-seq raw sequencing data undergo data quality control, the ATAC-seq data uses the MACS2 tool to identify chromatin accessible regions, and the RNA-seq data uses the salmon tool to quantify gene expression.

[0015] The transcription factor lists identified in the expressed genes were obtained based on the collection of reported rice transcription factor lists.

[0016] The above data quality control requirements are: the TSS enrichment score of ATAC-seq data is greater than 5, and the number of gene detections of RNA-seq data is greater than 15,000.

[0017] Furthermore, in step S2, the eCAAS method includes two levels of correlation analysis:

[0018] a. Transcription factor-chromatin accessible region (regulatory region) association analysis: A linear mixed-effects model was established, with transcription factor expression as the independent variable and the degree of chromatin accessible region openness as the dependent variable:

[0019] A ijk =β0+β1Y ik +Δγ jk +ε ijk

[0020] Among them, A ijk represents the accessibility level of chromatin accessible region j in sample i and time point k, β0 represents the fixed effect intercept term, β1 represents the fixed effect slope coefficient, and Y ik represents the expression level of transcription factor in sample i and time point k, γ jk is a random effect term, which obeys the multivariate normal distribution N(0,∑ RNA ), ε ijk represents the residual error term;

[0021] The variance-covariance matrix ∑ of the random effect term RNA Sample correlation calculation by RNA-seq data:

[0022]

[0023] Where n represents the total number of samples measured in the RNA-seq data, X g Represents the expression vector of gene g in all samples, is the average expression level of gene g in all samples, and G is the total number of genes;

[0024] b. Regulatory region-target gene association analysis: A linear mixed-effects model was established, with the degree of openness of the chromatin accessible region obtained in step a as the independent variable and the expression level of the gene closest to the chromatin accessible region as the dependent variable:

[0025] E ijk =α0+α1A jk +δ ik +∈ ijk

[0026] Among them, E ijkrepresents the expression level of target gene i in sample k and time point j, α0 is the fixed effect intercept term, α1 is the fixed effect slope coefficient, and A jk represents the chromatin accessibility level of the regulatory region at sample k and time point j, δ ik is a random effect term, which obeys the multivariate normal distribution N(0,∑ ATAC ), ε ijk represents the residual error term;

[0027] The variance-covariance matrix ∑ of the random effect term ATAC Sample correlation calculation by ATAC-seq data:

[0028]

[0029] Where n represents the total number of samples measured in the ATAC-seq data, and Z r represents the signal intensity vector of the rth chromatin region in all samples, The average signal intensity of the rth chromatin accessible region, R is the total number of regulatory regions.

[0030] The present invention also provides an application of a regulatory network constructed by the above method in identifying key transcription factors and the positioning of their regulatory regions, as well as in analyzing variety differences and formulating improvement strategies.

[0031] According to the actual situation, experimental verification of key transcription factors is carried out, including CUT&Tag chromatin immunoprecipitation experiment to verify the direct binding relationship between transcription factors and regulatory regions, EMSA gel migration experiment to verify the binding ability of transcription factors to DNA motifs, and dual luciferase reporter gene experiment to verify the functional effect of regulatory relationship.

[0032] The present invention also provides a method for identifying key transcription factors and the location of their regulatory regions. The method comprises the following steps: using the nitrogen-responsive transcriptional regulatory network constructed by the above method, and employing a network topology analysis method to calculate the regulatory importance score of each transcription factor; sorting the factors according to the importance scores, identifying the top 10% of the transcription factors as key transcription factors, and simultaneously accurately locating the position of the chromatin regulatory region corresponding to each regulatory relationship.

[0033] Furthermore, the formula for calculating the regulatory importance score of each transcription factor is as follows:

[0034]

[0035] Among them, IS TFk Transcription factor TF k Importance score, DC TFk is degree centrality, BC TFk is the betweenness centrality, CCTFk is the closeness centrality, α, β, and γ are weight coefficients, and their values ​​are 0.5, 0.3, and 0.2 respectively.

[0036] The present invention also provides a method for analyzing rice variety differences and identifying breeding targets, comprising the following steps:

[0037] 1) Obtaining the genome sequence of the target variety material to be analyzed and obtaining multi-omics time series data of the target variety material to be analyzed based on the method of step S1 of the above method;

[0038] 2) Determine candidate regulatory regions using the regulatory network constructed using the above method. Use a deep learning model to predict the chromatin accessibility of the candidate regulatory regions of the target variety to be analyzed, and perform differential analysis with the provided ATAC-seq data. Based on the degree of difference between the predicted chromatin accessibility of the target variety and the chromatin accessibility in the provided ATAC-seq data, distinguish the contributions of cis- and trans-acting effects, determine breeding targets based on the type of effect, and provide variety improvement strategies:

[0039] If the candidate regulatory region is identified as a trans-acting effector, a list of candidate upstream transcription factors recommended for variety improvement is output;

[0040] Alternatively, if it is a cis-acting effect, the location of this regulatory region can be directly output for variety improvement.

[0041] Furthermore, in step 2), a deep learning model is constructed based on an improved Basenji framework, including a convolutional neural network module and a multi-task learning framework, to predict chromatin accessibility at different time points. The loss function of the deep learning model is:

[0042]

[0043] Among them, L accessibility Represents the loss function value of chromatin accessibility prediction, which measures the difference between the predicted result and the true value. N is the total number of samples. is the predicted accessibility value of the i-th chromatin region, is the true accessibility value of the i-th chromatin accessible region;

[0044] Identify regulatory regions associated with any transcription factor in the target variety, predict the chromatin accessibility of this regulatory region based on the above deep learning model, and compare the predicted results with the chromatin accessibility results of the regulatory region in the ATAC-seq data provided by the target variety to distinguish cis- and trans-effects. The quantitative indicators are:

[0045]

[0046] Among them, CTE geneiis the cis-trans effect ratio of gene i, A pred,target The chromatin accessibility degree predicted by the model for the target species, A exp,target is the degree of chromatin accessibility obtained from the ATAC-seq data of the target species.

[0047] CTE-based genei The specific value of is as follows:

[0048] When CTE genei When ≥0.8, it was determined to be cis-effect dominant;

[0049] Or, when CTE genei When <0.5, it was determined to be dominated by the trans effect;

[0050] Or, when 0.5≤CTE genei When <0.8, it was determined to be a mixed effect.

[0051] Furthermore, in step 2), breeding targets are determined based on the effect type determination results, and variety improvement strategies are formulated:

[0052] For regulatory regions dominated by trans effects, output a list of upstream transcription factors associated with the regulatory region, and the expression difference of each transcription factor compared with the corresponding transcription factor in the reference variety; output in descending order based on the regulatory strength of the transcription factors with the regulatory region dominated by the trans effect;

[0053] For regulatory regions dominated by cis effects, accurate regulatory region coordinate information is output as gene editing targets;

[0054] For regulatory regions dominated by trans effects, output the transcription factors associated with the regulatory regions. The recommended transcription factors are ranked based on the regulatory strength score:

[0055] RS target =w1·|FC|+w2·(-log 10 (P value ))+w3·CS

[0056] Among them, RS target is the recommendation score of the transcription factor, FC is the maximum expression change fold of the transcription factor in the time series data, Pvalue is the eCAAS significance P value of the transcription factor and the selected regulatory region, CS is the interaction score of the homologous "upstream transcription factor-target" gene pair in other species such as Arabidopsis, w1, w2, and w3 are weight coefficients, which are 0.4, 0.4, and 0.2 respectively;

[0057] Alternatively, if there is no other species information, the CS weight is set to 0.

[0058] The present invention also provides a rice nitrogen response regulatory network analysis and breeding target identification system based on multi-omics data for implementing the above method, comprising:

[0059] Data preprocessing module, used to preprocess multi-omics time series data, including data quality control, format conversion and standardization functions;

[0060] A network construction module, used to construct transcriptional regulatory networks and identify key transcription factors, enabling efficient parallel computing of the eCAAS algorithm;

[0061] Effect differentiation module, used to analyze regulatory differences between varieties and distinguish effect types, integrating deep learning prediction models;

[0062] The target recommendation module is used to formulate improvement strategies and recommend molecular targets based on effect types, providing personalized improvement plans.

[0063] Beneficial effects of the present invention:

[0064] 1. Systematic identification capability of key transcription factors: Through the eCAAS method and network topology analysis, the present invention can systematically identify transcription factors that play a key role in the nitrogen response regulatory network.

[0065] 2. Ability to formulate personalized molecular modification strategies: This invention can develop targeted modification strategies for different gene types based on the differentiation of cis- and trans-effects. For genes dominated by trans-effects, upstream key transcription factors are provided as modification targets; for genes dominated by cis-effects, precise regulatory region locations are provided as gene editing targets. This personalized strategy-making capability significantly improves the targetedness and success rate of molecular modification, avoiding the blind trial and error of traditional methods.

[0066] 3. Precise Localization of Regulatory Regions: Compared to traditional methods that can only identify regulatory relationships, this method can precisely locate the chromatin regulatory region corresponding to each regulatory relationship, with a localization accuracy of 128bp. This precise localization capability provides important technical support for precise gene editing, transforming molecular modification from a "needle in a haystack" approach to targeted modification based on precise coordinates.

[0067] 4. Broad Application Value: The key transcription factors and regulatory regions identified by this invention can be widely applied in various fields, including variety breeding, precision gene editing, and molecular mechanism analysis. The established technical framework has good scalability and can be extended to analyze regulatory networks in other crops and other biological processes.

[0068] In summary: The present invention achieves systematic identification of key transcription factors through the innovative eCAAS method, and accurately distinguishes cis- and trans-effects through a deep learning model, thereby providing targeted molecular improvement strategies for different types of genes. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 This is a schematic diagram of the composition of the rice nitrogen response regulatory network analysis system;

[0070] Figure 2 A schematic diagram of the process for constructing the rice nitrogen response regulatory network, identifying its key transcription factors, analyzing cultivar differences, and developing improvement strategies;

[0071] Figure 3 Diagram of the experimental design for time-series multi-omics data collection;

[0072] Figure 4 Schematic diagram of building the regulatory network for the eCAAS approach. DETAILED DESCRIPTION

[0073] The present invention is further described in detail below with reference to specific embodiments so that those skilled in the art can understand.

[0074] Example 1: Rice nitrogen response regulatory network analysis and breeding target identification system based on multi-omics data

[0075] like Figure 1 The rice nitrogen response regulatory network analysis and breeding target identification system based on multi-omics data shown includes:

[0076] Data preprocessing module, used to preprocess multi-omics time series data, including data quality control, format conversion and standardization functions;

[0077] A network construction module, used to construct transcriptional regulatory networks and identify key transcription factors, enabling efficient parallel computing of the eCAAS algorithm;

[0078] Effect differentiation module, used to analyze regulatory differences between varieties and distinguish effect types, integrating deep learning prediction models;

[0079] The target recommendation module is used to formulate improvement strategies and recommend molecular targets based on effect types, providing personalized improvement plans.

[0080] Based on the rice nitrogen response regulatory network analysis and breeding target identification system, there are three main usage scenarios: Figure 2 ):

[0081] 1. A method for constructing a rice nitrogen response regulatory network based on multi-omics time series data, comprising the following steps:

[0082] S1. Multi-omics time series data acquisition ( Figure 3 ): Seedlings of a single or multiple rice varieties are placed in a nitrogen-deficient nutrient solution for nitrogen starvation pretreatment for 7 days, and then transferred to a normal nutrient solution for culture. Root samples are collected at ≥10 time points within 48 hours. The root samples at each time point are divided into two equal parts for ATAC-seq and RNA-seq sequencing to obtain chromatin accessibility and gene expression data for the root samples at each time point. The collection of the two data from all time points is the multi-omics time series data; among them, after data quality control of the original ATAC-seq and RNA-seq sequencing data, the ATAC-seq data uses the MACS2 tool to identify chromatin accessible regions, and the RNA-seq data uses the salmon tool to quantify gene expression.

[0083] The transcription factor lists identified in the expressed genes were obtained based on the collection of reported rice transcription factor lists.

[0084] The above data quality control requirements are: the TSS enrichment score of ATAC-seq data is greater than 5, and the number of gene detections of RNA-seq data is greater than 15,000.

[0085] S2. Regulatory network construction and identification of key transcription factors: using the eCAAS method for expression-chromatin accessibility association study ( Figure 4 ), using a linear mixed model to correlate the expression of transcription factors with the dynamic changes in chromatin accessible regions across the genome, and also using a linear mixed model to dynamically correlate the expression of chromatin accessible regions associated with transcription factors with the expression of the genes closest to them, to construct a nitrogen-responsive transcriptional regulatory network, that is, a regulatory network consisting of transcription factors, target chromatin accessible regions (regulatory regions), and target genes. The eCAAS method includes two levels of association analysis:

[0086] a. Transcription factor-chromatin accessible region (regulatory region) association analysis: A linear mixed-effects model was established, with transcription factor expression as the independent variable and the degree of chromatin accessible region openness as the dependent variable:

[0087] A ijk =β0+β1Y ik +γ jk +ε ijk

[0088] Among them, A ijk represents the accessibility level of chromatin accessible region j in sample i and time point k, β0 represents the fixed effect intercept term, β1 represents the fixed effect slope coefficient, and Y ik represents the expression level of transcription factor in sample i and time point k, γ jk is a random effect term, which obeys the multivariate normal distribution N(0,∑RNA ), ε ijk represents the residual error term;

[0089] The variance-covariance matrix of the above random effect terms ∑ RNA Sample correlation calculation by RNA-seq data:

[0090]

[0091] Where n represents the total number of samples measured in the RNA-seq data, X g Represents the expression vector of gene g in all samples, is the average expression level of gene g in all samples, and G is the total number of genes;

[0092] b. Regulatory region-target gene association analysis: A linear mixed-effects model was established, with the degree of openness of the chromatin accessible region obtained in step a as the independent variable and the expression level of the gene closest to the chromatin accessible region as the dependent variable:

[0093] E ijk =α0+α1A jk +δ ik +∈ ijk

[0094] Among them, E ijk represents the expression level of target gene i in sample k and time point j, α0 is the fixed effect intercept term, α1 is the fixed effect slope coefficient, and A jk represents the chromatin accessibility level of the regulatory region at sample k and time point j, δ ik is a random effect term, which obeys the multivariate normal distribution N(0,∑ ATAC ), ε ijk represents the residual error term;

[0095] Variance-covariance matrix of random effect terms∑ ATAC Sample correlation calculation by ATAC-seq data:

[0096]

[0097] Where n represents the total number of samples measured in the ATAC-seq data, and Z r represents the signal intensity vector of the rth chromatin region in all samples, The average signal intensity of the rth chromatin accessible region, R is the total number of regulatory regions.

[0098] 2. A method for identifying key transcription factors and their regulatory region locations is to use the nitrogen-responsive transcriptional regulatory network constructed by the above method, calculate the regulatory importance score of each transcription factor using a network topology analysis method; sort the transcription factors according to the importance score, identify the top 10% of the transcription factors as key transcription factors, and accurately locate the chromatin regulatory region corresponding to each regulatory relationship; wherein,

[0099] The formula for calculating the regulatory importance score of each transcription factor is as follows:

[0100]

[0101] Among them, IS TFk Transcription factor TF k Importance score, DC TFk is degree centrality, BC TFk is the betweenness centrality, CC TFk is the closeness centrality, α, β, and γ are weight coefficients, and their values ​​are 0.5, 0.3, and 0.2 respectively.

[0102] 3. A method for analyzing rice variety differences and formulating improvement strategies, comprising the following steps:

[0103] 1) Obtaining the genome sequence of the target variety material to be analyzed and obtaining multi-omics time series data of the target variety material to be analyzed based on the method of step S1 of the above method;

[0104] 2) The regulatory network constructed using the above method was combined to identify candidate regulatory regions. The chromatin accessibility of the candidate regulatory regions of the target species to be analyzed was predicted using a deep learning model, and differential analysis was performed with the provided ATAC-seq data. The deep learning model was built based on an improved Basenji framework, including a convolutional neural network module and a multi-task learning framework, to predict chromatin accessibility at different time points. The loss function of the deep learning model was:

[0105]

[0106] Among them, L accessibility Represents the loss function value of chromatin accessibility prediction, which measures the difference between the predicted result and the true value. N is the total number of samples. is the predicted accessibility value of the i-th chromatin region, is the true accessibility value of the i-th chromatin accessible region;

[0107] Identify regulatory regions associated with any transcription factor in the target variety, predict the chromatin accessibility of this regulatory region based on the above deep learning model, and compare the predicted results with the chromatin accessibility results of the regulatory region in the ATAC-seq data provided by the target variety to distinguish cis- and trans-effects. The quantitative indicators are:

[0108]

[0109] Among them, CTE genei is the cis-trans effect ratio of gene i, A pred,target The chromatin accessibility degree predicted by the model for the target species, A exp,target is the degree of chromatin accessibility obtained from the ATAC-seq data of the target species;

[0110] CTE-based genei The specific value of is as follows:

[0111] When CTE genei When ≥0.8, it was determined to be cis-effect dominant;

[0112] Or, when CTE genei When <0.5, it was determined to be dominated by the trans effect;

[0113] Or, when 0.5≤CTE genei When <0.8, it was determined to be a mixed effect.

[0114] Based on the degree of difference between the predicted chromatin accessibility of the target variety and the chromatin accessibility in the provided ATAC-seq data; distinguish the contribution of cis-acting effects and trans-acting effects, determine breeding targets based on the effect type, and provide variety improvement strategies:

[0115] If the candidate regulatory region is identified as a trans-acting effector, a list of candidate upstream transcription factors recommended for variety improvement is output;

[0116] Alternatively, if it is a cis-acting effect, the location of this regulatory region can be directly exported for variety improvement;

[0117] The specific methods are as follows:

[0118] For regulatory regions dominated by trans effects, output a list of upstream transcription factors associated with the regulatory region, and the expression difference of each transcription factor compared with the corresponding transcription factor in the reference variety; output in descending order based on the regulatory strength of the transcription factors with the regulatory region dominated by the trans effect;

[0119] For regulatory regions dominated by cis effects, accurate regulatory region coordinate information is output, and the haplotypes of the target varieties can be directly used for breeding or as gene editing targets;

[0120] For regulatory regions dominated by trans effects, output the transcription factors associated with the regulatory regions. The recommended transcription factors are ranked based on the regulatory strength score:

[0121] RS target =w1·|FC|+w2·(-log 10 (P value ))+w3·CS

[0122] Among them, RS target is the recommendation score of the transcription factor, FC is the maximum expression change fold of the transcription factor in the time series data, Pvalue is the eCAAS significance P value of the transcription factor and the selected regulatory region, CS is the interaction score of the homologous "upstream transcription factor-target" gene pair in other species such as Arabidopsis, w1, w2, and w3 are weight coefficients, which are 0.4, 0.4, and 0.2 respectively;

[0123] Alternatively, if there is no other species information, the CS weight is set to 0.

[0124] Example 2 Analysis of nitrogen response regulation network of Zhenshan 97 and Nipponbare

[0125] Based on the above method, the nitrogen response regulatory network of Zhenshan 97 and Nipponbare was constructed as follows:

[0126] The indica rice variety Zhenshan 97 (ZS97) and the japonica rice variety Nipponbare (NIP) were selected as research materials. After germination, seeds were cultured for three weeks in a standard rice nutrient solution composed of 1.44 mM NH₄NO₃, 0.3 mM NaH₂PO₄, 0.5 mM K₂SO₄, 1.0 mM CaCl₂, 1.6 mM MgSO₄, 0.17 mM Na₂SiO₃, 50 μM Fe(II)-EDTA, and 0.06 μM (NH₄)₆Moₐ₇O₄. 24 , 15μM H3BO3, 8μM MnCl2, 0.12μM CuSO4, 0.12μM ZnSO4, 29μM FeCl3, 40.5μM Citric acid, pH 5.5.

[0127] Plants were then transferred to a nutrient solution without NH₄NO₃ for 7 days of nitrogen starvation. At 9:00 AM on day 8, the plants were transferred to a complete nutrient solution containing 1.44 mM NH₄NO₃. Root samples were collected at 0, 3, 6, 9, 12, 15, 20, 30, 60, and 120 minutes, with two biological replicates per time point, each containing roots from two plants.

[0128] ATAC-seq experiments were performed strictly according to standard protocols: 0.2 g of fresh root tissue was minced rapidly in lysis buffer (15 mM Tris-HCl pH 7.5, 20 mM NaCl, 80 mM KCl, 0.5 mM spermidine, 3 mM DTT, 0.2% Triton X-100) on ice, and the single-nuclei suspension was filtered. 100,000 nuclei were sorted using flow cytometry and labeled by incubation with Tn5 transposase at 37°C for 30 minutes. DNA was purified and amplified by PCR to construct a sequencing library.

[0129] For RNA-seq experiments, total RNA was extracted using the OminiPlant RNA kit, and RNA quality was assessed by NanoDrop and agarose gel electrophoresis. Strand-specific sequencing libraries were constructed.

[0130] Both libraries were sequenced using an Illumina HiSeq X10 platform using 150bp paired-end sequencing. ATAC-seq generated an average of 22 million high-quality reads per sample, while RNA-seq generated an average of 64 million reads per sample, both meeting publication standards.

[0131] Bioinformatics analysis revealed 3,311 DEGs and 28,547 DARs in ZS97, and 4,275 DEGs and 23,350 DARs in NIP. Using the eCAAS method, 174 transcription factors differentially expressed in both varieties were analyzed, and a regulatory network encompassing these transcription factors and their target genes was successfully constructed.

[0132] Example 3 Identification and Verification of the Key Transcription OsLBD38

[0133] Based on the above method, the key transcription factor OsLBD38 was identified and verified as follows:

[0134] First, we extracted OsLBD38 expression at each time point and analyzed its association with the accessibility of 150,535 ACRs across the entire genome. Using a linear mixed-effects model, we calculated the regression coefficient β1 and the corresponding P value, identifying 1,579 ACRs that were significantly correlated with OsLBD38 expression (P < 0.05).

[0135] Mapping these ACRs to the nearest genes yielded 1,579 potential target genes. Searching for the LBD family DNA binding motif "GCGGCG" within these ACRs revealed 312 ACRs containing this motif, and the corresponding genes were considered direct target genes.

[0136] Network topology analysis calculated the importance score of OsLBD38. Its degree centrality (DC) was 1579 (number of connected target genes), betweenness centrality (BC) was 0.23 (frequency of occurrence in the shortest path in the network), and closeness centrality (CC) was 0.15 (the inverse of the average distance to other nodes in the network). Based on the formula IS = 0.5 × 1579 + 0.3 × 0.23 + 0.2 × 0.15 = 789.6, it ranked fourth among 174 transcription factors, identifying it as a key transcription factor.

[0137] Experimental validation results showed that, after constructing an OsLBD38 overexpression strain and conducting RNA-seq analysis, 736 of the 1,579 predicted target genes showed significant expression changes in the overexpression strain, a validation rate of 46.5%. Further CUT&Tag experiments revealed that OsLBD38 binding signals were significantly enriched in the promoter regions of predicted target genes, particularly those related to nitrogen metabolism, such as OsNIA1 and OsNiR2.

[0138] Importantly, the eCAAS method not only predicted regulatory relationships but also precisely located regulatory regions. For example, a key enhancer region was identified 3.2 kb upstream of the OsNiR2 gene. Chromatin accessibility in this region showed the highest correlation with OsNiR2 expression (r = 0.78, P < 1e-2). This precise location information provides a clear target for subsequent gene editing.

[0139] Example 4: Development of an Improved Strategy Based on Effect Differentiation

[0140] Using Nipponbare as a reference variety, the method for analyzing differences and formulating improvement strategies for the rice variety Zhenshan 97 based on the above method is as follows:

[0141] The nitrogen transporter gene OsNRT2.4 is an important nitrogen transporter. A deep learning model predicted a 0.3 difference in chromatin accessibility in the promoter region of this gene between Zhenshan 97 and Nipponbare, while the experimentally observed difference was 1.5. The CTE value (CTE) = 0.3 / 1.5 = 0.2 < 0.5, indicating a dominant trans-action. The key transcription factor OsLBD37 was significantly more expressed in Nipponbare than in Zhenshan 97 (1.8-fold difference, P < 0.001), affecting gene expression by regulating chromatin accessibility in the OsNRT2.4 promoter region. Therefore, a strategy for improving the nitrogen transport capacity of Zhenshan 97 should prioritize improving the expression level of OsLBD37. The recommended target is the promoter region of the OsLBD37 gene, rather than direct editing of the OsNRT2.4 gene.

[0142] Similarly, the NAD-dependent malate dehydrogenase gene, OsMDH1, is an important nitrogen assimilation-related gene. A deep learning model predicted a chromatin accessibility difference of 1.2 between Zhenshan 97 and Nipponbare in the promoter region of this gene, compared to an experimentally observed difference of 1.3. The CTE value (CTE) = 1.2 / 1.3 = 0.92 > 0.8, indicating a dominant cis-effect. Sequence alignment revealed a variety-specific polymorphic site in the promoter region of this gene, and the haplotype in Nipponbare could enhance gene expression. Therefore, a suggested improvement strategy is to introduce the haplotype in the OsMDH1 region of Nipponbare into Zhenshan 97 to enhance its nitrogen transport capacity. Furthermore, the promoter region of this gene in Zhenshan 97 could be used as a target for gene editing to create varieties with higher OsMDH1 expression.

[0143] In summary, this invention achieves systematic identification of key transcription factors through the innovative eCAAS method and precise differentiation of cis- and trans-actions through a deep learning model. This method can provide targeted molecular improvement strategies for different gene types, providing important technical support for the molecular design of nitrogen-efficient rice varieties. This technical solution boasts high precision, rapid efficiency, and wide applicability, and holds broad application prospects in the field of modern agricultural molecular breeding.

[0144] Although the above embodiments have been described in detail, they are only a part of the embodiments of the present invention, not all of them. People can also obtain other embodiments based on this embodiment without inventiveness, and these embodiments all fall within the scope of protection of the present invention.

Claims

1. A method for constructing a rice nitrogen response regulatory network based on multi-omics time series data, characterized in that: The following steps are involved: S1. Multi-omics time-series data collection: Seedlings of a single or multiple rice varieties are placed in a nitrogen-deficient nutrient solution for 7 days of nitrogen starvation pretreatment. The seedlings are then transferred to a normal nutrient solution for culture. Root samples are collected at ≥10 time points within 48 hours. The root samples at each time point are divided into two equal parts for ATAC-seq and RNA-seq sequencing to obtain chromatin accessibility and gene expression data for the root samples at each time point. The collection of the two data at all time points is the multi-omics time-series data; S2. Construction of regulatory network and identification of key transcription factors: The expression-chromatin accessibility association study eCAAS method was used to associate the dynamic changes of transcription factor expression with the genome-wide chromatin accessible regions through a linear mixed model. The chromatin accessible regions associated with transcription factors were dynamically associated with the expression of the genes closest to them through a linear mixed model to construct a nitrogen-responsive transcriptional regulatory network, that is, a regulatory network including transcription factors, target chromatin accessible regions and target genes.

2. The method according to claim 1, characterized in that In step S1, after the ATAC-seq and RNA-seq raw sequencing data undergo data quality control, the ATAC-seq data uses the MACS2 tool to identify chromatin accessible regions, and the RNA-seq data uses the salmon tool to quantify gene expression.

3. The method according to claim 1, characterized in that In step S2, the eCAAS method includes two levels of correlation analysis: a. Transcription factor-chromatin accessible region association analysis: A linear mixed-effects model was established, with transcription factor expression as the independent variable and the degree of chromatin accessible region openness as the dependent variable: A ijk =β0+β1Y ik +g jk +e ijk Among them, A ijk represents the accessibility level of chromatin accessible region j in sample i and time point k, β0 represents the fixed effect intercept term, β1 represents the fixed effect slope coefficient, and Y ik represents the expression level of transcription factor in sample i and time point k, γ jk is a random effect term, which obeys the multivariate normal distribution N(0,∑ RNA ), ε ijk represents the residual error term; The variance-covariance matrix ∑ of the random effect term RNA Sample correlation calculation by RNA-seq data: Where n represents the total number of samples measured in the RNA-seq data, X g Represents the expression vector of gene g in all samples, is the average expression level of gene g in all samples, and G is the total number of genes; b. Regulatory region-target gene association analysis: A linear mixed-effects model was established, with the degree of openness of the chromatin accessible region obtained in step a as the independent variable and the expression level of the gene closest to the chromatin accessible region as the dependent variable: E ijk =α0+α1A jk +d ik +∈ ijk Among them, E ijk represents the expression level of target gene i in sample k and time point j, α0 is the fixed effect intercept term, α1 is the fixed effect slope coefficient, and A jk represents the chromatin accessibility level of the regulatory region at sample k and time point j, δ ik is a random effect term, which obeys the multivariate normal distribution N(0,∑ ATAC ), ε ijk represents the residual error term; The variance-covariance matrix ∑ of the random effect term ATAC Sample correlation calculation by ATAC-seq data: Where n represents the total number of samples measured in the ATAC-seq data, and Z r represents the signal intensity vector of the rth chromatin region in all samples, The average signal intensity of the rth chromatin accessible region, R is the total number of regulatory regions.

4. An application of the regulatory network constructed by the method of claim 1 in identifying key transcription factors and the location of their regulatory regions, analyzing variety differences, and formulating improvement strategies.

5. A method for identifying key transcription factors and their regulatory region locations, characterized by: The method comprises the following steps: using a nitrogen-responsive transcriptional regulatory network constructed by the method of claim 1, and using a network topology analysis method to calculate the regulatory importance score of each transcription factor; sorting the importance scores, identifying the top 10% of the transcription factors as key transcription factors, and precisely locating the position of the chromatin regulatory region corresponding to each regulatory relationship.

6. The method according to claim 5, characterized in that: The formula for calculating the regulatory importance score of each transcription factor is as follows: Among them, IS TFk Transcription factor TF k Importance score, DC TFk is degree centrality, BC TFk is the betweenness centrality, CC TFk is the closeness centrality, α, β, and γ are weight coefficients, and their values ​​are 0.5, 0.3, and 0.2 respectively.

7. A method for analyzing rice variety differences and identifying breeding targets, characterized by: The following steps are involved: 1) obtaining the genome sequence of the target variety material to be analyzed and obtaining multi-omics time series data of the target variety material to be analyzed based on the method of step S1 of claim 1; 2) Determine candidate regulatory regions using the regulatory network constructed using the method of claim 1, predict the chromatin accessibility of the candidate regulatory regions of the target variety to be analyzed using a deep learning model, and perform differential analysis with the provided ATAC-seq data; based on the degree of difference between the predicted chromatin accessibility of the target variety and the chromatin accessibility in the provided ATAC-seq data; distinguish the contributions of cis-acting effects and trans-acting effects, determine breeding targets based on the effect type, and provide variety improvement strategies: If the candidate regulatory region is identified as a trans-acting effector, a list of candidate upstream transcription factors recommended for variety improvement is output; Alternatively, if it is a cis-acting effect, the location of this regulatory region can be directly output for variety improvement.

8. The method according to claim 7, wherein: In step 2), a deep learning model is constructed based on an improved Basenji framework, including a convolutional neural network module and a multi-task learning framework, to predict chromatin accessibility at different time points. The loss function of the deep learning model is: Among them, L accessibility Represents the loss function value of chromatin accessibility prediction, which measures the difference between the predicted result and the true value. N is the total number of samples. is the predicted accessibility value of the i-th chromatin region, is the true accessibility value of the i-th chromatin accessible region; Identify regulatory regions associated with any transcription factor in the target variety, predict the chromatin accessibility of this regulatory region based on the above deep learning model, and compare the predicted results with the chromatin accessibility results of the regulatory region in the ATAC-seq data provided by the target variety to distinguish cis- and trans-effects. The quantitative indicators are: Among them, CTE genei is the cis-trans effect ratio of gene i, A pred,target The chromatin accessibility degree predicted by the model for the target species, A exp,target is the degree of chromatin accessibility obtained from the ATAC-seq data of the target species; CTE-based genei The specific value of is as follows: When CTE genei When ≥0.8, it was determined to be cis-effect dominant; Or, when CTE genei When <0.5, it was determined to be dominated by the trans effect; Or, when 0.5≤CTE genei When <0.8, it was determined to be a mixed effect.

9. The method according to claim 8, characterized in that: In step 2), breeding targets are determined based on the effect type determination results, and variety improvement strategies are formulated: For regulatory regions dominated by trans effects, output a list of upstream transcription factors associated with the regulatory region, and the expression difference of each transcription factor compared with the corresponding transcription factor in the reference variety; output in descending order based on the regulatory strength of the transcription factors with the regulatory region dominated by the trans effect; For regulatory regions dominated by cis effects, accurate regulatory region coordinate information is output, and the haplotypes of the target varieties can be directly used for breeding or as gene editing targets; For regulatory regions dominated by trans effects, output the transcription factors associated with the regulatory regions. The recommended transcription factors are ranked based on the regulatory strength score: RS target =w1·|FC|+w2·(-log 10 (P value ))+w3·CS Among them, RS target is the transcription factor recommendation score, FC is the maximum expression change fold of the transcription factor in the time series data, Pvalue is the eCAAS significance P value of the transcription factor and the selected regulatory region, CS is the interaction score of the homologous "upstream transcription factor-target" gene pair in other species such as Arabidopsis, w1, w2, and w3 are weight coefficients, which are 0.4, 0.4, and 0.2 respectively; Alternatively, if there is no other species information, the CS weight is set to 0.

10. A rice nitrogen response regulatory network analysis and breeding target identification system based on multi-omics data for implementing the method according to any one of claims 1, 5, and 7, characterized in that: include: Data preprocessing module, used to preprocess multi-omics time series data, including data quality control, format conversion and standardization functions; A network construction module, used to construct transcriptional regulatory networks and identify key transcription factors, enabling efficient parallel computing of the eCAAS algorithm; Effect differentiation module, used to analyze regulatory differences between varieties and distinguish effect types, integrating deep learning prediction models; The target recommendation module is used to formulate improvement strategies and recommend molecular targets based on effect types, providing personalized improvement plans.

Images