s2fieldsofstudy: Biology Chemistry

text:
Neuropeptides, New Ligands of SARS-CoV-2 Nucleoprotein, a Potential Link between Replication, Inflammation and Neurotransmission

This work identifies new ligands of the nucleoprotein N of SARS-CoV-2 by in silico screening, which used a new model of N, built from an Alphafold model refined by molecular dynamic simulations. The ligands were neuropeptides, such as substance P (1-7) and enkephalin, bound at a large site of the C-terminal or associated with the N-terminal β−sheet. The BA4 and BA5 Omicron variants of N also exhibited a large site as in wt N, and an increased flexibility of the BA5 variant, enabling substance P binding. The binding sites of some ligands deduced from modeling in wt N were assessed by mutation studies in surface plasmon resonance experiments. Dynamic light scattering showed that the ligands impeded RNA binding to N, which likely inhibited replication. We suggest that the physiological role of these neuropeptides in neurotransmission, pain and vasodilation for cholecystokinin and substance P could be altered by binding to N. We speculate that N may link between viral replication and multiple pathways leading to long COVID-19 symptoms. Therefore, N may constitute a “danger hub” that needs to be inhibited, even at high cost for the host. Antivirals targeted to N may therefore reduce the risk of brain fog and stroke, and improve patients’ health.

Introduction
The nucleoprotein N is one of the four structural proteins of SARS-CoV-2 virus [1][2][3][4]. N is present in large number of copies. As such, N is one of the major targets for antibody development and has been widely used for COVID-19 detection in the present pandemic outbreak. N binds the long viral RNA genome and is associated in the ribonucleoprotein (RNP) complex required for viral replication [5,6]. N is also involved in the formation of new virions through its interactions with the membrane protein M. Besides these functions, N acts as a mediator of inflammation. N represses the host antiviral response (as RNA interference and RIG-I mediated interferon) [7,8]. N targets the stress granule protein G3BP1, an essential antiviral protein known to induce innate immune response [9]. Importantly, N is associated with long-term SARS-CoV-2-specific immune and inflammatory responses, since the frequency of N-specific interferon-γ-producing CD8 + T cells decline more rapidly in long-hauler COVID-19 patients [10]. N also activates endothelial cells dysfunction, leading to vasculopathy and coagulopathy observed in some COVID-19 patients [11]. N is thus an important target for the development of antivirals [12]. It is recognized that SARS-CoV-2 induces an early host inflammatory response that activates a cyclooxygenase-2 (COX-2) inflammatory cascade associated with NF-κB activation [13,14]. Indeed, in a mouse model, N protein promoted the expression of proinflammatory cytokines and triggered lung injury via NF-κB activation [11,15]. The use of anti-inflammatory drugs, especially non-steroidal anti-inflammatory drugs (NSAIDs), in the initial outpatient stage of COVID-19 appears to be a valuable therapeutic strategy [16]. We identified, using a structure-based approach a cyclooxygenase (COX) inhibitor, naproxen, as an antiviral against SARS-CoV-2 and Influenza A virus that combines anti-inflammatory properties [17][18][19][20]. We showed that naproxen bound to SARS-CoV-2 N-terminal domain (NTD) in vitro and reduced the viral load of infected cells. Naproxen protected lung cells against viral injury in a model of lung epithelium, in contrast with celecoxib or paracetamol [17]. NTD was also shown to bind AMP, NADPH and single-stranded RNA [3,12].
In this paper, we extend our previous work and identify by in silico screening new ligands of the full-length N protein (FL). We report that the nucleocapsid can sequester in a large cavity of its C-terminal various host neuropeptides involved in neurotransmission, vasodilation, inflammation and ligands perturbing cell metabolism. This capture, although likely decreasing viral replication, may contribute to perturbations in brain function and metabolism. We discuss these hypotheses in light of the recent literature. Our study sheds light on the need to block N by antivirals which may overcome some long COVID-19 symptoms [21,22].

Generation of Models of Full-Length Nucleoprotein N
Recent structures of both N amino-terminal domain (NTD) and C-terminal domain (CTD) have been solved by X-ray crystallography, Nuclear Magnetic Resonance (NMR) and electron microscopy [2,3,[23][24][25]. N also bears dynamic, disordered domains: an SR-rich domain between the NTD and CTD, as well as N-and C-terminal tails ( Figure 1A). The isolated NTD formed a monomer which bound RNA and other single-stranded nucleic acid fragments [2].
The N NTD structure presented a central antiparallel sheet of 3-5 β-strands ("the hand palm") with a characteristic extrusion loop ("basic finger") that can close on RNA upon binding [2,3]. The CTD also bound RNA and forms oligomers, mostly dimers [2]. The N protein also underwent liquid-liquid phase separation when mixed with RNA [1]. The full length (FL) N was reported to be a tetramer by size exclusion chromatography coupled to light scattering [2], but may also form a mixture of monomer and dimer [26]. However, the structure of the full-length N protein of SARS-CoV-2 is not yet available, partly because this recombinant protein is resistant to RNase treatments [26].
We present here a model of the full-length N protein initially generated by Alphafold and further optimized by MD simulations as schematized in Figure 1B, and described in the experimental section [27]. Figure 1B also shows the comparison of our model with the NMR structure of the NTD (PDB 7ACT [2]) shown in brown and the X-ray structure of the CTD (PDB 6WZO [3]) shown in blue, with root mean square deviations (RMSD) of 4.0 and 4.7 Å, respectively (see also Supplementary Figure S1 for the superimposition of the structures and the RMSD of two dynamic trajectories). Figure 1C shows that the protein presented a large cavity at its C-terminal domain, highlighted by a star, which could bind large ligands.
While the NTD presented a number of suitable sites for ligands binding, these sites were small compared to the large cavity found at the C-terminal (Figure 2A), also detected in the isolated CTD X-ray structure [4,28]. Moreover, as this model of N was obtained based on the initial Wuhan strain of SARS-CoV-2 (wt FL), we also built the models of two Omicron variants of N to make sure that the detected cavity in wt FL would also bind ligands in more recent variants of N. Of the two Omicron variants of N, BA4 carried five mutations-P13L, P151S, R203K, G204R, S413R-and a deletion-31-33del-and BA5 had the same modifications but no mutation in P151 [29]. Figure 2A shows modifications of the cavities in the two variants, which only differed by the mutation of one residue. The latter cavity in BA5 was shifted apart from the β-sheet. The RMSD of the main chain atoms with respect to the wt ones were 1.57 Å and 1.28 Å in BA4 and BA5, respectively. large ligands.
While the NTD presented a number of suitable sites for ligands binding, these sites were small compared to the large cavity found at the C-terminal (Figure 2A), also detected in the isolated CTD X-ray structure [4,28]. Moreover, as this model of N was obtained based on the initial Wuhan strain of SARS-CoV-2 (wt FL), we also built the models of two Omicron variants of N to make sure that the detected cavity in wt FL would also bind ligands in more recent variants of N. Of the two Omicron variants of N, BA4 carried five mutations-P13L, P151S, R203K, G204R, S413R-and a deletion-31-33del-and BA5 had the same modifications but no mutation in P151 [29]. Figure 2A shows modifications of the cavities in the two variants, which only differed by the mutation of one residue. The latter cavity in BA5 was shifted apart from the -sheet. The RMSD of the main chain atoms with respect to the wt ones were 1.57 Å and 1.28 Å in BA4 and BA5, respectively.  [3] just below are shown in brown, and the CTD and the X-ray structure of the CTD (PDB 6WZO) just below are shown in blue, the SR-rich motif is presented in green, the other linkers in grey; the RMSD of the model compared to 7ACT is 4.7 Å and with one unit of the dimer of 6VYO [30] is 6.2 Å, the RMSD of the model with 6WZO [2] and 7DE1 [24] are 4.0 and 5.1 Å, respectively; the superimposition of the model and the experimental structures is shown in Supplementary Figure S1; (C) FL is shown as a surface with each domain colored with the same color code as depicted in A; note the large, mainly hydrophobic cavity, highlighted by the white star, at the C-terminal; see also   [3] just below are shown in brown, and the CTD and the X-ray structure of the CTD (PDB 6WZO) just below are shown in blue, the SR-rich motif is presented in green, the other linkers in grey; the RMSD of the model compared to 7ACT is 4.7 Å and with one unit of the dimer of 6VYO [30] is 6.2 Å, the RMSD of the model with 6WZO [2] and 7DE1 [24] are 4.0 and 5.1 Å, respectively; the superimposition of the model and the experimental structures is shown in Supplementary Figure S1; (C) FL is shown as a surface with each domain colored with the same color code as depicted in A; note the large, mainly hydrophobic cavity, highlighted by the white star, at the C-terminal; see also   Figure S1 for a 2D representation of the RMSF). The position of the deletion of the 31-33 residues is highlighted, the other mutations are shown in CPK. The RMSF follows a color code from blue (less flexible) to red (most flexible), note, in wt, BA4 and BA5 variants, the modifications of the fluctuations of the helix close to the R203K, S204R mutations located in close vicinity to the basic finger. These data highlight short-and long-range effects of the mutations on the protein flexibility, with a large enhancement in BA5 only missing the P151S mutation as compared to BA4 that may increase replication in the former variant.
Associated with these changes in the cavity, we measured the protein flexibility of the three proteins, estimated by the root mean square fluctuations (RMSF). The model of wt N FL suggested that the -helix close to the basic finger in the NTD was quite mobile (when the protein was not engaged in interaction with ligands), with flexibility of some linkers as part of the serine-rich domain and the C-tail (see also Supplementary Figure  S1).
The R203K, G204R mutations in both variants tilted the edge of the basic finger that binds RNA towards the linker that carries these two positively charged mutated residues. Interestingly, the flexibility of the-helix edge(R209) in close proximity to the R203K, G204R mutations was enhanced in the BA5 variant as compared to BA4, which may help to grasp the RNA more quickly and/or more efficiently. The data suggested subtle, longrange effects of the mutations (here the mutation status of residue 151 influencing the mobility of R209) that may modify replication. The N variant carrying the R203K, G204R mutations increased the infectivity, fitness and virulence of SARS-CoV-2 [30]. Altogether, the variants presented a modified cavity as compared to wt N, which nevertheless remained quite large and very flexible.

