1 Introduction

Mechanosensitive (MS) ion channels represent pivotal constituents in mechanosensory transduction, playing a fundamental role in the perception of sensory stimuli within living organisms. These mechanosensitive ion channels, coalescing to form mechanical receptors, are strategically positioned within the sensory neuron terminals intricately nestled within the epidermal layer. Their cardinal function lies in the reception and transduction of mechanical signals emanating from the external milieu, ultimately orchestrating their conversion into electrical signals. This intricate process, duly recognized as mechano-electrical transduction (MET), engenders alterations in neuronal excitability, thereby facilitating the transmission of vital information to higher-order neural centers. ¹

Among these, NOMPC, a transient receptor potential (TRP) superfamily MS channel identified in C. elegans, Drosophila, and zebrafish, assumes a central role in various mechanosensing-related activities, encompassing soft touch perception, auditory function, and locomotion. ^2–6 A salient characteristic distinguishing NOMPC from its counterparts is its possession of a unique intracellular domain comprising twenty-nine ankyrin repeats, stacking one by one into a distinctive spring-like conformation. This exceptional architecture, which was proposed to be associated with microtubules within cellular contexts, designates NOMPC as the only known standalone tether-gating MS channel. ⁷ Consequently, in contrast to the extensively discussed force-from-lipid MS ion channels, ^8,9 an in-depth investigation into the specific molecular mechanism of NOMPC gating stands poised to offer a paradigmatic insight into the realm of the tethered MS ion channels.

In previous studies, it was proposed that compression of the intracellular ankyrin repeat (AR) domain or a pushing force exerted on the AR domain from the intracellular side, could open the channel. ^10,11 In the meantime, the TRP domain, widely recognized as crucial in the gating process of TRP channels, ¹² was observed to rotate clockwise (CW, looking from the intracellular side) under the torque generated by the pushing force. ¹¹ Therefore, it was still ambiguous whether the compressive force or the twisting torque was the direct mechanical cause behind channel gating. In this work, we determine the critical force component for the gating of NOMPC by using all-atom molecular dynamics (MD) simulations. Our results suggest that the twisting force is the key to the channel gating, while the compression-twist coupling of the AR domain ensures that a compression of the AR domain can activate the channel as well.

2 Results

2.1 A membrane-parallel torsion force was observed on the TRP domain during the push-to-open gating process

For a comprehensive exploration of the atomistic gating details through molecular dynamics (MD) simulations with an affordable computational cost, we employed a divide-and-conquer strategy as depicted in Fig. 1(a). Essentially, the membrane-protein simulation system was divided into two subsystems. System I was composed of the transmembrane (TM) region, intracellular linker helices (LH) domain, and ankyrin repeat 29 (AR29). Within the TM region, there exist six transmembrane helices (S1-S6) in each chain. Four S6 helices of each chain collectively constitute the central pore of NOMPC, and the TRP domain is sandwiched between the TM domain and the LH domain. System II included the LH domain and twenty-nine ankyrin repeat (AR) units. The four spring-like AR chains of the tetrameric protein were organized into a bundle structure, thought to be a key element in mechanical perception. ^13–15

Figure 1:

Building upon the precedent push-to-open model elucidated by Wang et al., ¹¹ we performed steered molecular dynamics (SMD) simulations for system I. In this context, a pushing/compressive force was applied to the ankyrin repeat 29 (AR29), pointing from the intracellular to the extracellular side (trajectories: CI-1/2). Consistent with our previous findings, the transmembrane central pore of NOMPC exhibited a discernible dilation, accompanied by the upward tilt and CW rotation of the TRP domain (Fig. S1). The upward tilt was probably caused by the membrane-normal component of the applied mechanical force, while the CW rotation was generated by the membrane-parallel component that led to a torque pointing to the extracellular side, raising the question that which of these two components was the direct driving force for the channel opening.

To address this pivotal question, we quantitatively analyzed the force distribution in the above simulations CI-1/2, so that we could decompose the force into the membrane-normal and membrane-parallel components to investigate which one could open the channel. In particular, we used the force distribution analysis (FDA) method ¹⁶ to quantify the interactions between pairwise residues along the force transmission pathway within system I (Fig. 1(a)). According to our simulation condition and trajectory analyses, the force transmission initiated from the AR29 domain, traversed through the LH domain, and then proceeded to influence the geometrically adjacent yet non-covalently bonded TRP domain, ultimately guiding the S6 domain and instigating the pore opening. Our analysis showed that the strongest interactions between these domains were predominantly located between a few pairs of residues at the interfaces (Fig. S2), including A1118-N1148, Y1109-D1142 on the AR29-LH interface, and K1244-E1571, D1236-R1581, L1256-S1577 on the LH-TRP interface. Their proximity in space made their Coulomb and Lennard-Jones interactions stronger upon the conformational changes. This feature facilitated the stability of hydrogen bonds, particularly the hydrogen bonds on K1244-E1571, D1236-R1581, and L1256-S1577, with occupancy of 57%, 71%, and 24%, respectively. This was consistent with the hydrogen bond pairs identified in previous studies when a pushing force was applied. ¹¹

To provide a clear depiction of the net mechanical force distribution stemming from the externally applied mechanical force, we deducted the FDA results of the force-free simulation trajectories (trajectory: FI-1/2) from the simulation trajectories with pushing/compressive force (trajectory: CI-1/2), and defined the residual force distribution as the net-FDA results (Fig. S3). Analyzing the net mechanical force distribution between the LH and TRP domains using the net-FDA method, we found that, in addition to the membrane-normal force component, the total net mechanical force transmitted from the LH domain to the TRP domain also exhibited a membrane-parallel torsion force component, generating a torque pointing towards the extracellular direction perpendicular to the membrane (Fig. 1(b)-(c)). Moreover, the transmitted force was primarily concentrated on the E1571 and R1581 residues of the TRP domain (Fig. 1(c)).

In brief, the FDA results at the interface between the transmembrane and intracellular domains elucidated a discernible pathway of force transfer, mostly through a few residues located at the AR-LH-TRP interface. Both membrane-normal and membrane-parallel force components were observed on the TRP domain, which provided the basis to further investigate which component is the key factor in channel gating. As previous studies showed that a CW rotation of the TRP domain may be related to the channel opening, it was of particular interest to study the effect of the torsion force on the channel structure.

2.2 The membrane-parallel torsion force can open NOMPC, but not the membrane-normal pushing force on the TRP domain

Motivated by the above net-FDA results, we aimed to separate the contributions of the membrane-normal and membrane-parallel force components to the gating process and pinpoint whether one or both of these components played an indispensable role in gating. To achieve this, we needed to determine the point of application, the magnitude, and the direction of the exerted mechanical force. As the TRP domain is sandwiched between the AR and TM domains and its conformational change was proposed to be directly related to the channel gating, it is a natural choice to select the TRP domain to apply the decomposed forces in the MD simulations. Notably, our FDA in the last section showed that the predominant forces propagated from the LH domain to the TRP domain were distributed on two key residues, E1571 and R1581 (Fig. 1(b)-(c)). Therefore, we devised simulations to individually apply membrane-normal and membrane-parallel forces on these two residues of the TRP domain, and the magnitude and direction of the exerted forces were equivalent to those derived from the aforementioned net-FDA results (Fig. 1(b)-(c), Fig. 2(a)-(b)). As the force transmission occurred mainly via side-chain interactions along the interface between the LH and TRP domains, we selected the side chains of the E1571 and R1581 as the point of application of the exerted forces.

Figure 2:

Our simulations revealed that the membrane-parallel torsion force alone was sufficient to open the gate in all three MD simulation repeats (trajectories: MPI-1/2/3, Fig. 2(c)), whereas the membrane-normal pushing component did not trigger pore opening within the simulation time (trajectories: MNI-1/2/3, Fig. 2(d)). According to the pore size evolution in the trajectories MPI-1/2/3 (Fig. 2(c)), it became evident that the most constricted region along the pore, the gate situated around the hydrophobic I1554, exhibited an initial radius measuring less than 1.0 Å, as calculated by the HOLE program. ¹⁷ This configuration represented a closed architecture. Under the influence of the membrane-parallel torsion force, the radius of lower constriction was significantly dilated to above 2.0 Å after 600 ns. This dilation signified that water molecules could smoothly pass through the channel (Fig. S4), a phenomenon often indicative of the first step of channel opening. This is highly similar to our previous observation when a pushing force was applied to the AR domain. ¹¹ In contrast, in the trajectories MN-1/2/3 with the membrane-normal force on the TRP domain, no pore dilation was observed (Fig. 2(d)). In fact, the channel appeared to become even more constricted at the upper constriction region (filter), as depicted in Fig. 2(d). This indicated that the membrane-normal force component on the TRP domain was probably not able to open the channel, at least no such sign was observed within our simulation time.