Characterization of Recombinant N, N-NTD and Its Mutants and N Full Length N-FL
The recombinant NTD protein was expressed as previously described [17]. NTD purified as a monomer deduced from SEC-MALS analysis [4]. The full-length protein (FL) was expressed and purified according to a similar protocol than that for NTD. However, the purified FL found in the soluble fraction remained contaminated by bacterial nucleic acids as attested by an absorbance ratio 260 nm/280 nm ranging between 0.9 to 1.3, despite the use of benzonase and RNase [26]. Instead, we used the insoluble fraction that was denatured with 6 M urea and then renatured in Tris-NaCl buffer. In these conditions, the absorbance ratio 260 nm/280 nm decreased to 0.47, consistent with no significant contam- Associated with these changes in the cavity, we measured the protein flexibility of the three proteins, estimated by the root mean square fluctuations (RMSF). The model of wt N FL suggested that the α-helix close to the basic finger in the NTD was quite mobile (when the protein was not engaged in interaction with ligands), with flexibility of some linkers as part of the serine-rich domain and the C-tail (see also Supplementary Figure S1).
The R203K, G204R mutations in both variants tilted the edge of the basic finger that binds RNA towards the linker that carries these two positively charged mutated residues. Interestingly, the flexibility of the α-helix edge (R209) in close proximity to the R203K, G204R mutations was enhanced in the BA5 variant as compared to BA4, which may help to grasp the RNA more quickly and/or more efficiently. The data suggested subtle, long-range effects of the mutations (here the mutation status of residue 151 influencing the mobility of R209) that may modify replication. The N variant carrying the R203K, G204R mutations increased the infectivity, fitness and virulence of SARS-CoV-2 [30]. Altogether, the variants presented a modified cavity as compared to wt N, which nevertheless remained quite large and very flexible.

Characterization of Recombinant N, N-NTD and Its Mutants and N Full Length N-FL
The recombinant NTD protein was expressed as previously described [17]. NTD purified as a monomer deduced from SEC-MALS analysis [4]. The full-length protein (FL) was expressed and purified according to a similar protocol than that for NTD. However, the purified FL found in the soluble fraction remained contaminated by bacterial nucleic acids as attested by an absorbance ratio 260 nm/280 nm ranging between 0.9 to 1.3, despite the use of benzonase and RNase [26]. Instead, we used the insoluble fraction that was denatured with 6 M urea and then renatured in Tris-NaCl buffer. In these conditions, the absorbance ratio 260 nm/280 nm decreased to 0.47, consistent with no significant contamination by bacterial nucleic acids. The FL was then tested by DLS for its ability to bind RNA (TAR-polyA). The size of the native protein with bacterial nucleic acids contamination was 14.5 ± 1.0 nm in volume (data not shown). As shown in Figure 1D, the size of the protein depleted from nucleic acids contaminants was 12.0 ± 0.4 nm in volume (14.0 ± 0.3 nm in intensity) that increased to 15.1 ± 0.3 nm (20.3 ± 0.3 nm in intensity) upon addition of RNA ( Figure 1D), showing that the protein is functional. By comparison, the size of the NTD monomeric protein was 6.0 ± 0.3 nm and increased to 15.7 ± 0.5 nm in intensity in the presence of nucleic acids. The size of NTD could not be fitted by a spherical model to obtain its value in volume. To further characterize FL protein, the apparent melting temperature (Tm) of the protein alone was determined by DLS as a function of temperature. FL presented a first transition at 43 ± 1 • C followed by a denaturation at about 54 • C, (Supplementary Figure S2), consistent with the hypothesis that the protein could be (at least partly) dimeric at 25 • C, with a transition to random coil upon heating at 54 • C, a transition also determined by circular dichroism [31]. Preliminary CD experiments showed that the spectrum of FL is dependent upon the protein concentration as previously reported.

Virtual Drug Screening on NTD and FL
We performed virtual screening to identify ligands binding to the large C-terminal cavity of the full-length protein ( Figure 2A) using the data base of the Sigma-Aldrich catalog. We also screened for ligands binding to NTD. We first used the Libdock software of Discovery studio. Then, docking of the interesting hits was further repeated using CDocker. Finally, MD simulations were performed on the best hits. Tables 1 and 2 summarize the mean-full hits both in terms of docking and biological function.

Ligands involved in metabolism
Dihydrofolate (DHF) and tetrahydrofolate (THF) are, respectively, the substrate and product of the enzyme dihydrofolate reductase (DHFR).
DHF binding to NTD: DHF interacted with NTD in its monomeric and dimeric forms ( Figure 3A,B, respectively). It was anchored by arginine R149 (or R107 in NTD dimer) forming a -cation complex with one of the aromatic rings of DHF, a  complex with H145 and further stabilized by a number of H bonds or electrostatic interactions with G44, N77, N150 (NTD monomer) and R92, S105 (NTD dimer). DHF binding on one face of the protein induced a compaction of the basic finger binding RNA on the other face, suggesting long range or allosteric effects. This compaction involved, in particular, R88 stabilized by salt bridges with D98 and E118. As R88, together with R92 and R107, have been all

Ligands involved in metabolism
Dihydrofolate (DHF) and tetrahydrofolate (THF) are, respectively, the substrate and product of the enzyme dihydrofolate reductase (DHFR).
DHF binding to NTD: DHF interacted with NTD in its monomeric and dimeric forms ( Figure 3A,B, respectively). It was anchored by arginine R149 (or R107 in NTD dimer) forming a -cation complex with one of the aromatic rings of DHF, a  complex with H145 and further stabilized by a number of H bonds or electrostatic interactions with G44, N77, N150 (NTD monomer) and R92, S105 (NTD dimer). DHF binding on one face of the protein induced a compaction of the basic finger binding RNA on the other face, suggesting long range or allosteric effects. This compaction involved, in particular, R88 stabilized by salt bridges with D98 and E118. As R88, together with R92 and R107, have been all

Ligands involved in metabolism
Dihydrofolate (DHF) and tetrahydrofolate (THF) are, respectively, the substrate and product of the enzyme dihydrofolate reductase (DHFR).
DHF binding to NTD: DHF interacted with NTD in its monomeric and dimeric forms ( Figure 3A,B, respectively). It was anchored by arginine R149 (or R107 in NTD dimer) forming a -cation complex with one of the aromatic rings of DHF, a  complex with H145 and further stabilized by a number of H bonds or electrostatic interactions with G44, N77, N150 (NTD monomer) and R92, S105 (NTD dimer). DHF binding on one face of the protein induced a compaction of the basic finger binding RNA on the other face, suggesting long range or allosteric effects. This compaction involved, in particular, R88 stabilized by salt bridges with D98 and E118. As R88, together with R92 and R107, have been all

Ligands involved in metabolism
Dihydrofolate (DHF) and tetrahydrofolate (THF) are, respectively, the substrate and product of the enzyme dihydrofolate reductase (DHFR).
DHF binding to NTD: DHF interacted with NTD in its monomeric and dimeric forms ( Figure 3A,B, respectively). It was anchored by arginine R149 (or R107 in NTD dimer) forming a -cation complex with one of the aromatic rings of DHF, a  complex with H145 and further stabilized by a number of H bonds or electrostatic interactions with G44, N77, N150 (NTD monomer) and R92, S105 (NTD dimer). DHF binding on one face of the protein induced a compaction of the basic finger binding RNA on the other face, suggesting long range or allosteric effects. This compaction involved, in particular, R88 stabilized by salt bridges with D98 and E118. As R88, together with R92 and R107, have been all

Ligands involved in metabolism
Dihydrofolate (DHF) and tetrahydrofolate (THF) are, respectively, the substrate and product of the enzyme dihydrofolate reductase (DHFR).
DHF binding to NTD: DHF interacted with NTD in its monomeric and dimeric forms ( Figure 3A,B, respectively). It was anchored by arginine R149 (or R107 in NTD dimer) forming a -cation complex with one of the aromatic rings of DHF, a  complex with H145 and further stabilized by a number of H bonds or electrostatic interactions with G44, N77, N150 (NTD monomer) and R92, S105 (NTD dimer). DHF binding on one face of the protein induced a compaction of the basic finger binding RNA on the other face, suggesting long range or allosteric effects. This compaction involved, in particular, R88 stabilized by salt bridges with D98 and E118. As R88, together with R92 and R107, have been all

Ligands involved in metabolism
Dihydrofolate (DHF) and tetrahydrofolate (THF) are, respectively, the substrate and product of the enzyme dihydrofolate reductase (DHFR).
DHF binding to NTD: DHF interacted with NTD in its monomeric and dimeric forms ( Figure 3A,B, respectively). It was anchored by arginine R149 (or R107 in NTD dimer) forming a -cation complex with one of the aromatic rings of DHF, a  complex with H145 and further stabilized by a number of H bonds or electrostatic interactions with G44, N77, N150 (NTD monomer) and R92, S105 (NTD dimer). DHF binding on one face of the protein induced a compaction of the basic finger binding RNA on the other face, suggesting long range or allosteric effects. This compaction involved, in particular, R88 stabilized by salt bridges with D98 and E118. As R88, together with R92 and R107, have been all