Thus, our investigations substantiated that a torque on the TRP domain, pointing to the extra-cellular side, was the key to driving the channel open. On the basis of the original push-to-open or compress-to-open model, here we provide further details on the force analysis and channel opening, arguing that it is the membrane-parallel torsion force that acts as the direct mechanical cause to open the channel. It should be noted though, this does not mean the push-to-open or compress-to-open model was wrong, as our simulations showed that a pushing force on the LH domain can generate a torsion force on the TRP domain, ¹¹ so does the compression of the AR domain. ¹⁰ Apparently, the unique structure of NOMPC ensures that a straightforward pushing force or compression of the AR spring can spontaneously generate a torsion force to open the channel. Therefore, the torsional mechanical properties of the AR domain became more intriguing.

2.3 The torsional mechanical properties of the AR domain and force transmission pathway

To study the torsion property of the AR spring, we performed all-atom steered molecular dynamics (SMD) simulations for system II, in which we enforced a gradual clockwise (CW) rotation to the lower end of the AR domain (AR1-AR8) at a rate of 0.1 ° ns⁻¹ along the membrane-normal axis of NOMPC while keeping the upper end (LH domain) position-restrained (Fig. 3(a), trajectory: RII). From the trajectory, we extracted structural snapshots at an interval of 2.5° during the twisting process and relaxed these structures for 100 ns with both ends position-restrained. This protocol was employed to further relax the twisted AR spring to improve the accuracy of the calculated torsional mechanical parameters (Fig. S5). Please refer to the Methods section for more details.

Figure 3:

With the above simulation data, we calculated the torsional coefficient of the AR bundle, using the equation:

where M represented the torque applied to both ends, c denoted the torsion coefficient of the AR domain bundle, and ϕ was the average torsion angle of AR1-8 with respect to its initial position. As shown by the green plots in Fig. 3(b), the fitted torsion coefficient turned out to be (2.40 ± 0.06) × 10³ kJ mol⁻¹ rad⁻¹. Such a torsion coefficient roughly corresponds to a similarlysized cylindrical volume of uniform polyethylene or rubber material. ¹⁸

Interestingly, we observed a gradual increase in the compressive deformation of the AR spring as the torsional angle increased in our simulations (Fig. 3(b)). This implies a coupling between the twist and compression of the AR domain. To quantitatively describe this compression-twist coupling property, we defined the compression-twist coupling coefficient , which represented the compressive deformation induced by a unit torsional deformation. For the AR8-AR29 of NOMPC, the fitted value of the compression-twist coupling coefficient k_tc was (1.67 ± 0.14) nm rad⁻¹ within 10 degree torsion. Indeed, it was also noticed that a torque could be generated during the application of compressive forces in previous investigations. ^10,11 In light of these findings, we believe that the unique structure of the AR and LH domains bolster NOMPC’s robustness in the sensation of external stimuli, i.e., the compression-twist coupling property of the AR bundle ensures that both the twist and compression of the intracellular domain can generate a torque to rotate the TRP domain, facilitating the opening of the pore.

To pinpoint the force transmission pathway along the AR spring while it is twisted, we utilized the same net-FDA method as above to analyze the force distribution during the twisting process of the trajectory RII. The net-FDA results of each (i + 1)^th AR unit exerted from the adjacent i^th AR unit unveiled a spiral force transmission route, rooted in the helical structure inherent to the AR spring (Fig. 4(a)). Strong pairwise interactions predominantly resided in bonded residue pairs, while only a small fraction was associated with non-bonded residue pairs (Fig. 4(b)). Clearly, these bonded residue pairs, which are integral to the spiral backbone of the AR spring, were critical for transferring external forces. The non-bonded residue pairs were primarily situated at the adjacent turns within the AR α-helices and the interfaces between AR units. To clearly illustrate the details, we have magnified the local structure of AR13-15 (Fig. S6). As depicted in Fig. 4(c)-left, a representative non-bonded net-FDA of the local region of AR13-15 was shown. The data around the diagonal line represented the pairwise residues located at the adjacent turns between α-helices, while the wing-shaped data positioned further away from the diagonal indicated pairwise residues situated at the interfaces between AR units. Apparently, although the bonded residue pair interactions are predominant, the non-bonded interactions, both within the same AR unit and across different AR units, are important for the force transmission too.