Ligands involved in metabolism
Dihydrofolate (DHF) and tetrahydrofolate (THF) are, respectively, the substrate and product of the enzyme dihydrofolate reductase (DHFR).
DHF binding to NTD: DHF interacted with NTD in its monomeric and dimeric forms ( Figure 3A,B, respectively). It was anchored by arginine R149 (or R107 in NTD dimer) forming a -cation complex with one of the aromatic rings of DHF, a  complex with H145 and further stabilized by a number of H bonds or electrostatic interactions with G44, N77, N150 (NTD monomer) and R92, S105 (NTD dimer). DHF binding on one face of the protein induced a compaction of the basic finger binding RNA on the other face, suggesting long range or allosteric effects. This compaction involved, in particular, R88 stabilized by salt bridges with D98 and E118. As R88, together with R92 and R107, have been all

Ligands involved in metabolism
Dihydrofolate (DHF) and tetrahydrofolate (THF) are, respectively, the substrate and product of the enzyme dihydrofolate reductase (DHFR).
DHF binding to NTD: DHF interacted with NTD in its monomeric and dimeric forms ( Figure 3A,B, respectively). It was anchored by arginine R149 (or R107 in NTD dimer) forming a -cation complex with one of the aromatic rings of DHF, a  complex with H145 and further stabilized by a number of H bonds or electrostatic interactions with G44, N77, N150 (NTD monomer) and R92, S105 (NTD dimer). DHF binding on one face of the protein induced a compaction of the basic finger binding RNA on the other face, suggesting long range or allosteric effects. This compaction involved, in particular, R88 stabilized •

Ligands involved in metabolism
Dihydrofolate (DHF) and tetrahydrofolate (THF) are, respectively, the substrate and product of the enzyme dihydrofolate reductase (DHFR).
DHF binding to NTD: DHF interacted with NTD in its monomeric and dimeric forms ( Figure 3A,B, respectively). It was anchored by arginine R149 (or R107 in NTD dimer) forming a π-cation complex with one of the aromatic rings of DHF, a π−π complex with H145 and further stabilized by a number of H bonds or electrostatic interactions with G44, N77, N150 (NTD monomer) and R92, S105 (NTD dimer). DHF binding on one face of the protein induced a compaction of the basic finger binding RNA on the other face, suggesting long range or allosteric effects. This compaction involved, in particular, R88 stabilized by salt bridges with D98 and E118. As R88, together with R92 and R107, have been all involved in RNA binding to NTD by NMR and further confirmed by mutation studies [2], it is likely that DHF binding strongly reduced RNA binding to N NTD.  Table 2) Although DHF formed stronger polar interactions in FL than DHF in NTD, it did not form multiple hydrophobic and π-π interactions as in NTD.