Figure 4:

Furthermore, we noticed that there is a correlation between the distribution of non-bonded net-FDA and the existence/strength of hydrogen bonds, as demonstrated by similar patterns in Fig. 4(c). For clarity, we divided these two 737× 737 matrices (ResID of AR9-AR29: 400-1136) into 21× 21 subgroups (grouped by AR unit). Then, we calculated the Pearson correlation coefficient (PCC) between the non-bonded net-FDA and the occupancy of hydrogen bonds of each subgroup, resulting in a 21× 21 PCC matrix (Fig. 4(d)). Relatively high PCCs in the sub-diagonal grids of Fig. 4(d) (inter-AR), and even higher PCCs in the diagonal grids (intra-AR), are observed (inset of Fig. 4(d)). This underscored the crucial role of the hydrogen bonding network in facilitating force transmission. Therefore, we concluded that the force transmission route primarily consisted of bonded residues that formed the spiral backbone of the AR chains, with supplementary support from the hydrogen bonding network formed by adjacent residues in the same or neighboring AR unit.

3 Discussion

The unique structure of NOMPC and its fascinating AR spring have sparked great curiosity. Building upon the push-to-open and compress-to-open mechanism, here we propose a twist-to-open model, as illustrated in Fig. 5. At the core of this model was the idea that a twisting force on the TRP domain is the key to the opening of the pore. This is not to say the compress-to-open or push-to-open models are wrong, but that a compression/pushing on the AR domain will lead to a twist force on the TRP domain, and this twist force is the direct cause of the pore opening. The spiral structure of each AR chain, as well as the hydrogen bonding network formed among spatially adjacent residues, enables the AR bundle to convert the compression deformation to a twist deformation, thus exerting a membrane-parallel torsional force on the LH domain. In addition, the structure around the LH domain can also convert a membrane-normal pushing force to a membrane-parallel torsional force, as shown in our previous study. ¹¹ Consequently, the compression of the AR domain, or a pushing force from the intracellular side, can be converted to a torque on the TRP domain, which points to the extracellular side. We believe this torque is the direct mechanical cause that drives the TRP domain to rotate clockwise (looking from the intracellular side) to open the channel.

Figure 5:

The highly conserved TRP domain serves as a structural characteristic of the TRP protein superfamily. In our simulations of NOMPC, we observed a clockwise (CW) rotation of the TRP domain during the gating process of the proteins. Interestingly, such rotation was also observed in other TRP channels. For example, the recently resolved cryo-EM structures of the thermo/ligand-gated mm-TRPM8 (PDB: 8E4P, 8E4N, 8E4M, 8E4L; ¹⁹), thermo/ligand-gated rat-TRPV1 (PDB: 7L2H, 7L2M, 7L2N, 7L2L, 7MZ5; ²⁰), and ligand-gated homo-TRPV4 (PDB: 6BBJ, 7AA5; ²¹) showed similar gating motion of the TRP domain. Based on these findings, we believe that the rotation of the TRP domain could be a common feature in the gating mechanism of these TRP proteins, and the unique structure and compression-twist coupling property of the intracellular domain of NOMPC ensure its robustness in the perception of multiple forms of mechanical stimuli: not only twist, but also compression.

Although we calculated the detailed interactions between residue pairs of NOMPC, it should be noted that, due to the limitations of the FDA method, our force analysis lacked the entropic force, which may also play a significant role in the force transmission along the AR spring. Nonetheless, this won’t affect our twist-to-open model, since the entropic effect was present in our MD simulations. Another limitation of the present study is that the simulation time is very short compared to the physiological timescale. However, we found that the mechanical stimuli are transferred at a fast speed along AR springs in simulations, at around 1.8 nm ps⁻¹. ¹¹ Therefore, the mechanical stimuli were nearly instant for our simulation time, which allowed us to determine the trend or the direction of the conformational changes, even though the observation of the complete gating process is still challenging. In this context, we think that the dilation of the pore under a torsion force on the TRP domain is a piece of strong evidence that twist is the key to the channel opening.

4 Methods

4.1 The all-atom simulation Systems