DHF in FL (
AICAR (Table 1) AICAR is an agonist of AMPK; it bound through electrostatic interactions at a site at the basic finger where RNA binds, involving K61, R89, R93, R95, K102 and R107 via its charged and polar groups.
Lauroyl-CoA in FL (Tables 1 and 2) The long aliphatic chain of lauroyl CoA was recognized by hydrophobic interactions with a series of five leucines and isoleucines residues and H356, F286 while its polar moiety forms H-bonds with S187 and S188 and K261 and R259 residues.

•
Ligands involved in the COX-prostaglandins pathway: The NSAIDs as naproxen bound the NTD [18]. They also can bind the C-terminal cavity although their size was not fitted to the large cavity. Table 2 describes the binding of acetamine at the C-terminal. Additionally, prostaglandins A2 and F2 and other eicosanoids bound FL with low affinity, further indicating the importance of this pathway in the host response to viral infection.

•
Neuropeptides: Substance P (1-7) in FL bound at a similar site than observed with other neuropeptides, enkephalin, cholescystokinin and lauroyl coenzyme A (Figures 4-6 and Table 2). Substance P (1-7) formed marked hydrophobic interactions with F286, Y298, I304 and Y333 as observed for enkephalin, and also binds S188 and R189 (Figure 4). This binding also resulted in a compaction of the RNA-binding finger, stabilized through interactions involving R95, D98, K100 and T205 although no direct interaction with the NTD was observed, including longrange effects between different sub-domains of the protein. Other fragments of substance P (4-11) also bound FL at a similar site and interaction energy than the P (1-7) fragment. shown in green that formed multiple stacking and hydrophobic interactions, further stabilized by electrostatic and polar interactions; (C) consequently, the effect of the NTD mutations on the SPR signal was mainly observed with Y111A and W52A, with little effect seen on the other mutations. Substance P (1-7) also bound the N FL protein as measured by SPR.
To further test whether substance P (1-7) may bind the BA4 and BA5 variants of N, we introduced the relevant mutations and deletion in the wt FL-substance P (1-7) model. Figure 5 shows that substance P (1-7) tightly bound to both wt and BA5 variant of N with similar interactions. Figure 5C also highlights the decrease of flexibility of N due to peptide binding in the structural element involved in RNA binding, in agreement with the competition between the ligand and RNA for binding to N (see below Table 3).  10.3 ± 0.6 (50%) 8.0 ± 1.0 (76%) 9.5 ± 0.6 (98%) 7.2 ± 1.1 (99%) FL 2 µM TPPS 2 7.6 ± 0.6 (25%) 7.3 ± 0.6 (99%) FL 2 µM   Enkephalin in NTD ( Figure 6): Enkephalin bound at a close-by site from that of DH via electrostatic interactions with R88 and R92. In addition, enkephalin extended toward the central sheet forming hydrophobic and polar interactions with Y111 and Y112. Thi binding induced some structuration of the linker carrying H145 and R149 in a short he lix.  This binding induced some structuration of the linker carrying H145 and R149 in a short α−helix.
Enkephalin in FL ( Figure 6A,C): Enkephalin interacted with FL through hydrophobic and π-π interactions with F274, F286, F298 and Y333 surrounding the ligand, stabilized by K261 and S318. Even though enkephalin bound FL through multiple hydrophobic interactions in the C-terminal cavity, it "closed" the basic finger on the N-terminal by H-bonds involving R95 with G97, G99 and T205, again suggesting long-range effects within the N structure. This hypothesis is confirmed by competition studies of enkephalin with RNA binding to N monitored by DLS (Paragraph 2-5).
Although we did not test directly the binding of enkephalin in the BA4 and BA5 variants of N FL, we anticipate that enkephalin would accommodate even more easily the variants than substance P as it can adopt an extended or folded conformation.
Cholescystokinin similarly bound FL through extensive hydrophobic interactions in the C-terminal cavity and also extended towards the N-terminal β-sheet via R107 and Y109. This binding mode was similar to that observed for substance P ( Table 2).

• Hemin in FL:
Hemin bound FL through hydrophobic interactions with V270, L291, I304, I337 and π-π stacking with F274, in agreement with [32]. The stability of the complex was lower than that observed with the neuro/vasopeptides ( Table 2).

SPR Study of DHF and Substance P (1-7) Interactions with NTD WT and Mutants, and FL
To comparatively assess the binding of DHF and substance P (1-7) to N, we generated mutants of the NTD recombinant protein by targeting some of the residues involved in ligand binding as detailed in experimental Section 4.1.3. The effect of such mutations on DHF binding to NTD are shown in Figure 3. Mutations H145A and R149A strongly reduced the signal of DHF as compared to DHF binding to NTD (monomer) wt as expected from the modeling. In addition to these mutations, we also observed strong effects of the mutations R92A, Y111A and R88A. These residues corresponded to those involved in DHF binding in NTD dimer. As compared to DHF signal in NTD, binding of DHF to the full-length protein was smaller. The relative signal ratio (DHF in NTD)/(DHF in FL) was dependent on the buffer used. It is possible that the oligomeric status of FL was different in HEPES buffer supplemented with 0.05% Tween, compared to Tris buffer. FL was reported to be either tetrameric or a mixture of monomer and dimer [2,26]. Monomer or dimer form of N FL, in particular in the C-terminal provided more accessible site(s) for ligand binding than the tetrameric form did. DHF K D for binding to NTD was about 1 µM, deduced from the variation of the signal as a function of DHF concentration in the range 0.3-10.0 µM.
Substance P (1-7) binding to NTD and FL was also observed as shown in Figure 4C. This large ligand bound NTD and exhibited relatively small effects of the mutations with notable exceptions of W52A and Y111A. This suggested that substance P (1-7) bound on the β-sheet at the conserved sequence of the five aromatic residues W 108 YFYY 112 . Substance P (1-7) binding to N therefore seemed driven by hydrophobic interactions, consistent with the sequence of the peptide. It was also in line with both the binding of the peptide to FL seen by SPR and the multiple hydrophobic contacts of this peptide in the C-terminal cavity of the FL protein suggested by modeling.

FL Interactions with the Ligands by DLS; Competition with Single-Stranded DNA or RNA Binding
To further address whether the ligands may modify the oligomeric status of the FL protein, we performed dynamic light scattering experiments. FL (2 µM) presented a main peak at 12.0 ± 0.6 nm in volume (14 nm in intensity); the size of the complex decreased to 10.3 and 7.8 ± 0.3 nm upon addition of 2 µM and 6 µM enkephalin, respectively ( Figure 6D). As shown in Table 3, the same trend was observed upon addition of naproxen, indomethacin, substance P (1-7) or a porphyrin ZnTPPS with a larger decrease observed with substance P (1-7), ZnTPPS and indomethacin. This suggested that ligand binding decreased the oligomerization state of the protein. We then tested whether ligand binding to FL could be competitive with nucleic acids binding. Table 3 indeed shows that the substance P (1-7) or enkephalin ( Figure 6E) are competitive with nucleic acids binding to FL as the size of the RNA-FL or DNA-FL complex about 14 nm observed without ligand always decreased by addition of the ligands.

FL Interactions with ZnTPPS
In our in silico screening, we identified hemin as a potential ligand binding to FL C-terminal site as recently proposed; hemin binding to N was shown to reduce viral replication [32]. Figure 7A shows that the binding site of hemin colored in red. Here, we chose a water-soluble, negatively charged porphyrin with properties of a photosensitizer, able to produce singlet oxygen and ROS in the perspective of a potential antiviral which could sensitize infected tissues containing N. ZnTPPS has a Soret band at 420 nm and two Q bands, being a metalated porphyrin. Upon addition of FL, the Soret band decreased and red-shifted to 431 nm with the presence of an isosbestic point at 426 nm, while the first Q band shifted from 557 to 561 nm. The data were repeated at various concentrations of ZnTPPS, yielding a K D = 0.40 ± 0.10 µM ( Figure 7B,C).

Discussion
In this work, we identified new ligands of N in the perspective of drug repurposing, with perspectives in basic understanding of N functions.
The model of FL, the full-length nucleoprotein of the wt sequence we built, was generated by Alphafold, with subsequent refinements using molecular dynamics simulations. The model clearly showed that the protein is very flexible, nevertheless the simulations successfully produced a folded model, presenting a large cavity at the C-terminal of the protein able to bind ligands with a broad distribution of molecular weights, similar to the cavity found in the X-ray structure of the isolated C-terminal protein. By introducing the mutations and deletion found in the Omicron BA4 and BA5 variants of N, we also detected a large cavity at the C-terminal, although modified as compared to the wt one. It is interesting to note the increased flexibility of the BA5 variant as compared to that of the BA4 variant, in particular in a linker carrying the 203-204 mutations just opposite the basic finger that bound RNA. This suggested a better fitness to replication in the BA5 variant because of a better/faster adaptation to the viral RNA as compared to that of the BA4 variant [30].
The ligands we identified by in silico screening are likely to decrease viral replication as: (i) the ligands competed with RNA binding; (ii) the ligands reduced FL (oligomer) size. Since N is known to oligomerize upon binding to RNA, both mechanisms are expected to interfere with viral replication. Four ligands we identified in our screening were indeed recently shown to decrease viral replication: naproxen, indomethacin, DHF and hemin [17,18,34,35]. The DLS data were consistent with a competition of DHF, naproxen, indomethacin, enkephalin, substance P (1-7) and ZnTPPS with RNA binding to FL.
Besides its roles in viral transcription and replication, N is involved in immunity via its control in host interferon release, in cytoskeleton rearrangement. The literature further suggests additional roles of N that our data may further highlight, as detailed below.

N as a Mediator of Inflammation and Its Effect in Long-Haul COVID-19
A sustained inflammation that extended well beyond clearance of the primary infection was observed in long-hauler patients and also found in a hamster model of SARS-CoV-2 infection [15,36]. The expression of prostaglandin receptors was increased in patients with severe COVID-19, as part of the cytokine storm developed by these patients [37]. Accordingly, our virtual screening identified prostaglandin E2 and other eicosanoids as N FL binders, these metabolites belonging to the COX-arachidonic -eicosanoids (including prostaglandin) pathways of inflammation and metabolism as an (early) host response triggered by the viral infection, in agreement with recent reports [16,38]. Treatment of patients hospitalized for mild and moderate COVID-19 with the COX inhibitor indomethacin helped the patients to recover [39]. Naproxen inhibited replication and reduced inflammation in a model of reconstituted lung epithelium [17]. In silico studies suggested that naproxen and acetaminophen may not only bind the N protein but also SARS-CoV-2 main protease, the RBD domain of spike and the polymerase [40].
The cytokine storm could also be considered as a possible driving factor for the expansion of neuropathies after severe COVID-19 infection, contributing to the chronic pain that appeared after acute infection recovery. The cytokine storm was at least partly mediated by (lung-resident and brain-penetrating) macrophages blocked in an M1 state that released pro-inflammatory cytokines as IL6 [41][42][43]. The neuropeptides identified in this work were known to have a role in inflammation: substance P is a mediator of inflammation and pain [41,44,45], while enkephalin had the opposite role to release pain [46,47]. In dopaminergic neurons, subpicomolar levels of substance P activated NADPH Oxydase-2 (NOX2) to increase ROS concentrations and subsequently synergistically activated the MAPK and NF-κB pathways, contributing to a potentiated pro-inflammatory cytokine production [47,48]. In addition, substance P and cholestocystokinin were vasoactive peptides [47]. We hypothesize that the sequestration by N of the identified ligands of N could have consequences for the host in terms of endothelial dysfunction, micro-thrombus in the brain and in other microvessels, with increased risk of stroke. Indeed, in the brain of SARS-CoV-2 infected non-human primates, the occurrence of a neuronal injury, brain micro-hemorrhages and hypoxia and rare viruses in brain-associated endothelium where N co-localized with the von Willebrand factor were observed [21]. It was also proposed that substance P and bradykinin, were likely to drive microvascular permeability, and be responsible for a phenomenon called «vasoactive peptide storm» as part of the development of COVID-19 pathology [48]. An agonist of cholecystokinin A (CCK-A) showed benefit in reduction of inflammation, among those hospitalized with moderate COVID-19 [49].
Therefore, we propose the hypothesis that the identified neuropeptides (and/or their fragments) could contribute to long-term inflammation in the CNS via N and possibly other mechanisms.

N Implication in Immunity, with Possible Long-Term Neurological Effects and Putative Viral Latency
N implication in immunity of the host against SARS-CoV-2 was expressed at multiple levels.
(i) T cells: Patients with persistent symptoms over 4 months following COVID-19 onset presented a lower frequency of CD8 + T cells expressing CD107a, a marker of degranulation, in response to Nucleocapsid (N) peptide pool stimulation, and a more rapid decline in the frequency of N-specific interferon-γ-producing CD8 + T cells [10]. (ii) Opioid peptides: Immune system and neuronal system cross-talk; this signaling is mediated by various molecules such as opioid peptides such as enkephalin. Enkephalins can impact lymphocytes proliferation, antibody synthesis. Enkephalin can enhance the release of pro-inflammatory cytokines like IL6 [50]. The cross-talk can take place between cytokines as CCL2 and opioid peptides and alter nociceptive synaptic transmission [51]. Thus, more work is required to test whether the enkephalin peptide identified here as a ligand of N may increase pain and immune disorders in the context of COVID-19, despite its physiological pain release function. (iii) Vasoactive peptides: Among the possible routes through which SARS-CoV-2 can invade the CNS, SARS-CoV-2 can directly invade the vagus nerve and retrograde into the CNS, or indirectly stimulate the enteric nervous system through immune pathways. Cholecystokinin could participate in this process. In addition, peripheral nerves may spread SARS-CoV-2 into the brain through the retro-neural route, including the olfactory nerve, trigeminal nerve, glossopharyngeal nerve and vagus nerve. Substance P is the main neuropeptide, neuromodulator and neuro-hormone of the trigeminal ganglion (TG), associated with nociception and inflammation through its receptor, the neurokinin-1. As observed with other viruses such as herpes or HIV, it was hypothesized that SARS-CoV-2 virus infection might become latent if it is acting through TG [41,45], involving substance P action.

N as a Mediator of Perturbed Metabolism via Its Ligands, with Possible Long-Term Effects
DHF and THF linked N to DHFR, an important enzyme that participates to DNA synthesis, being coupled to methionine metabolism. Additionally, folate exerted a protective role in the cardiovascular system [52,53]. Moreover, folate levels were usually low in patients with viral infection [54,55]. The reduction of viral replication by DHF/THF came at the expense of large perturbations of the host purine metabolism, nicely described in [35]; it is likely that the agonist of AMPK, AICAR would also have a dual effect on replication and lipid metabolism. The screening also identified Lauroyl-coA as a ligand of N, potentially linking N to beta oxidation of fatty acids and to triacylglycerol biosynthesis.

Zntpps as a Prototype of a Photoactive N Ligand for Photodynamic Therapy (PDT)
The binding site of hemin was deduced from docking [32], with a putative NLS signal (amino acids 258-268) at the edge of this site [33]. This site shared a number of predicted residues involved in the binding of the neuropeptides and DHF. Binding of ZnTPPS was demonstrated by spectroscopic changes of the porphyrin. ZnTPPS is one of the sensitizers that can release singlet oxygen and ROS upon light/ultrasound activation. Such compounds have been useful in the treatment of some cancers and proposed as antiviral treatment to inactivate SARS-CoV-2 [56].

Heme Sequestering by N and Its Potential Effect on NO and ROS Signaling
Dysfunctional, low hemoglobin levels, observed in patients with COVID-19 were possibly linked to well-characterized brain hypoxia [22,57]. Here, we further suggested that folate binding to N could also be responsible for an altered iron metabolism via perturbation in folate levels [35]. Moreover, we speculate that heme binding to FL could increase disruption of labile hemes from proteins such as guanylate cyclase, with potential effects in NO signaling, reduced oxygen/hypoxia sensing and/or impediment of their gaseous ligand binding such as O 2 , NO and CO, with all the known consequences of altered NO/O 2 signaling in the cardiovascular system (eNOS), in neurotransmission, control of serotonin levels (nNOS) and in immunity (iNOS). The following X-ray structures extracted from the Protein Data Bank (PDB) have been used: N-terminal domain of SARS-CoV-2: 6VYO [58], 7ACT [2].

Materials and Methods
A computational model of the full-length nucleoprotein of SARS-CoV-2 was generated with the ColabFold [59] version of AlphaFold 2 [27] using an API hosted at the Södinglab based on the MMseqs2 server [60] for the multiple sequence alignment creation. 300 sequences were aligned on uniprot entry A0A6C0T6Z7 as query sequence and no template (419 residues). Predicted lDDT is highest in the NTD and CTD regions that were previously determined by experimental methods (PDB entries 7act, 6wzo).
MD simulations were performed based on Alphafold and initial minimization and equilibration steps for 21 ns in implicit solvent. The wt protein was prepared with Discovery Studio 2022 (Dassault Systèmes BIOVIA, Discovery Studio Modeling Environment, Release 2022, (San Diego, CA, USA): Dassault Systèmes 2022) the models were first protonated at pH 7.4, typed with CHARMM36 force field and the protein was placed in an orthorhombic box solvated with a pre-equilibrated solvent containing TIP3 waters, with a 15Å minimum distance from boundary. Two trajectories were generated with the NAMD 2.13 engine and CUDA acceleration as implemented in Discovery Studio 2022, by varying the Random Seed Number. The MD simulations were run for 100 ns each under constant pressure (NPT ensemble) at 300 K using a Langevin thermostat and piston, for temperature and pressure control, respectively, and long-range electrostatic interactions were computed using Particle Mesh Ewald (PME) algorithm [61]. The trajectory was analyzed with the Analyze Trajectory Protocol to compute the RMSD (Root Mean Square Deviation) and RMSF (Root Mean Square Fluctuation).

Model of N Variants
This wt model was then used to generate the BA4 and BA5 variants of N with the MODELER algorithm [62] as implemented in Discovery Studio version 2022 (Dassault Systèmes BIOVIA, Discovery Studio Modeling Environment, Release 2022, San Diego: Dassault Systèmes 2022). For each of the two mutants BA4 and BA5, a 100 ns trajectory was also generated with the same protocol as described for the WT protein.
All the cavities were obtained with the Define Site from Receptor Cavities of Discovery Studio 2022 as shown in Figure 2.

Ligand Screening and Docking to wt N
In silico screening of ligands that bind the protein structures was performed using the Sigma catalog data base and assessed in a two-step protocol. For this screening, we first used the docking program Libdock of Discovery Studio version 2021, with the exclusion of ligands with molecular weights above 2000 and according to the binding site(s) defined by cavity detection. A site-docking was then performed using Cdocker with energy minimization (Discovery Studio version 2021) to more precisely identify the ligand binding mode. For each protein, the most representative pose of the ligand was selected. The resulting protein-ligand complexes were finally refined by a molecular dynamics simulation using the CHARMm force field [63] and the standard dynamic cascade protocol of Discovery studio version 2021. This protocol started with a first minimization of 1000 steps using the Steepest Descent algorithm and a RMSD gradient of 1 Å, followed by a second minimization of 2000 steps using the Adopted Basis Newton-Raphson algorithm and an rmsd gradient of 0.1 Å. The third step involved heating from 50 K to 300 K, with a fourth step of equilibration during 1 ns and a fifth step, production. The time of the production step was initially set at 10 ns, but extension to 20 ns was applied if the ligand was not stable. Three replicas were carried out for each complex. For each trajectory, the displacement of the ligands was studied by RMSD calculation. The representative structure (i.e., with the smallest average rmsd from all other structures of the cluster) of the largest cluster of each complex was selected.

Substance P (1-7) Bound to BA5 Variant
An homology model of BA5 in complex with substance P (1-7) was built with MOD-ELER, using the wt-substance P (1-7) as a template. For each of the two mutants BA4 and BA5 in complex with the substance P (1-7), a 100 ns trajectory was also generated with the same protocol as described for the wt protein.

Proteins Expression and Purification
The N-terminal domain of the nucleoprotein of SARS-CoV-2 (residues 50-173) with an N-terminal His6 tag was cloned in pET-28a vector (pET28a-His6-NP-NTD). Heterologous expression at 15 • C for 16 h in E. coli BL21 bacterial strain (DE3 (NEB, Evry, France). The recombinant protein (15 kDa) found in the soluble fraction was purified on a Ni 2+ -NTA affinity column and SP sepharose ion exchange chromatography, and presented a single band revealed by SDS PAGE. The yield was 1.5 mg of pure protein per liter of culture.
The expression of the FL protein was performed as described above. The insoluble fraction of the pellet was resuspended in lysis buffer containing 6 M urea for an hour at room temperature, the insoluble part was removed by centrifugation over 10 min. The purification involved an affinity Ni 2+ column, with a first washing at 3 M urea, followed by washings in a 20 mM Tris-HCl pH = 7.9 and 1 M NaCl, and then increasing concentrations of imidazole (10, 20 and 50 mM). Elution was performed with 500 mM imidazole. Imidazole was then removed from the solution by dilution-ultrafiltration (centricon molecular weight cutoff 10,000, Millipore). The protein was then left at 4 • C for three days for renaturation in a 20 mM Tris-HCl buffer at pH = 7.9, 100 mM NaCl.

SPR Experiments
The recombinant FL protein was tagged with six His for purposes of identification and affinity purification, and the same principle was exploited for capture of the recombinant proteins on the sensor surface. Poly-histidine is a commonly used tag that can chelate with Ni 2+ ions in complex with nitrilotriacetic acid (NTA, Biacore, Uppsala, Sweden), providing a convenient approach for capturing tagged constructs on Sensor Chip NTA. The FL protein was injected at a flow rate of 30 µL/minute; the ligand was injected after the protein capture; then, the protein or protein-ligand complex was removed from the surface by applying 350 mM EDTA, with a very good reproducibility. The experiments were usually performed in 20 mM HEPES buffer (Biacore), 100 mM NaCl and 0.05% Tween. The FL protein was captured at a concentration of 0.1 µM as the signal decreased at higher FL concentrations, presumably because the protein adopts higher oligomeric state(s). The NTD wt and mutant proteins were captured at a concentration of 1 µM for shorter times. In both FL and NTD proteins, signals of 5000 to 8000 RU were captured on the surface. The NTD proteins were always used as freshly prepared solutions. The ligands concentration range was usually in the range 0.3-10 µM.

Dynamic Light Scattering
The experiments were performed on a Malvern nanosizer apparatus. The temperature was set at 20 • C, and 10 scans with a duration of 10 s each were acquired in duplicate for each time and sample. The size distribution in the intensity of the scattered light was obtained using the Cumulants method from the instrumental software, yielding the hydrodynamic diameter. The N NTD concentration was in the range of 40-60 µM, the FL concentration was in the range of 1 to 10 µM in 20 mM Tris buffer at pH = 7.9 containing 100 mM NaCl. The melting experiments were performed at a heating rate of 1 • C/min over the range 30 • C to 65 • C in sealed disposable cuvettes. The apparent melting temperature, Tm, was determined by the first derivative of the melting curve.

Conclusions
The abundant N protein potentially could disrupt many signaling pathways, through sequestration of signaling molecules as the neuropeptides identified here or other useful metabolites at its large C-terminal cavity. An increased flexibility of N seemed to emerge in the BA5 N variant as compared to the BA4 one, which may both speed up replication and enhance ligand binding, with possible cooperative or long-range effects of the mutations/deletion. Although these neuropeptides usually signal in the pico-to nanomolar concentration range in physiological conditions, their concentrations in the blood, in the CNS and in TG can increase locally, in particular in COVID-19 patients [64][65][66][67][68][69][70], in line with the hypothesis suggested by this study. N may affect the cell metabolism via folate and AMPK and neurotransmission via neuro/vasoactive peptides. The physiological function of these ligands or their fragments would probably be altered consecutive to their association with N.
Altogether, the host likely senses the N protein as a "danger hub"; any possible mean seems to be undertaken to neutralize N. Therefore, instead of letting valuable neuropeptides be sequestered by N, blocking N with antivirals may decrease symptoms associated with long COVID-19, which accounts for 5-20% of the patients [22,71]. Based on the biological functions of the neuropeptides, we speculate that N could be involved in pain, neurological brain fog and inflammation and possibly immune imbalance and latency. Current NSAIDs as naproxen and indomethacin are readily available and bind to N; extended ligands could be also adequate to fill the C-terminal cavity of N. The hypothesis raised in this discussion offers a comprehensive picture that may link multiple long-haul COVID-19 symptoms through N, which is yet to be demonstrated by clinical data.


GmARP is Related to the Type III Effector NopAA to Promote Nodulation in Soybean (Glycine max)

Type III effectors secreted by rhizobia regulate nodulation in the host plant and are important modulators of symbiosis between rhizobia and soybean (Glycine max), although the underlying mechanisms are poorly understood. Here, we studied the type III effector NopAA in Sinorhizobium fredii HH103, confirming its secretion into the extracellular environment under the action of genistein. The enzyme activity of NopAA was investigated in vitro, using xyloglucan and β-glucan as substrates. NopAA functions were investigated by the generation of a NopAA mutant and the effects of NopAA deficiency on symbiosis were analyzed. Soybean genes associated with NopAA were identified in a recombinant inbred line (RIL) population and their functions were verified. NopAA was confirmed to be a type III effector with glycosyl hydrolase activity, and its mutant did not promote nodulation. Quantitative trait locus (QTL) analysis identified 10 QTLs with one, Glyma.19g074200 (GmARP), found to be associated with NopAA and to positively regulate the establishment of symbiosis. All these results support the hypothesis that type III effectors interact with host proteins to regulate the establishment of symbiosis and suggest the possibility of manipulating the symbiotic soybean–rhizobia interaction to promote efficient nitrogen fixation.

INTRODUCTION
Soybean is an important food crop, could provide rich plant protein and oil for human. With the increasing demand for protein and oil, more and more attention has been paid to soybean agricultural production . In agricultural production, the yield of soybean is increased by applying a large amount of industrial nitrogen fertilizer (Wang et al., 2020a). However, excess application of nitrogen fertilizers has had adverse effects on the environment, leading also to economic and health concerns (Jordan et al., 2011;Chen et al., 2017;Guo and Wang, 2021). The use of symbiotic nodule-forming soybean maintains soil fertility and sustainable crop production, leading to reduced requirements for chemical N fertilizers in agricultural systems (Goyal et al., 2021). Many native Rhizobium species have been isolated and inoculated into soybean agricultural production to reduce the amount of applied N fertilizer . Rhizobia inoculations benefit the production of soybean as they promote plant growth by producing hormones, antibiotics, vitamins, and siderophores (Roy et al., 2020); in addition, rhizobia can enhance disease resistance (Plett et al., 2021;Zboralski et al., 2022).
The establishment of symbiosis system is a complex process, which requires signal exchange between soybean and rhizobia (Roy et al., 2020). In rhizobia, six secretion systems, namely, type I to VI systems, have been identified and classified according to their function (Freiberg et al., 1997;Krehenbrink and Downie, 2008;Pukatzki et al., 2009;Tampakaki, 2014;van Ulsen et al., 2014). Each secretion system plays a different role in symbiotic nitrogen fixation, with the type III secretion system (T3SS) being the most widely studied (Kambara et al., 2009;Miwa and Okazaki, 2017). The T3SS secretes nodulation outer proteins (Nops) into host cells which play key roles in the establishment of symbiosis. Several Nops have been identified in Sinorhizobium fredii HH103, including NopD, NopL, NopM, NopP, and NopT (Ausmees et al., 2004;Bartsev et al., 2004;Dai et al., 2008). NopD and NopM have been identified by LC-MALDI and LC-ESI (Rodrigues et al., 2007), and protein domain prediction has shown that the C-terminal region of NopD contains a domain homologous to the ubiquitin-like protease Ulp1, thus suggesting that NopD may combine with SUMO (small ubiquitin-like modifier)-conjugated proteins, leading to the removal of the SUMO conjugate in HH103-infected soybean (Xiang et al., 2020;Ratu et al., 2021). NopM is an E3 ubiquitin ligase and may target the MAPK pathway in host plants during nodule formation (Xin et al., 2012). NopL and NopP are substrates for soybean kinases (Skorpil et al., 2005;Zhang et al., 2011a;Ge et al., 2016). The NopL mutant was found to inhibit nodule formation, while the NopP mutant increased nodule numbers (Wang et al., 2018;Zhang et al., 2018;Zhang et al., 2021a). These reports indicate that the various type III effectors have different roles in symbiosis. The analysis and identification of type III effector proteins in soybean could provide a good basis for understanding the functions of these proteins in symbiosis. The release of the soybean (Glycine max) reference genome (Lam et al., 2010;Schmutz et al., 2010) has allowed the identification of many novel genes based on quantitative trait loci (QTLs) (Qi et al., 2014a;Lu et al., 2017). The Rj2 protein was cloned and found to recognize T3SS proteins . GmNNL1, an R gene that interacts with the rhizobial effector NopP, was identified by GWAS (Zhang et al., 2021a). NopD promotes HH103 infection of soybean by direct or indirect regulation of GmPP2C, and different haplotypes of GmPP2C were found to vary their expression patterns in response to NopD (Wang et al., 2020b). GmRPK was found to be related to NopL and to promote HH103 infection . Although many proteins related to type III effectors have been identified in soybean, much work is still needed to uncover the mechanisms by which these proteins function.
In this study, we aimed to identify the genes associated with NopAA in regulating nodulation. We identified NopAA as a glycosyl hydrolase that hydrolyzed xyloglucan and β-glucan into glucose. Nodule tests in 13 soybean natural varieties were used to study the effects of NopAA mutation, indicating that NopAA positively influenced symbiosis. We identified 10 QTLs related to NopAA using the Dongnong594 and Charleston germplasms. Further haplotypic, genetic, and transgenic analyses showed that Glyma.19g074200 (GmARP) responded to NopAA to regulate nodulation. These findings thus provide useful information on the role of NopAA in nodulation, and the identification and study of soybean host genes can provide support for the effective utilization of symbiotic nitrogen fixation in soybean.

Strains, Vectors, and Primers
S. fredii HH103 and Escherichia coli DH5α and BL21 (DE3.1) were used in this study (Supplementary Table S2). HH103 was cultured in or on TY medium supplemented with rifampicin at 28°C, and all E. coli strains were cultured in or on LB media supplemented with appropriate antibiotics at 37°C. The plasmids used for the construction of mutants and the primers used for PCT and qRT-PCR are listed in Supplementary Table S1. Antibiotics were used at the following concentrations: 50 μg/ ml rifampicin; 50 μg/ml kanamycin; 34 μg/ml chloramphenicol; 50 μg/ml spectinomycin; 50 μg/ml streptomycin.

Construction of the NopAA Mutant
The NopAA mutant was constructed by triparental mating as previously described  with the kanamycinresistant DNA fragment inserted into the coding sequence of NopAA by homologous recombination. The mutant was confirmed by PCR and Southern blotting (Chen and Kuo, 1993).

qRT-PCR and Analysis of Nodulation Outer Proteins
HH103, the NopAA mutant, and the TtsI mutant were grown in YM medium (Yeast Mannitol Broth) in the presence or absence of genistein, and extracellular proteins were extracted as previously described (Jim é nez-Guerrero et al., 2015). Extracted RNA was subjected to qRT-PCR analysis by a Roche LightCycler LC480 (Roche, Switzerland), and the relative expression levels of NopAA were normalized using 16S rRNA. All sample collections were analyzed using three independent biological replicates as well as three technical replicates.

Enzyme Activity Assays
The NopAA genes were then independently subcloned into pET28b as BamH I-Sal I fragments to create pET28b-NopAA. This vector was used for active protein purifications as described (Xin et al., 2012). The His-NopAA protein was purified from BL21 (DE3) E. coli after induction with 0.1 mM isopropylthio-β-galactoside using Ni-NTA resin beads (Qiagen, Hilden, Germany). The purified protein was identified by western blotting with an anti-His antibody and its enzyme activity was measured using pure xyloglucan and beta-glucan (Megazyme, Bray, Ireland) as substrates. Activity assays were performed by glucose as the standard and a iMark Microplate Absorbance Reader (Bio-Rad, United States) as previously described (Grishutin et al., 2004).

Agrobacterium Infiltration Assays
Agrobacterium-mediated transformation of tobacco leaves was performed as described (Wang et al., 2011). A. tumefaciens EHA105 containing pSoy1-NopAA, the empty vector pSoy1, Frontiers in Genetics | www.frontiersin.org May 2022 | Volume 13 | Article 889795 and pSoy1 35S-NopM were transformed into leaves of 30-day-old Nicotiana benthamiana. The OD 600 of the A. tumefaciens culture was adjusted to 0.2 with buffer (Ma et al., 2015). Infiltrated leaves were harvested each day for analysis of electrical conductivity by a conductivity meter (Lei-ci, China) (Wang et al., 2020a).

Nodulation Tests
The RILs population were planted in Xiangyang Farm, Harbin, Heilongjiang Province (Harbin, latitude 45°450″N, longitude 126°380″E) in 2016, and nodulation tests were performed in 2017. 100 soybean germplasms used in the study came from different ecoregions were planted in Xiangyang Farm in 2019, and nodulation tests were performed in 2020. Sterilized soybean seeds (15 seeds from each variety) were grown under greenhouse conditions in sterilized vermiculite containing a low-nitrogen nutrient solution . All plants were grown in a greenhouse at 25°C with a photoperiod of 16 h. At the Vc stage (unifoliolate leaves expand), all the plants were inoculated with HH103 and the NopAA mutant. Nodulation in the germplasms was evaluated by nodule number (NN) and nodule dry weight (NDW) (28 days post-inoculation, dpi). Three different biological replicates and 15 plants for each replicate were conducted and the t-test was used to detect the statistical significance of differences in NN and NDW.

Observation of Rhizobial Infection Events
DN594 and Charleston were inoculated with GFP-conjugated HH103 and NopAA mutant strains. One-centimeter samples of infected lateral roots were collected at 1 dpi, comprising ten lateral roots per plant, and infection in the entire plant was evaluated using confocal fluorescence microscopy (Zeiss LSM700, Germany). The number of infection events was then multiplied by the length of the roots. Infection events included foci, infection threads in epidermal cells (ITs), and infection threads extending into the cortical cells (rITs) . Ten independent lateral root segments from each plant were used for each biological replicate and three biological replicates were used for each condition.

Quantitative Trait Locus Mapping and Identification of Genes in Quantitative Trait Locus Regions
Charleston and DN594 were used to construct the RIL population, using a genetic map that had been constructed in a previous study (Qi et al., 2014b). We used WinQTL Cartographer and composite interval mapping methods to map QTLs related to nodulation in the RIL population. QTLs for NN and NDW were defined by LOD scores greater than 3.0 (Brensha et al., 2012;Wang et al., 2020b). The reference genome of Williams 82 was used to identify genes within the QTL regions, and candidate genes were annotated according to the gene annotation information of Williams 82 (Xin et al., 2016).

Analysis of the Expression Patterns of Candidate Genes Related to NopAA
qRT-PCR was performed to analyze the expression patterns of candidate genes in Charleston and DN594 inoculated with HH103, the NopAA mutant, and the control. Roots were harvested at several time points and rapidly submerged in liquid nitrogen, after which the total RNA was extracted by the TRIzol method. cDNA was synthesized using the PrimeScript ™ RT Reagent Kit (Takara Biotech Co., Beijing,

Hairy-Root Transformation of Soybean
The A. rhizogenes strain K599 containing pSoy10-GmARP-GFP, pSoy10-GFP, pB7GWIWG2-GmARP-DsRed, and pB7GWIWG2-DsRed were used for hairy-root transformation as described previously (Kereszt et al., 2007). GFP and DsRED were used as markers to detect positive hairy roots using LUYOR-3415RG. Hairy roots were inoculated with HH103, the NopAA mutant, and the control to identify the influence of gene overexpression and interference on nodulation. All nodules were harvested at 28 dpi and transcript levels of GmARP were measured using qRT-PCR. Three different biological replicates and 15 plants for each replicate were conducted and the t-test was used to detect the statistical significance of differences in NN and NDW.

Haplotype Analysis
Haplotype analysis of GmARP was performed in 100 soybean natural varieties. The genomic sequence of GmARP, including the coding sequence and the 3000-bp promoter sequence, was obtained from the resequencing of the genomes of the 100 varieties, and the genomic sequences were subjected to local BLAST analysis to obtain SNP information (Wang et al., 2020a). Haploview 4.2 software (Cambridge, MA, United States) was used for analysis using the Haps Format module, and GraphPad Prism 8 software was used to analyze the correlations between the GmARP haplotypes and nodule traits in the 100 natural varieties .

NopAA is Induced by Genistein and Secreted by Rhizobia
A phylogenetic analysis was performed to study the evolutionary relationships of NopAA in different rhizobial species. This showed that the NopAA tree was separated into three branches (blue, green, and pink regions in Figure 1A). The blue region contained only four Sinorhizobium species, the green region contained only three Bradyrhizobium species, and the pink region contained only two Mesorhizobium species ( Figure 1A). The phylogenetic analysis showed that there is obvious diversity in NopAA among different rhizobial genera. Previous RNA-seq results showed that NopAA might be induced by genistein; thus, to further study whether NopAA was a type III effector, qRT-PCR and nodulation analyses were performed. The qRT-PCR results showed that genistein increased NopAA expression significantly in HH103 compared with rhizobia not Frontiers in Genetics | www.frontiersin.org May 2022 | Volume 13 | Article 889795 treated with genistein but did not induce NopAA expression in the NopAA mutant or the TtsI mutant regardless of the presence or absence of genistein ( Figure 1B). After extraction of the nodulation outer proteins, NopAA was detected in the supernatants after genistein treatment but not in the supernatants of the NopAA mutant or TtsI mutant outer proteins induced by genistein ( Figure 1C). These results indicated that NopAA could be induced by genistein and TtsI and could be secreted into the extracellular milieu of the rhizobia.

NopAA has Glycosyl Hydrolase Activity and Does not Induce Cell Death in Tobacco
To determine whether NopAA had glycosyl hydrolase activity, 6 × His-tagged NopAA was expressed and purified by Ni-affinity chromatography. The molecular weight of the His-NopAA fusion protein was~28 kDa ( Figure 1D). Xyloglucan and β-glucan were used as substrates to analyze hydrolase activity, showing that NopAA hydrolyzed both xyloglucan and β-glucan directly to sugars ( Figure 1D). The Phytophthora sojae apoplastic effector PsXEG1 and NopAA belong to the glycosyl hydrolase family, which is known to induce cell death in tobacco leaves. To study the effect of NopAA on tobacco leaves, we infiltrated NopAA into tobacco leaves using Agrobacterium tumefaciens EHA105. In contrast to PsXEG1, NopAA did not trigger cell death ( Figure 1E). These results indicated that NopAA had glycosyl hydrolase activity and did not induce cell death in tobacco.

NopAA Mutant Does not Promote Nodulation in Several Soybean Germplasms
In this study, we selected 13 soybean germplasms to study the effect of NopAA on nodulation, including five cultivated, six landrace soybean, and two wild soybean germplasms (Supplementary Figure S1). Nodulation tests showed that while the NopAA mutant significantly reduced the nodule numbers in 11 of the germplasms, the nodule numbers did not differ significantly in Dongnong594 (DN594) and Qingdou compared with the wild strain HH103. Inoculation with the NopAA mutant led to lower nodule dry weights in 12 of the 13 germplasms except for DN594, compared with the wild strain HH103. Nodulation tests of DN594 and Qingdou showed that the NopAA mutant did not influence nodule numbers, nor did it alter the change in the dry weight of DN594. These results indicated that NopAA did not negatively influence symbiosis in different types of soybean germplasms, and that differences in the genetic backgrounds of the germplasms might be the explanation for the variations in nodule traits associated with NopAA. This suggested that these genetic differences could be used to identify genes or loci associated with symbiosis.

The Influence of NopAA on Dongnong594 and Charleston
Based on the results of the nodule tests on various soybean germplasms, DN594 and Charleston were selected for further investigation of the role of NopAA in symbiosis. Compared with the wild strain HH103, the NopAA mutants caused a significant reduction in both the nodule number (NN) and nodule dry weight (NDW) in Charleston but did not affect DN594 (Figures 2A,B). In the mature nodule, the infection zone and infection cells did not differ significantly between the two germplasms after inoculation with HH103, and inoculation with the NopAA mutant did not change this (Figure 2A). Because of the differences in the nodule traits of DN594 and Charleston, we further investigated the infection events after inoculation with GFP-labeled HH103 and the NopAA mutant.
The results showed there were more infection events in curled root hairs in Charleston compared with DN594 3 days after inoculation with GFP-labeled HH103 compared with the wild strain, while the GFP -labeled NopAA mutant induced fewer infection events in Charleston but had no significant effects on DN594 ( Figure 2C). This implied that the type III effector NopAA might influence nodulation through involvement in infection events and that the genetic differences between DN594 and Charleston might result in different symbiotic phenotypes, with some of the genetic differences being associated with NopAA.

Quantitative Trait Locus for Nodule Number and Nodule Dry Weight Related to NopAA in Soybean Recombinant Inbred Line Populations
The soybean RILs were derived from the cross and continuous self-cross of DN594 and Charleston and were used to identify genes related to NopAA. Nodulation tests were performed on 150 RILs after inoculation with the wild strain and the NopAA mutant, with the results showing that the NopAA mutant could significantly reduce both the NN and NDW in whole RILs compared with HH103, further suggesting that NopAA plays a positive role in symbiosis. WinQTL Cartographer was used to identify the QTLs underlying nodule-related traits using a composite interval mapping method. Two main-effect QTLs underlying NDW and five QTLs underlying NN were identified with the background inoculated with S. fredii HH103 (Table. 1). Two QTLs, QNDW3-1 and QNDW7-1 were associated with NDW and explained 4.40% and 2.42% of the phenotypic variation, respectively. Five QTLs, QNN9-1, QNN13-1, QNN13-2, QNN16-1, and QNN17-1. were associated with NN, explaining 1.84%, 0.99%, 0.39%, 5.89%, and 6.13% of the phenotypic variation, respectively. A total of 33 candidate genes were found in these seven loci (Supplementary Table S3).
In the background inoculated with the S. fredii NopAA mutant, two QTLs underlying NDW were located on chromosome 19, while only one QTL underlying NN was identified on chromosome 19 ( Table 2). The QTLs QNDW19-1 and QNDW19-2 explained 0.15% and 2.59% of the phenotypic variation, respectively QTLs QNDW19-2 and QNN19-1 were in the same region, with very small R 2 values. Twenty genes were identified in these loci (Supplementary Table S3). As no comparable region was identified in the backgrounds of S. fredii HH103 and the S. fredii NopAA mutant and as QNDW19-2 and QNN19-1 were in the same region, we considered that this QTL region might contain the key candidate genes for interaction with NopAA. Gene functional annotation showed that gene Glyma.19G074000 encoded a nodulin protein, while Glyma.19G073900 encoded a DNA lyase, and Glyma.19G074200 was co-expressed with genes in the rootspecific co-expression subnetwork. Glyma.19G073800 and Glyma.19G074100 had no specific functional annotations.

qRT-PCR Verification of Candidate Genes
The candidate genes were verified by qRT-PCR to confirm their association with NopAA, using the roots of both DN594 and Charleston, as well as non-inoculated roots (Figure 3). Among these genes, Glyma.19G074200 was confirmed to be upregulated with an almost two-fold change after inoculation with S. fredii HH103 but not with the S. fredii NopAA mutant 6 h after inoculation. This supported the association between Glyma.19G074200 (named GmARP) and NopAA. The remaining genes were not significantly affected by the wildtype and mutant rhizobial inoculations.

The Effects of GmARP Silencing and Overexpression on Nodulation
To further investigate the role of GmARP in soybean nodulation, the Agrobacterium rhizogenes strain K599 carrying pB7GWIWG2(II)-Dsred-GmARP (for RNA interference), pSoy10-GmARP-GFP (for overexpression), and the corresponding empty vectors were used for the transformation of soybean transgenic hairy roots. Silencing and overexpression were confirmed by qRT-PCR after transformation (Supplementary Figure S2). Reduced NN and NDW values were seen after GmARP silencing compared with the controls inoculated with HH103 and the NopAA mutant, respectively ( Figure 4A). In the silenced hairy roots, the NopAA mutant did not induce significantly different NN or NDW compared with HH103 ( Figure 4B). In the overexpressed hairy roots, GmARP overexpression resulted in increased NN and NDW values compared with the controls inoculated with HH103 and the NopAA mutant, respectively ( Figure 4B). We found that overexpression of GmARP led to elevated NN and NDW in comparison with the control hairy roots inoculated with the NopAA mutant. However, both the NN and NDW after inoculation with the NopAA mutant were still lower than HH103 in the overexpressed hairy roots, suggesting that overexpression could partially compensate for the loss of Note: Chr., chromosome; LOD, log-of-odds; QTL, quantitative trait loci; LOD score cutoff of major QTLs was determined by permutation tests (1,000 times; p < 0.05).
Frontiers in Genetics | www.frontiersin.org May 2022 | Volume 13 | Article 889795 FIGURE 3 | Analysis of the relative expression levels of candidate genes in the QTL interval by qRT-PCR. Relative expression of genes was calculated by the 2 −ΔΔCt method, and the relative expression levels of each gene were normalized using GmUNK1 (Glyma.12g020500). Data are means ± SE of three replicates. NopAA in HH103 during symbiotic formation. These results indicated that GmARP may be involved in signal transduction after recognition of NopAA by the soybean, leading to the regulation of nodulation.

Haplotype Analysis Suggests That GmARP Regulates Nodule Number and Nodule Dry Weight in Soybean
Based on the observed differences in the recognition of NopAA, haplotype (Hap) analysis of GmARP was performed in 100 soybean natural varieties using Dnasp5.0 software. Using the results of the resequencing of soybean natural varieties, eight Haps were identified. Two Haps containing more than 10 accessions were considered the dominant Haps, and only one SNP and one indel were observed in the promoter regions of Hap1 and Hap2 (Figures 5A,B). What is even more remarkable is that there was no difference in the NN and NDW values in the Hap1 accessions inoculated with HH103 compared with the NopAA mutant; however, in the Hap2 accessions, the NopAA mutant induced both lower NN and NDW than HH103 ( Figure 5C). To identify whether the difference in the promoter region caused the expression difference between Hap1 and Hap2, qRT-PCR was used to analyze the expression of GmARP at 12 h post-inoculation with HH103 and the NopAA mutant in some varieties of Hap1 and Hap2. The expression patterns suggested that the relative expression level of GmARP did not differ between several Hap1 accessions (Heihe13 and Suinong15) regardless of the presence of NopAA in HH103 ( Figure 5D). However, in the Hap2 accessions, HH103 induced significant expression of GmARP compared with the NopAA mutant. The haplotype analysis further confirmed that GmARP was associated with the type III effector NopAA to regulate nodulation in soybean.

DISCUSSION
In this study, we confirmed that NopAA was a type III effector and could hydrolyze xyloglucan and β-glucan directly into glucose. Nodulation tests indicated that NopAA did not have a negative effect on nodulation, and the soybean gene GmARP was found to be associated with NopAA in the soybean RIL population. There are few studies on NopAA function and the mechanisms involved, and few address the question of its role in symbiosis. Our study of NopAA Frontiers in Genetics | www.frontiersin.org provides an essential foundation for understanding the role of NopAA and its associated signaling pathways in symbiotic nitrogen fixation.
The type III effectors of rhizobia have obvious characteristics, that is, they are induced and regulated by flavonoids and TtsI, and the most important is that they can be secreted into the extracellular of rhizobia (Teulet et al., 2022). Similar to the reported type III effectors, such as NopD, NopL and NopM (Zhang et al., 2011b;Xin et al., 2012;Wang et al., 2020a), NopAA can also be induced by flavonoids and its expression can be regulated by TtsI, these confirmed NopAA is a type III effector. In the previous study, several type III effectors can act as enzymes and have biochemical activities (Teulet et al., 2022), and elucidating the biochemical functions of type III effectors, particularly those that are conserved among different species, has greatly enhanced our understanding of the mechanisms underlying bacterial pathogenesis Zboralski et al., 2022). By analysis of NopAA amino acid sequence, we found that NopAA belonged to the Glycoside hydrolase 12 (GH12) family and had the potential to hydrolyze polysaccharides (Zhang et al., 2021b). Nodulation requires the assistance of cell-wall degrading enzymes, such as hemicellulases, pectinases, polygalacturonases, glucanases, cellulases, xyloglucanases, and pectinases, which allow root penetration by beneficial microorganisms, such as Rhizobium, Frankia, and the arbuscular mycorrhizal (Mateos et al., 1992;Zhang et al., 2021b;Tsyganova et al., 2021). GH12 is a cell wall-degrading enzyme that, along with the bacterial effector Cif and the Phytophthora sojae apoplastic effector PsXEG1, have been identified as playing pivotal roles in pathogen survival and systemic infection processes (Ma et al., 2015;Xia et al., 2020). Identifying the host protein interactors could elucidate the mechanisms underlying host-microbe interactions. NopAA was identified as a glycoside hydrolase belonging to the GH12 family. In Phytophthora sojae, PsXEG1 is also a glycoside hydrolase 12 family member with glycoside hydrolase activity and hydrolases the same substrates as NopAA. Because PsXEG1 is an effector of P. sojae, the replication of a P. sojae virus could be promoted using a PsXLP1 decoy that was similar to PsXEG1 but without its enzymatic activity. Additionally, PsXLP1 can protect PsXEG1 from GmGIP1 binding in vitro and in planta (Ma et al., 2015;Ma et al., 2017). We did not find any proteins similar to NopAA in the S. fredii HH103 genome, suggesting that HH103 does not infect soybeans using the same decoy pattern as P. sojae by NopAA. The investigation of the effect of NopAA on tobacco leaves did not show any evidence of cell death. Rhizobia may be involved in soybean immune signaling pathways during infection and the establishment of the symbiotic relationship. However, the establishment of an effective symbiosis does not require as strong or prolonged an immune response as occurs against pathogens. These results differ significantly from the actions of PsXEG1 in tobacco leaves. Thus, we propose that the role of NopAA may be different from that of PsXEG1. In this study, xyloglucan and β-glucan were found to be hydrolyzed to glucose; both xyloglucan and β-glucan are important components of cellulose and hemicellulose in plant cell walls, and their hydrolysis could promote the entry of rhizobia into host cells. These results are consistent with the observations of infection sites as fewer infection sites were observed with the NopAA mutant, indicating that NopAA could degrade the cell walls of soybean root hair cells, thus promoting infection and the formation of the infection thread.
A total of 10 QTLs were identified in this study, two of which were related to NDW, and five were related to NN with HH103, with two QTLs related to NDW and only one QTL related to NN seen with the NopAA mutant. No overlapping QTL loci were found when comparing inoculated HH103 and NopAA mutant, because NopAA mutation obviously changed the NN and NDW of DN594, Charleston and RIL populations, so no overlapping QTL could be identified, this result was similar to previous studies of other effectors, such NopL and NopT Ni et al., 2022). Of the 10 identified QTLs, two were found to overlap with the reported symbiotic loci Ns1-3 (Hwang et al., 2013) and Nw1-9 (Nicol á s et al., 2006), respectively. These results indicated the accuracy of the QTL mapping and genetic analyses in the DN594 and Charleston RIL populations. The NopAA mutant was derived from the wild strain HH103 and was found to induce fewer NNs and lower NDW, as well as to generate fewer infection sites. Analysis of the expression patterns of several symbiotic marker genes showed that the NopAA mutant could also change the expression pattern of genes, such as PR1. It is possible that the interaction network might change without the involvement of NopAA. These results suggested that the soybean response to the NopAA mutant differed significantly from its response to the wild strain, leading to the observed QTL differences between HH103 and the NopAA mutant.
In the overlapped QTLs, Glyma.19G074200 was found to be associated with NopAA using gene annotation and qRT-PCR. Glyma.19G074200 is a root-specific co-expression gene belonging to the uncharacterized protein At3g27210 protein family. In Arabidopsis, the At3g27210 protein interacts with AtRKL1 (Tarutani et al., 2004) and participates in the response to pathogen infection. Soybean RNA-seq results confirmed that Glyma.19G074200 was involved in symbiosis during Bradyrhizobium japonicum infection, although the function of the gene was not determined (Libault et al., 2010). In soybean, the function of this gene has not been studied, in this study, overexpression of GmARP promoted nodule formation after inoculation with either the wild strain or the NopAA mutant but, after GmARP silencing, no significant differences in either NN or NDW between plants inoculated with HH103 or the NopAA mutant were observed. These results suggest that GmARP is a positive regulator of nodule formation and that it mediates NopAA signaling in plants, the mechanism of the synergistic regulation of GmARP and NopAA in symbiosis remains to be further studied. We demonstrated that NopAA had glycosyl hydrolase activity and may promote rhizobia infection by hydrolyzing the host cell wall, and that host recognition of NopAA activated the expression of symbiotic genes. When HH103 was used for infection, NopAA induced GmARP expression to complete the establishment of symbiosis. We have two hypotheses concerning the mechanism of this activated expression: the first is that NopAA induces the expression of GmARP by activating its downstream signals while the second is that the products resulting from NopAA hydrolysis of cell wall polysaccharides induce signaling changes in the host cells, which further activate the expression of GmARP. However, these hypotheses require further experimental verification.

CONCLUSION
We demonstrated that NopAA is a type III effector with glycosyl hydrolase activity. Investigation of the effects of the NopAA mutant on nodulation showed that the mutant did not promote nodulation. Using QTL mapping and analysis of soybean genes associated with NopAA, we found that the soybean protein GmARP positively regulated nodulation through its association with NopAA. This study provides an essential reference for analyzing the function of NopAA and provides support for efficient utilization of symbiotic nitrogen fixation in soybean agricultural production.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
JHW, QC, MY, and DX conceived the study and designed and managed the experiments. CM and SM performed trials and collected data. JXW, HZ, HF, and YW completed statistical analyses of phenotypic data. CL provided plant lines. JXW, QC, MY, DX, and CM all contributed to writing the paper.