To set up the molecular dynamics (MD) simulations, we first prepared the simulation systems using the cryo-EM structure of NOMPC (PDB: 5VKQ). ⁷ CHARMM-GUI was used to generate the configuration and topology of the simulation systems. ²² The parameter files were created for the CHARMM36m force field. ^23,24

The system I included the TM domain, the LH domain, and the AR29 of NOMPC. The PPM server was used to reorient the NOMPC structure to ensure that the TM domain of NOMPC was well located in a lipid bilayer. ²⁵ The protein was embedded in a POPC bilayer. Then, we solvated the protein-membrane complex in a water box using the TIP3P water model. Counterions were added to neutralize the system, ensuring physiological ionic strength by adding Na⁺ and Cl⁻ ions to achieve a concentration of 0.15 M. The details of simulation system I are shown in Table 1.

Table 1:

The system II included the LH domain and the AR domain (AR1-AR29) of NOMPC. The protein structures were immersed in a rectangular box of TIP3 water molecules and 0.15 M NaCl. The details of simulation system II are shown in Table 1.

4.2 All-atom MD simulations

All of the MD simulations were performed with GROMACS 5.1.2. ²⁶ The REDUCE program in AMBER was used to add hydrogens to the original PDB files and determine the protonation state of the histidine residues. ^27,28

For system I, we used the steepest descent algorithm to achieve energy minimization. Then, following up with a two-stage equilibration, a 0.4 ns NVT equilibration simulation with harmonic restraint was applied to the protein molecule (force constants of 4000 kJ mol⁻¹ nm⁻² on the back-bone and 2000 kJ mol⁻¹ nm⁻² on the side chains), and a 20 ns NPT equilibration simulation with gradually decreased restraint (from 2000 to 100 kJ mol⁻¹ nm⁻² on the backbone and from 1000 to 50 kJ mol⁻¹ nm⁻² on the side chains). During the equilibration processes, planar restraints were used to keep the positions of lipid head groups along the membrane-normal direction.

For system II, we used the steepest descent algorithm to achieve energy minimization. Then, also followed by a two-stage equilibration, a 0.2 ns NVT equilibration simulation with harmonic restraint was applied to the protein molecule (force constants of 400 kJ mol⁻¹ nm⁻² on the backbone and 40 kJ mol⁻¹ nm⁻² on the side chains), and a 10 ns NPT equilibration simulation with harmonic restraint (force constants of 400 kJ mol⁻¹ nm⁻² on the backbone and 40 kJ mol⁻¹ nm⁻² on the side chains).

The simulation temperature of the system was set to 300 K. The time step was 2 fs. The cubic periodic boundary condition was used during the simulations, and the van der Waals interaction was switched off from 10 to 12 Å. The long-range electrostatic interactions were calculated using the particle mesh Ewald (PME) method. ²⁹

4.3 Steered MD simulations

Steered molecular dynamics (SMD) simulations were utilized to probe the mechanical gating mechanisms of the ion channels. SMD was performed by attaching a virtual spring to a specific group of atoms, and a constant velocity or constant force was applied to the spring to pull the molecule along a predefined reaction coordinate. The force exerted on the molecule was recorded throughout the simulation. To ensure reproducibility, multiple independent SMD simulations were performed in this research.

The SMD simulation detail of the compressive trajectories of system I (CI-1/2) and free trajectories of system I (FI-1/2) can be found in previous work. ¹¹

For trajectories MPI-1/2/3 (membrane-parallel force component on system I Repeat-1/2/3) and MNI-1/2/3 (membrane-normal force component on system I Repeat-1/2/3) of system I (Fig. 2(c)-(d)), we separately exerted membrane-parallel (MPI-1/2/3) or membrane-normal (MNI-1/2/3) external force on the side chain of E1571 and R1581, as shown in Table 2. At the first 200 ns, we applied an increasing force by using the constant-rate pulling method. Then we applied the calculated constant force by using the constant-force pulling method in the last 600 ns. To prevent the protein from rotation as a whole, we used the enforced rotation method, ³⁰ setting the rotation speed as zero for the S1-S4 helices in the transmembrane region of the protein. Distance restraints were applied between the acceptors and donors on the side chains of residue pairs 1564-1568 and 1567-1571, to maintain the hydrogen bonds and prevent structural damage in the α-helix under the applied mechanical force.

Table 2:

For trajectory RII (Rotation of system II), we performed a three-step simulation to better study the torsional mechanical properties of the AR domain, as shown in Table 3.

Table 3:

In the first step, we aimed to relax the structure by fixing the relative rotational angle between both ends along the rotational axis, which is the Z axis of system II. We choose to manipulate AR1-AR8 as a whole, considering the fact that the cryo-EM structure of AR1-AR7 was not resolved but modeled. Therefore, it was prudent to exert and transmit torque starting from the more reliably resolved structure, i.e., from AR8 upward, while preserving the integrity of the entire protein structure. With the help of the enforced rotation function of GROMACS, ³⁰ we used the pm rotational potential on the bottom terminal of the AR domain (AR1-AR8) and the up terminal (LH domain) with a zero angular velocity for 20 ns. The pm rotational potential can constrain the group of atoms to rotate within the rotational plane as required, while allowing free movement in the direction perpendicular to the rotational plane. This means that the entire AR bundle will be subjected to torsion, while being able to stretch and compress freely.

In the second step, we clockwise rotated the AR1-8 along the Z axis at an angular velocity of 0.1 ° ns⁻¹, while keeping the LH domain still (Fig. 3(a)) with pm rotational potential. The difference in angular velocity between the bottom and up terminal accumulated a specific rotation angle within 10° (Fig. S5, trajectory: RII-0).

In the third step, when rotating at every 2.5° in trajectory RII-0, we extracted structural snapshots and held both ends of the AR domain for 100 ns to get a better equilibrated structure (Fig. S5, trajectory: RII-1 for the 2.5°-rotational angle, RII-2 for 5°, RII-3 for 7.5°, and RII-4 for 10°), and the readouts of torques on both ends gradually increased with opposite directions. Specifically, we used the trajectory RII-1/2/3/4 to calculate the torsional mechanical properties of the AR bundle and the trajectory RII-0 to study the force transmission pathway along the AR chain. When the torsional angles of the AR bundle are 2.5i° (i = 1, 2, 3, 4), we analyze the 100 ns-long RII-i trajectories. For the torques at both ends of the AR, we collect torque data on the LH domain every 2 ps. Under equilibrium conditions, the torques at the top and bottom of the AR bundle are equal in magnitude but opposite in direction. To assess the length change of the AR bundle, we measure the distance between the centroids of the LH domain and AR8 unit.

4.4 Force distribution analysis

Force distribution analysis (FDA) focuses on computational efficiency and low-memory usage and makes possible time series analysis of pairwise forces variation in long MD simulations and for large molecular systems. ¹⁶ It can compute atomic pairwise forces for all types of bonded interactions as well as Coulomb and Lennard-Jones potential by rerunning the trajectories. By summing up all the atomic pairwise forces within the two residues of interest, it can also calculate the pairwise force of the interaction between residues r_i and r_j, which is labelled as . Similarly, if we focus on the interaction between two domains, we can sum up all the pairwise forces of the residues from these two domains.

In the calculation, we made an input file which contained all the parameters that FDA needed to run. We set onepair = detailed and type = all so that all pairwise forces of different types between the same pair of residues are stored separately. The option residuebased was set as pairwise_forces_vector so that the residue-based pairwise forces were written to output files as vector values. The forces for every frame were written by setting time_averages_period = 0.

Then, forces sampled by MD simulations under two conditions (with and without perturbation) were compared. We analyzed the interactions between residues in the trajectories with the external force (CI-1/2), labelled as , and without the external force (FI-1/2), labelled as , respectively. Then the average net force distribution induced by the compressive external force, which was named net-FDA, , was obtained by subtracting the FDA result of trajectories without the external force from the FDA result of trajectories with the external force,

In the end, we averaged the vector net-FDA results of the four subunits by rotationally symmetrized around the protein center.

5 Competing interests

The authors declare no competing interests.

Acknowledgements

We thank Dr. Yang Wang at Wenzhou Institute, University of Chinese Academy of Sciences, for his help in setting up the initial simulation systems and for his valuable discussions. This work was supported by the General Program of the National Natural Science Foundation of China (Grant No. 32071251), and the Science Fund for Creative Research Groups of the National Natural Science Foundation of China (Grant No. T2321001).

Supplementary Information

Supplementary Figures

Figure S1:

Figure S2:

Figure S3:

Figure S4:

Figure S5:

Figure S6: