extfieldsofstudy: Medicine Biology

s2fieldsofstudy: Biology Chemistry

text:
Regulation mechanisms of the dual ATPase in KaiC

Significance KaiC, a core clock protein in the cyanobacterial circadian clock system, hydrolyzes adenosine triphosphate (ATP) at two distinct sites in a slow but ordered manner to measure the circadian timescale. We used biochemical and structural biology techniques to characterize the properties and interplay of dual-adenosine triphosphatase (ATPase) active sites. Our results show that the N-terminal and C-terminal ATPases communicate with each other through an interface between the N-terminal and C-terminal domains in KaiC. The dual-ATPase sites are regulated rhythmically in a concerted or opposing manner dependent on the phase of the circadian clock system, controlling the affinities of KaiC for other clock proteins, KaiA and KaiB.

of KaiC ATPase, the mechanisms of its activation and inactivation remain unknown. In this study, we measured the contributions of C1-ATPase and C2-ATPase separately by introducing amino acid substitutions in KaiC. Then, we attempted to estimate and parameterize the degree to which C1-ATPase is activated and inactivated by the C2 domain. Furthermore, we performed a crystal structure analysis of the activated and inactivated forms of KaiC. Our results suggest that the components (C) Measured rates of ADP production. Gray, white, and hatched bars represent the ATPase activities of full-length KaiCs (contributions from C1 and C2), full-length KaiCs carrying the E318Q substitution (contribution from C1), and truncated KaiCs (KaiC1) consisting solely of the C1 domain (C1-ATPase basal ; C1 basal activity without the C2 domains). ATPase activities were assayed by measuring ADP production for 10 to 16 h. In the case of KaiC-ST (0.20 mg/mL) + KaiA (0.12 mg/mL), data are given as apparent activities estimated by the ADP production for 10 h. The ATPase activities were analyzed using one-way ANOVA (F 4,28 = 2.714, P = 2.30 × 10 À17 ) and the Bonferroni post hoc t test. *Insignificant (P > 0.05/10); **significant (P < 0.005/10); ***significant (P < 0.0005/10). (D) Calculated rates of ADP production in the C1 (C1-ATPase full ; red bars) and C2 domains (C2-ATPase full ; blue bars) of the full-length KaiCs. C1-ATPase full indicates ATPase activity of fulllength KaiCs carrying the E318Q substitution (white bars in C). C2-ATPase full indicates ATPase activity of full-length KaiCs (gray bars in C) minus C1-ATPase full (white bars in C). *Insignificant (P > 0.05): t test of C2-ATPase full . (E) The b values defined as the ratio of C1-ATPase full (red bars in D) to C1-ATPase basal (hatched bars in C) (details are in Eq. 1 and Materials and Methods). The b values were analyzed using one-way ANOVA (F 3,26 = 2.975, P = 3.95 × 10 À17 ) and the Bonferroni post hoc t test. *Insignificant (P > 0.05/6); **significant (P < 0.005/6); ***significant (P < 0.0005/6). (F) Potential correlation between b value and KaiB-KaiC affinity. The fraction of KaiB forming KaiB-KaiC complexes was determined by native PAGE analysis (SI Appendix, Fig. S3A). C1/C2-ATPase and b values for a series of mutants (black circles) are compiled in SI Appendix, Table S3. of the dual ATPase communicate with each other through C1-C2 interaction to control the assembly and disassembly cycle of KaiA and KaiB.

Results
Contributions of C1-and C2-ATPases. The steady-state C1/C2-ATPase activity of fully dephosphorylated KaiC (KaiC-ST) was 13.6 ± 1.4 d À1 at 30°C (gray bars in Fig. 1C). To estimate the contribution of C1-ATPase in full-length KaiC (C1-AT Pase full ), we inactivated the kinase and C2-ATPase by replacing the catalytic glutamate 318 with glutamine (13)(14)(15)(16)22). The resultant KaiC E318Q -ST exhibited a measured activity of 15.8 ± 0.6 d À1 (white bars in Fig. 1C), which was comparable with the activity of KaiC-ST. Given that this measured activity of KaiC E318Q -ST reflects the C1-ATPase full activity of KaiC-ST (red bars in Fig. 1D), the contribution of C2-ATPase full can be estimated as the difference in the measured activities between KaiC-ST and KaiC E318Q -ST, which was as small as À2.2 ± 1.6 d À1 (blue bars in Fig. 1D). On the other hand, the ATPase activity of a truncated KaiC (KaiC1) consisting solely of the C1 domain (11), which reflects the basal activity of the C1-ATPase (C1-ATPase basal ) without the influence of the C2 domain, was 10.4 ± 0.5 d À1 (hatched bars in Fig. 1C).
The regulatory effect of the C2 domain on C1-ATPase full was evaluated using a parameter (b value) defined as the ratio of C1-ATPase full to C1-ATPase basal by the following equation: where C1-ATPase full and C1-ATPase basal correspond to the ATPase activities of full-length KaiC with the E318Q substitution and truncated KaiC consisting solely of the C1 domain, respectively. The b value was 1.52 ± 0.10 for KaiC-ST (Fig.  1E), reflecting the up-regulation of its C1-ATPase by the C2 domain. These results suggest that the C1-ATPase full activity in KaiC-ST is ∼50% up-regulated by the C2 domain (Fig. 1E), while the contribution of the C2-ATPase full in KaiC-ST is negligibly small (Fig. 1D). It must be noted that we did not rely on an alternative strategy of inactivating C1-ATPase by corresponding E77Q/E78Q mutations (13)(14)(15)(16) to estimate the contribution of C2-ATPase, as KaiC E77Q,E78Q -ST and even KaiC1 E77Q,E78Q exhibited residual C1-ATPase activity as high as 3.4 ± 1.4 and 1.9 ± 0.6 d À1 , respectively ( Fig. 1C and SI Appendix, Supplementary Text).
We found that KaiC-AA is rather similar to a KaiAstimulated state of KaiC-ST. KaiA was previously shown to activate the ATPase activity of KaiC (12). Consistent with that study, the addition of KaiA to KaiC-ST resulted in the activation of the dual ATPase. The C1-ATPase full activity of KaiAstimulated KaiC-ST was 18.2 ± 1.7 d À1 (Fig. 1D), and it was ∼80% up-regulated relative to KaiC-ST alone (t test: **significant, P < 0.005/10 in Fig. 1C and b = 1.75 ± 0.18 in Fig.  1E). The contribution of C2-ATPase full in KaiC-ST became very obvious in the presence of KaiA, with enhanced activity of 15.2 ± 1.7 d À1 for the initial 10 h (Fig. 1D). These observations suggest that the KaiA stimulation of C1/C2-ATPase of KaiC-ST is mimicked in KaiC-AA even without KaiA.  Table S2). Dark-colored circles correspond to the mean W1 positions from independent C1-C1 interfaces (pale-colored circles).
Structural Basis for Up-and Down-Regulation of C1-ATPase full .
ATP hydrolysis includes a reaction step in which a lytic water molecule (W1) attacks the phosphorus atom (P γ ) of the terminal γ-phosphate of ATP. Thus, the relative positioning of W1 and P γ is one of the crucial factors determining the reaction efficiency, as reported previously (11). To investigate the W1 position, we sought to crystallize KaiC-AA (SI Appendix, Table  S1) and to compare the results with those for KaiC-ST and KaiC-pSpT, which we previously reported (23). We could identify a reasonable electron density for W1 in all cases ( Fig. 2 A-C and SI Appendix, Table S2), except for KaiC-pSpT. Fortunately, however, another 2.2-Å crystal structure of thermophilic KaiC from Thermosynechococcus elongatus (TeKaiC) was determined in the TeKaiC-pSpT state (SI Appendix, Table S1), with a clear electron density for W1.
In TeKaiC-pSpT ( Fig. 2A), W1 was hydrogen bonded to the carbonyl oxygen atom of F200 and R227 N η (F199 and R226 in KaiC) and sequestered at the position most unfavorable for S N 2-type ATP hydrolysis, with a W1-P γ distance of 5.1 Å and a W1-P γ -O 3β angle of 141° (Fig. 2D). By contrast, W1 was located at a less unfavorable position in KaiC-ST (Fig. 2B), with a shorter W1-P γ distance (4.0 Å) and a larger W1-P γ -O 3β angle (147°) (Fig. 2D). W1 was found in positions more suitable for ATP hydrolysis in KaiC-AA ( Fig. 2 C and D). The descending and ascending orders of the W1-P γ distance and the W1-P γ -O 3β angle (Fig. 2D), respectively, matched the increasing order of the C1-ATPase full activity (KaiC-EE[pSpT] < KaiC-ST < KaiC-AA) ( Fig. 1 C and D). These results confirm that the relative positioning of W1 and P γ is one of the factors impacting up-and down-regulation of C1-ATPase full .
Structural Basis for Up-and Down-Regulation of C2-ATPase full .
The electron densities for W1s involved in C2-ATPase could not be identified in the crystal structures of KaiC-pSpT, TeKaiC-pSpT, KaiC-ST, or KaiC-AA, possibly due to their dynamic nature. On the other hand, we found unique variation patterns in a tunnel-shaped void in the C2-C2 interface (Fig.  3), through which bulk water molecules migrate to the C2-ATPase active site.
In KaiC-ST, the distal side of the γ-phosphate of C2-ATP is stereochemically crowded by several residues (E318, S431, and T432) and consequently, does not have enough space to stably accommodate any water molecules ( Fig. 3 A, Inset). Although a tunnel-shaped void was confirmed near the distal side, it was undulated ( Fig. 3A) and too narrow (Fig. 3D) to freely take up water molecules (radius of ∼1.4 Å) from the inner-radius side of the hexamer. The corresponding water tunnels were not uniform within the KaiC-pSpT hexamer (SI Appendix, Fig. S1), but the shortest one suggests more effective uptake by KaiC-pSpT than by KaiC-ST ( Fig. 3 B and D). By contrast, the distal side of KaiC-AA was more relaxed. A wider straight channel  Fig. S2) are shown by green, blue, magenta, and gray lines, respectively. Note that the diversity of tunnel shapes at each protomer interface is more pronounced in KaiC-pSpT (SI Appendix, Fig. S1). Open circles correspond to the point where the tunnel radius first falls below the radius of a water molecule (1.4 Å; dashed line in D). The traces of KaiC-AA and the A431S/A432T model of KaiC-AA are terminated at the path length of ∼10 Å because their tunnels are shorter by ∼5 Å than those in other KaiCs. In KaiC-ST, the first half of the tunnel is undulated (thick black line in A, Inset) and then narrowed by bottleneck residues down to a radius of ∼0.8 Å, far smaller than the radius of the water molecule. Thus, the tunnel of KaiC-ST is not wide enough to freely take up water molecules from the inner-radius side. By contrast, the tunnel without a bottleneck that exhibits straight extension in KaiC-AA is ready to transfer water molecules from the inner-radius side to the active site of C2-ATPase. The modeling of A431S/A432T using the KaiC-AA template results in a narrowing of the tunnel radius at the path length of around 3 Å, but neither the total length of the tunnel nor the radius near the origin are affected at all.
( Fig. 3C) connected the space in the distal side to the exterior of the hexamer and appeared to allow effective migration of water molecules into the active site of C2-ATPase (Fig. 3D). Furthermore, the side chain of catalytic E318 was positioned at one end of the water tunnel and could immediately activate the migrated water molecules (Fig. 3 C, Inset). The descending and ascending orders of the path length and width of the tunnel-shaped voids, respectively, approximately matched the increasing order of the C2-ATPase full activity (KaiC-ST < KaiC-EE[pSpT] < KaiC-AA) ( Fig. 1 C and D). These observations provide a structural explanation for why C2-ATPase is down-regulated in KaiC-ST but up-regulated in KaiC-EE and KaiC-AA.
To evaluate local effects of the S431A/T432A substitutions that intrinsically reduce the volume of the side chains, we constructed an A431S/A432T model using the main chain structure of KaiC-AA as a template (SI Appendix, Fig. S2). While a partial decrease in the tunnel radius at the path length of around 3 Å was observed (gray line in Fig. 3D), neither the total length of the tunnel nor the radius at the origin were affected at all, and the position where the tunnel radius first became below the radius of the water molecule was sufficiently closer to the origin than that of KaiC-ST and KaiC-pSpT (open circles in Fig. 3D). The efficient uptake of the water molecules into C2-ATP in KaiC-AA can be attributed to the local reduction of the side-chain volume by the dual alanine substitutions and to other nonlocal structural changes that shorten the tunnel length and enlarge the radius around the origin. Quaternary Structural Changes. A region upstream of the dual-phosphorylation sites in the C2 domain called phosphor switch (PSw) (23) adopts coiled and helical conformations in KaiC-pSpT (EE) and KaiC-ST, respectively (Fig. 4A). Our previous study (23) suggested that this fold switch of PSw influences the quaternary structure of the C1 ring through the C1-C2 interface. When superimposed using the C1 domains of KaiC-pSpT and KaiC-ST as references (Fig. 4C), one of the two adjacent C1 domains (right subunit in Fig. 4C) of KaiC-pSpT shifted by ∼1 Å away from the C2 domain, but conversely, the other C1 domain (left subunit in Fig. 4C) shifted toward the C2 domain. These observations indicate that during the helix to coil transition of PSw from KaiC-ST to KaiC-pSpT, every C1 domain gently tilts and then leans over a neighboring C1 domain (Fig. 4C). This tilting and sliding movement causes W1, which is sequestered by F199 and R226 (F200 and R227 in TeKaiC), to move one position further toward the ATP molecule bound to the C1 domain (Figs. 2D and 4E).
PSw of KaiC-AA was found to adopt a loop conformation (Fig. 4A), providing the tunnel-shaped void in the C2-C2 interface to allow uptake of bulk water molecules to the C2-ATPase active site (Fig. 3 C and D). Compared with KaiC-ST, the coiled PSw of KaiC-AA was shifted away from the nearest C2-ATP, similar to what was observed in KaiC-pSpT (Fig. 4A). However, the effects of these PSw shifts on two adjacent α-helices (α 8 and α 7 ) were slightly but notably different between KaiC-AA and KaiC-pSpT. In KaiC-AA, the short helix (α 9 ) located downstream of PSw was shifted to the outerradius side of the hexamer, and the resultant steric hindrance between I433 (α 9 ) and I397 (α 8 ) pushed the N terminus of α 8 toward the C1 direction by ∼1 Å (Fig. 4B). Concomitantly, α 7 placed in parallel to α 8 was pushed in the same direction. The side chain of E357 (α 7 ) was hydrogen bonded to P237 and R218 in the C1 domain, and it was inserted into the C1 side as α 7 changed its orientation. These positional shifts of PSw, α 9 , α 8 , and α 7 in KaiC-AA assist the movement of every C1 domain to the outer-radius side by ∼1 Å (Fig. 4D), finally shifting W1 in every C1-C1 interface into a position more suitable for ATP hydrolysis (Fig. 2D) and thus, further activating C1-ATPase (t test: ***significant, P < 0.0005/6; **significant, P < 0.005/6 in Fig. 1E).
Coiled and shifted PSw also influenced the C-terminal region of KaiC-AA. Amino acids positioned from 488 to 497 form a loop (A loop) (24) in both KaiC-ST and KaiC-pSpT, but the loop is deformed in KaiC-AA to allow its KaiA-binding C terminus (498 to 519) (25) to be more exposed to solvent (Fig. 4A). This observation is consistent with those of previous studies (26,27) that reported a high affinity of KaiA to KaiC-AA, implying that the structure of KaiC-AA reflects the KaiAstimulated or KaiA-postbound state of KaiC-ST. Consistent with the expectation that the A loop is pulled down by KaiA binding to the C terminus of KaiC-ST (24), a key hydrogen bond between G492 and E444 is disrupted in KaiC-AA (Fig. 4A).
Relationship between the Regulatory State of C1-ATPase full and KaiB Affinity. It has been shown that KaiB binds to the C1 domain of KaiC (28)(29)(30), and the affinity and kinetics of the KaiB-KaiC interaction have been investigated in detail (31,32). To evaluate the functional significance of the observed quaternary structural changes of the C1 ring (Fig. 4 C-E), we studied the affinity of KaiC to KaiB with native polyacrylamide gel electrophoresis (PAGE) analyses. A mixture of equimolar amounts of KaiB and KaiC-ST (5 μM on a monomer basis) after incubation at 30°C for 30 h revealed that 50% of KaiB bound to KaiC-ST (>97% dephosphorylation) to form a KaiB-KaiC-ST binary complex (SI Appendix, Fig. S3A). Interestingly, as summarized in Fig. 1F, the KaiB affinity of a series of KaiC mutants revealed a correlation with the b value; the affinity of KaiC for KaiB decreased as the b value increased (SI Appendix, Fig. S3A and Table S3). This correlation suggests that the C1-ring status responsible for KaiB binding switches in the range of b values between ∼1.3 and 1.7. While the b value of KaiC-ST was in the transition region, the affinity of KaiC-AA for KaiB was below the present limit of detection by native PAGE, as reported previously (6,33,34). Astonishingly, as long as the rules for up-and down-regulation are satisfied, KaiC could dramatically switch even from extremely lowaffinity (KaiC-AA) to high-affinity forms (KaiC E357Q -AA). In this context, KaiC-AA is a poor mimic of KaiC-ST in terms of not only structure (Figs. 2 and 3) but also, KaiB binding (Fig.  1F). Our results clearly indicate that the C1-ATPase full of KaiC-ST is not up-regulated enough to completely disrupt the interaction between KaiB and KaiC-ST (Fig. 1F).
The complete disassembly of KaiB from KaiC requires the KaiA-assisted activation of the dual ATPase. Immediately after the addition of KaiA (t = 0 h) to the KaiB-KaiC-ST complex that was preformed for 48 h (black line in Fig. 5A), KaiA was trapped on KaiB interacting with the C1 domain of KaiC-ST and then eluted from a size-exclusion chromatography (SEC) column as the KaiA-KaiB-KaiC-ST ternary complex (t = 0) (Fig. 5A). Within 1.5 h at 30°C, the peak intensity of the ternary complex decreased, whereas that of free KaiB increased (Fig. 5 A, Inset), and the dual-phosphorylation site of KaiC-ST remained unchanged (Fig. 5C). By contrast, more ADP was produced by the KaiA/KaiB/KaiC-ST mixture than by the KaiB/KaiC-ST mixture within 1.5 h, peaking ∼4 h after the KaiA addition (Fig. 5D). A small amount of free KaiA (Fig. 5E and SI Appendix, Fig. S3B) likely binds the C-terminal loop of KaiC-ST and then activates ATPase to disrupt the interaction between KaiB and KaiC-ST (31), as a KaiC mutant (KaiC Δ505 ) lacking the C terminus (506 to 519) (19, 24, 33) was assembled into a KaiA-KaiB-KaiC Δ505 -ST complex that rarely disassembled (Fig. 5B). These results are qualitatively consistent with the earlier observation (35) that KaiB and KaiC-EE are prevented from forming the complex in the presence of an excessive amount of KaiA. Notably, the KaiA-assisted disassembly shown in Fig. 5 has the appearance of an autocatalytic reaction often found in oscillators; a KaiA molecule released from a ternary complex can react with another ternary complex to promote the release of another KaiA molecule (SI Appendix, Fig. S4F).
To examine this hypothesis, we conducted a preliminary mathematical model analysis of the KaiA-assisted disassembly. Our present model provides only a qualitative assessment because it does not take into consideration KaiB rebinding to the KaiC hexamer (SI Appendix, Extended Methods). Nevertheless, important aspects of our experimental observations were reflected by the model when optimized parameter values were used. First, the disassembly of the KaiA-KaiB-KaiC-ST complex preceded the accumulation of phosphorylated KaiC (Fig. 5  A and C and SI Appendix, Fig. S4A). Second, free KaiB accumulated ∼4 h before free KaiA increased (Fig. 5E and SI  Appendix, Fig. S4B). Third, the concentration of the KaiA-KaiC complex peaked at ∼4 h (SI Appendix, Fig. S4A), when the time-dependent increase in KaiC-SpT abundance was maximized (Fig. 5C). Fourth, the present model (SI Appendix, Fig. S4D) successfully predicted the KaiA-dependent acceleration of disassembly and phosphorylation dynamics (SI Appendix, Fig. S4E). Fifth, the KaiA-KaiB-KaiC Δ505 -ST complex, which was rarely disassembled and phosphorylated (Fig.  5B), could be reproduced when the binding rates of KaiA to the C2 domain of KaiC Δ505 -ST were 30-to 50-fold lower than those to the C2 domain of KaiC-ST (SI Appendix, Fig. S4C). These findings qualitatively support the KaiA-assisted and autocatalytic disassembly of the KaiA-KaiB-KaiC-ST complex through the activation of C1/C2-ATPase.
Although our results may emphasize a unidirectional causality from the C2 domain to the C1 domain, we also obtained several clues for an allosteric reversal from the C1 domain to the C2 domain. As listed in SI Appendix, Table S3, the R217A substitution in KaiC1 little affected the C1-ATPase basal activity (9.0 ± 0.6 d À1 with R217A and 10.4 ± 0.5 d À1 without R217A). However, in the full-length KaiC, the effect on the C2-ATPase full activity became pronounced (8.5 ± 1.9 d À1 with R217A and À2.2 ± 1.6 d À1 without R217A). The tendency of the allosteric reversal is further pronounced in the C2-AT Pase full activity under the KaiC-AA background (195 ± 25.7 d À1 with R217A and 5.0 ± 2.9 d À1 without R217A). It is clear that there exists a bidirectional communication between the C1 and C2 domains.
In this context, it would be important to mention the effect of the E318Q mutation, which was systematically introduced to estimate the contribution of the dual-ATPase activities. Selective suppression of C2-ATPase by the E318Q mutation would result in a longer-lived C2 domain in the ATP-bound state. The analysis to separate the C1/C2 contributions is based on the assumption that such a C2-ATP-bound state affects neither the structure nor activity of the C1 domain. While our experimental data indicate that the E318Q substitution has little effect on the overall C1/C2-ATPase activity (Fig. 1C), the long-lived C2-ATP-bound state could have some effects on the structure of the C1 domain.
The present results showed that the properties of KaiC-AA were more similar to those of KaiC-ST stimulated by KaiA than to those of unstimulated KaiC-ST. Even in the absence of KaiA, C1-ATPase full and C2-ATPase full of KaiC-AA were both activated to a similar level as when KaiA was added to KaiC-ST (Fig. 1D). In contrast to KaiC-ST, PSw of KaiC-AA was coiled (Fig. 4A) to form a tunnel-shaped cavity at the C2-C2 interface for bulk water uptake (Fig. 3C), and at the same time, W1 was relocated to a suitable position for ATP hydrolysis (Fig.  2D) by the quaternary structural change of the C1 ring (Fig.  4E). These results suggest that KaiC-AA retains some of the structural and functional features of KaiC-ST stimulated by KaiA.
Assuming that KaiC-AA can be treated as a functional ( Fig.  1 C and D) and structural (Figs. 2D and 3D) mimic of KaiC-ST stimulated by KaiA, we can infer the results of KaiC-ST activation triggered by KaiA by tracing the conformational change from KaiC-ST to KaiC-AA starting from the KaiA binding site (Fig. 6). The folded A loop hooks the KaiAbinding C terminus (498 to 519) of KaiC-ST, which is in a less exposed configuration (Fig. 6A), by a hydrogen bond formed between the main-chain carbonyl oxygen atom of G492 and the main-chain amide of E444 (Fig. 4A). This sequence of events begins with KaiA binding to the C-terminal tail of KaiC-ST (25) and drawing out the A loop located on the N-terminal side of the tail (24) (Fig. 6B). This induces breaking of the hydrogen bond between G492 and E444, as confirmed in KaiC-AA (Fig. 4A). The destabilization of the A loop causes the helix to coil transition of PSw adjacent to the A loop (Fig. 6B). Eventually, the tunnel-shaped void at the C2-C2 interface is expanded, and water molecules are taken up through it to the C2-ATPase active site where E318 and C2-ATP are located (Fig. 3). The effect of the structural change in PSw is not confined to the C2 ring but extends to the C1 ring through the C1-C2 interface. The short α 9 -helix downstream of PSw moves along its helical axis to the outer-radius side of the hexamer, pushing the α 8 -and α 7 -helices up to the C1 side (Fig. 6C). These perturbations cause a change in the quaternary structure of the C1 domain so that the ring structure expands slightly (Fig. 6D); W1 that reacts with C1-ATP is sequestered at the C1-C1 interface via hydrogen bonds between F199 and R226 (F200 and R227 in TeKaiC) and is relocated to a position much more suitable for ATP hydrolysis (C1-ATPase) by the change in the C1-C1 interface.
By contrast, the mechanism of inactivating C1-ATPase will be inferred by tracing the conformational change from KaiC-ST to KaiC-pSpT or to KaiC-pST (23). The helix to coil transition of PSw from KaiC-ST to KaiC-pSpT is coupled to the rearrangement of polar residues constituting electrostatic interactions in the C1-C2 interface (23), which causes the tilt and slide movement of the C1 domains to facilitate the quaternary structural change of the C1 ring (Fig. 4C). This update of the C1-C1 interface results in W1 sequestration at a position much less suitable for C1-ATP hydrolysis (Fig. 4E).
At the same time, it is necessary to mention the possible side effects of the S431A/T432A substitutions. While the changes in the tertiary and quaternary structures (Figs. 3 and 4), as well as the change in the arrangement of W1 (Fig. 2), are intrinsic observations, we cannot rule out the possibility that the expansion of the local cavity around A431/A432 in KaiC-AA is an artificial volume reduction due to the dual alanine substitutions (Fig. 3D). Although not limited to this study, amino acid variants designed from a limited number of 19 choices do not always show ideal behaviors reflecting only the intention of designers and are sometimes accompanied by unexpected side effects due to amino acid substitutions. The validity of our proposal and the need for its update must ultimately be tested by elucidating the complex structure of KaiC-ST and KaiA.
Our results indicate that the activation of dual-KaiC ATPases contributes to a dawn-phase autocatalytic disassembly of the ternary night complexes (Fig. 6 and SI Appendix, Fig.  S4F) and is crucial for resetting subjective night and then pushing the system forward along the unidirectional cycle. Although the structural details of C2-ATPase are still unclear because we were unable to identify the W1s for C2-ATPase, it is simply amazing that it possesses a dynamic range as wide as 0 to 200 d À1 (Fig. 1D and SI Appendix, Table S3). Given that the activity level of C2-ATPase is suppressed to zero in KaiC-ST, its autoinhibitory mechanism and its physiological significance are the next important research targets, which will hopefully lead to a better understanding of the structural and functional roles of the KaiC dual ATPase in the circadian clock system.
All crystals were obtained by the vapor diffusion method as reported previously (23). The TeKaiC-pSpT crystal in the P2 1 2 1 2 1 space group was obtained in solutions containing 80 mM acetic acid and 1.0 M sodium acetate and frozen with 25% (wt/vol) polyethylene glycol 8000. The KaiC-AA crystal in the R3 space group was obtained in solutions containing 100 mM tris(hydroxymethyl)aminomethane (Tris)ÁHCl (pH 7.0), 1 M KCl, 0.8 M sodium/potassium tartrate, 0.95 M sodium acetate, and 5 mM adenylyl-imidodiphosphate.
Data Collection and Structure Determination. Preliminary X-ray diffraction experiments were conducted using FRX-Synergy (RIGAKU), and the final dataset was collected on beamline BL44XU at SPring-8 (Harima, Japan) at 100 K under a cryostream. Diffraction images were recorded with an MX-300HE (Rayonix) or PILATUS (Eiger) detector and processed using the HKL2000 (37) and XDSGUI (38) software. Initial phases were determined by molecular replacement methods using structures previously deposited as 2GBL (39) or 4O0M (40) in Protein Data Bank (PDB) and MOLREP (41). Refmac5 (42) and COOT (43) were used to refine and model the crystal structures, respectively. Graphic representations were drawn using PyMOL (Schr€ odinger). Tunnel-shaped voids were analyzed with a probe radius of 0.6 Å using CAVER 3.0 (44). Statistics of data collection and refinements are listed in SI Appendix, Table S1.
Biochemical Assays of Kai Proteins. The four phosphorylation states of KaiC were detected by sodium dodecyl sulfate-polyacrylamide gel electrophoresis and then quantified by LOUPE software (45). The number of ATP molecules hydrolyzed into ADP molecules per KaiC monomer per unit time was quantified at 30°C as previously described (18,19). The ATPase activities were analyzed using one-way ANOVA and the Bonferroni post hoc t test.
Determination of b Values as a Measure of Up-and Down-Regulation of C1-ATPase full . In this study, the ATPase activity of a truncated KaiC consisting solely of the C1 domain (11) was defined as the basal activity (C1-ATPase basal ) of C1-ATPase. The ATPase activities of full-length KaiC and its mutants, each of which carries the E318Q substitution that inactivates both kinase and C2-ATPase, were measured to estimate the contribution of C1-ATPase (C1-ATPase full ) (SI Appendix, Table S3). The b value was then calculated as the ratio of C1-ATPase full to C1-ATPase basal using Eq. 1. The b values were analyzed using one-way ANOVA and the Bonferroni post hoc t test.
Data Availability. Atomic coordinates and structure factors have been deposited in PDB (ID codes 7DY1 [TeKaiC] and 7DYE [KaiC-AA]). All other data are included in the manuscript and/or SI Appendix.


Upregulation of CCNB2 and Its Perspective Mechanisms in Cerebral Ischemic Stroke and All Subtypes of Lung Cancer: A Comprehensive Study

Cyclin B2 (CCNB2) belongs to type B cell cycle family protein, which is located on chromosome 15q22, and it binds to cyclin-dependent kinases (CDKs) to regulate their activities. In this study, 103 high-throughput datasets related to all subtypes of lung cancer (LC) and cerebral ischemic stroke (CIS) with the data of CCNB2 expression were collected. The analysis of standard mean deviation (SMD) and summary receiver operating characteristic (SROC) reflecting expression status demonstrated significant up-regulation of CCNB2 in LC and CIS (Lung adenocarcinoma: SMD = 1.40, 95%CI [0.98–1.83], SROC = 0.92, 95%CI [0.89–0.94]. Lung squamous cell carcinoma: SMD = 2.56, 95%CI [1.64–3.48]. SROC = 0.97, 95%CI [0.95–0.98]. Lung small cell carcinoma: SMD = 3.01, 95%CI [2.01–4.01]. SROC = 0.98, 95%CI [0.97–0.99]. CIS: SMD = 0.29, 95%CI [0.05–0.53], SROC = 0.68, 95%CI [0.63–0.71]). Simultaneously, protein-protein interaction (PPI) analysis indicated that CCNB2 is the hub molecule of crossed high-expressed genes in CIS and LC. Through Multiscale embedded gene co-expression network analysis (MEGENA), a gene module of CIS including 76 genes was obtained and function enrichment analysis of the CCNB2 module genes implied that CCNB2 may participate in the processes in the formation of CIS and tissue damage caused by CIS, such as “cell cycle,” “protein kinase activity,” and “glycosphingolipid biosynthesis.” Afterward, via single-cell RNA-seq analysis, CCNB2 was found up-regulated on GABAergic neurons in brain organoids as well as T cells expressing proliferative molecules in LUAD. Concurrently, the expression of CCNB2 distributed similarly to TOP2A as a module marker of cell proliferation in cell cluster. These findings can help in the field of the pathogenesis of LC-related CIS and neuron repair after CIS damage.

INTRODUCTION
The etiology of cerebral ischemic stroke (CIS) is complex and diverse. The traditional high-risk factors include hypertension, hyperlipidemia, and diabetes, among others . Recently, there have been increasing reports of CIS caused by malignant tumors (Salazar-Camelo et al., 2021;Wu et al., 2021b). Tumors disrupt organ function and affect the stability of the body's internal environment. Furthermore, the radiotherapy and chemotherapy drugs and surgery used in tumor treatment also significantly impact the body (Sooklal and Siddiki, 2021). Some retrospective studies have reported a risk of CIS in patients with malignant tumors, such as prostate cancer, head and neck tumors, colorectal cancer, and cervical cancer, and most of these patients have no traditional CIS risk factors (Aspberg et al., 2020;Ciolli et al., 2021), suggesting that the occurrence of CIS in some patients with malignant tumors is unrelated to the pathogenesis of traditional cerebrovascular diseases but is directly or indirectly related to tumors, which is called tumor-related CIS. Lung cancer (LC) is one of the tumors with the highest incidence worldwide, and its death rate is increasing. Even in some developing countries, LC is the malignant tumor most frequently occurring. Studies have confirmed the association between LC patients and increased risk of subsequent stroke (Chen et al., 2011;Anai et al., 2022). However, the detailed and incisive mechanism is still unclear. Therefore, clarifying the common pathogenesis of LC and CIS will be of great significance to the studies on the correlation between CIS and LC by providing them with a new theoretical perspective and actively promoting interdisciplinary research.
Recent studies have found that an uncontrolled cell cycle is a significant cause of unordered cell proliferation and carcinogenesis. Meanwhile, numerous experimental animal models confirmed the abnormal expression of cyclin in neurons after ischemia. Therefore, CCNB2 (CCNB2) was selected as our research focus. CCNB2 belongs to type B cell cycle family protein, which is located on chromosome 15q22, and it binds to cyclin-dependent kinases (CDKs) to regulate their activities (Wu et al., 2021a;Xia T. et al., 2021). Different cyclins have different spatiotemporal specificities at specific stages of the cell cycle (Gupta R. et al., 2021). CCNB2 participates in the G2/M phase transformation in the eukaryotic cell cycle by activating CDKs (Raina et al., 2021). The research of Meyer Da showed that the up-regulation of protein kinase CDK5 is the main reason for neuronal damage during CIS. Restraining abnormal CDK5 during CIS can protect dopamine neurotransmission, and obstruct excitatory toxicity (Meyer et al., 2014). Meanwhile, extensive studies have shown that overexpression of CCNB2 is associated with the development and deterioration of colorectal, lung, and breast cancers (Ma et al., 2021;Moradpoor et al., 2021;Wang et al., 2021;Toolabi et al., 2022). Consequently, CCNB2 can be used as a potential introduction to study the common pathways of the two diseases.
To describe the diversity of cell types in body or disease development and to reveal the functions and mechanisms of various types of cells, the gene expression differences between individual cells need to be explored. Single-cell transcriptome sequencing (scRNA-seq) technology analyzes transcriptome sequencing data at the single-cell level through unbiased, highthroughput, and high-resolution methodology (Ranzoni et al., 2021). Presently, the scRNA-seq method is extensively used in the fields of tumor, immunity, neural development, and animal genetics, among others (Fan et al., 2020). Therefore, scRNA-seq is an effective tool for mining molecular mechanisms.
In summary, the objective of this study is to evaluate the expression and clinicopathological significance of CCNB2 in CIS and LC using high-throughput datasets from public databases and to know preliminarily about the biological mechanism of CCNB2 in CIS through molecular classification and scRNAseq analyzes.

Differential Expression of Cyclin B2 Between Cerebral Ischemic Stroke and Lung Cancer Samples Compared With the Control Group
Two microarray matrixes containing 12 samples from First Affiliated Hospital of Guangxi Medical University [three lung adenocarcinoma (LUAD) vs. three non-LUAD, three lung squamous cell carcinoma (LUSC) vs. three non-LUSC] were used for detecting mRNA expression of CCNB2 (purchased from Fanpu Biotech, Guilin, China) and were included in our study. Afterward, the keywords "CIS" and "lung cancer" were used for data retrieval in the public databases of ArrayExpress 1 , Sequence Read Archive ( 2 SRA), and Gene Expression Omnibus (GEO 3 ). The inclusion criteria were as follows: (1) samples containing mRNA expression of CCNB2 in CIS, LC, and control tissues; (2) sample is a homo sapiens; (3) CIS and LC patients did not receive any treatment; and (4) sample is derived from serum, bodily fluids, or tissue of the patient. To illustrate the expression level of CCNB2 in the experimental and control groups as well as the discrimination ability of the samples in the experimental group, the included datasets were evaluated using integration analysis. The standard mean deviation (SMD) and its 95% confidence interval were calculated by extracting the mean ± SD value in each study. In addition, I 2 > 50% and p < 0.05 indicate that the study cohort has significant heterogeneity and thus the randomeffect model should be used; otherwise, the fixed model should be used. A test of publication bias was conducted to measure the heterogeneity of the cohort. p > 0.05 indicates no significant publication bias in the study cohort. Subsequently, a diagnostic accuracy analysis was conducted to judge the ability of CCNB2 to distinguish the experimental group from the control group, and the best cut-off value was determined through sensitivity (SEN) and specificity (SPEC). Furthermore, True Positive (TP), False Positive (FP), False Negative (FN), and True Negative (TN) were obtained to summarize the results of diagnostic analysis and calculate the SEN, SPEC, negative likelihood ratio (NLR), positive likelihood ratio (PLR), and diagnostic odds ratio (DOR).
Finally, the summary receiver operating characteristic (SROC) curve was drawn to obtain the area under the curve (AUC) of SROC according to previous reports He et al., 2021;Huang W. Y. et al., 2021;Liang et al., 2021a,b).

Expression of the Cyclin B2 Protein in Lung Tissue Compared With Lung Adenocarcinoma and Lung Squamous Cell Carcinoma
From the Human Protein Atlas (HPA) database (Colwill and Gräslund, 2011) 4 , the data of the CCNB2 protein expression in lung and LUAD as well as LUSC tissues were obtained. The detailed IHC procedure was conducted as described in the original reference (Colwill and Gräslund, 2011). The expression level was evaluated according to the intensity and density of immunohistochemical staining. For the scores of staining intensities, "Not detected" was 0 points, "Low intensity" was 1 point, "Medium intensity" was 2 points, and "High intensity" was 3 points. For the scores of staining densities, "Negative density" was 0 points, "Weak density" was 1 point, "Moderate density" was 2 points, and "Strong density" was 3 points. To obtain the aggregate score, the intensity and density scores were multiplied to divide the protein expression into 0-9 points.

Protein-Protein Interaction Network and Enrichment Analysis of Intersected Genes
To identify the potential interaction relationships between genes from the protein level, a PPI network was formed based on 123 intersection genes using the Search Tool for the Retrieval of Interacting Genes (STRING). In addition, according to the KEGG result, another PPI network was also constructed by STRING. To further explore the hub genes with the most potential among the selected genes, Cytoscape v3.7.2 was used to screen the hub genes. To further excavate the potential functions and related molecular mechanisms of the intersection of differentially expressed genes (DEGs) in LUAD, lung squamous cell carcinoma (LUSC), small cell lung cancer (SCLC), and CIS, we entered the intersection genes into R (Version 4.0.4) to conduct gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. GO terms and the KEGG signal pathway with p < 0.05 and FDR < 0.1 were identified. Similarly, we used R (Version 4.0.4) to search for the significantly enriched pathways from the Reactome database, and the pathways were selected when p < 0.05.

Regulatory Network of Cyclin B2 in Cerebral Ischemic Stroke
Pearson correlation analysis was conducted to analyze CNNB2 and other gene expressions in the CIS patient cohort. The top 2,000 significant genes were defined as the related genes of CCNB2. Gene co-expression network analysis can effectively identify functional co-expression gene modules related to complex human diseases. Multiscale embedded gene co-expression network analysis (MEGENA) adopts the network embedding paradigm in the topology field. Moreover, gene expression data contain information about the relationships between genes. Consistency clustering is an unsupervised clustering analysis that can identify interrelated genes and verify the rationality of clustering through a resampling-based method. It has a potentially high return in the prognostic analysis and treatment of cancer. Afterward, the R package ConsensusClusterPlus was adopted for the classification of clustering to obtain the K value with the best stability.

Mutation Analysis of Cyclin B2 in Lung Cancer and Clinical Parameter Analysis
A cell may encounter more than 700,000 DNA damages every day. If these damages cannot be repaired or are incorrectly repaired, they may lead to DNA nucleotide substitution, including single nucleotide variation (SNV), small or large insertions and deletions (indels), copy number variation (CNV), and chromosome rearrangement (GCR). Moreover, DNA replication, transcription, and recombination will destroy DNA stability and further increase the mutation load of the genome. Therefore, somatic mutations can reflect individual environmental exposure and the DNA repair process. R package TCGAbiolinks (Version 2.22.1) was used to download the mutation information of LC patients, and Mutect2 was used to analyze the mutation of CCNB2 in LC. In addition, datasets containing the pathological grade, age, and survival rate of patients were searched to analyze the clinical significance of the CCNB2 mRNA expression.

Single-Cell Transcriptome Sequencing (scRNA-seq) Analysis of Cyclin B2 in Brain Organoids and Lung Adenocarcinoma
The GSE184409 and GSE189357 dataset of GEO contains six brain organoids and nine LUAD tissues for dissecting the cells' transcriptional traits. The subpopulations of cells were obtained through linear dimensionality reduction principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) clustering. The marker genes and disease difference molecules in the cell subpopulations were analyzed and screened using the Wilcox rank-sum test. The marker genes of each population obtained in this dataset were compared with the marker molecules collected in the articles. According to the relative expression of the single-cell genes after integrating and correcting the batch effect, the differentially expressed gene set was selected as the variable to conduct the trajectory construction function and construct the single-cell development trajectory diagram in Monocle 2 software package, which was displayed with 2D visual results.

Statistical Analysis
In this study, alongside the statistical analysis explained above, when the CCNB2 expression presented a normal distribution, R 4.0.3 was used for Student's t-test. If the CCNB2 expression was abnormally distributed, a non-parametric test was used to compare the differential expression of CCNB3 in the experimental and control groups. GraphPad Prism 8.0 was employed to draw a box chart showing the distribution of the CCNB2 expression in the experimental and control groups. The AUC was used to evaluate the ability of the CCNB2 expression to distinguish between the experimental and control groups. STATA 14.0 was employed to calculate the SMD and 95% confidence interval, systematically reflecting the expression of CCNB2. sROC was calculated to illustrate the ability of CCNB2 to distinguish the experimental group from the control group. The Kaplan-Meier plotter was used to analyze the effect of the CCNB2 mRNA expression on the survival rate of patients. Using the median value, patients were divided into two groups according to the high and low expressions of the CCNB2 mRNA, a Kaplan-Meier curve was drawn, and the hazard ratio (HR) was calculated to judge the predictive ability of the CCNB2 mRNA expression on the patients' survival rate.

Evidence From Multiple High-Throughput Data Suggests That Cyclin B2 Is Significantly Up-Expressed in the Cerebral Ischemic Stroke and Lung Cancer Samples
Overall, 103 high-throughput datasets related to the study, including the CCNB2 expression, were bring into (Figures 1, 2 and Supplementary Table 1) while it failed to collect tissue samples of CIS patient and bodily fluids samples of LC patient. Each dataset was strictly under the corresponding CIS or LC diagnosis guidelines, and the experimental design of each included study was approved by the relevant review committee. In the included datasets, P < 0.05 indicated great heterogeneity, and thus the random-effect model was used. Integration analysis indicated that CCNB2 was up-regulated in the CIS and LC samples (Figure 3). Moreover, the funnel plot was distributed symmetrically, which indicated no significant publication bias in our study (Supplementary Figure 1A). Based on calculation of SEN and SPEC (Supplementary Figures 2, 3), the results of the SROC discriminative accuracy and the likelihood ratio reflected the fine identification ability of CCNB2 in CIS and LC (Figure 4 and Supplementary Figure 1B).

Verification of Upregulation of Cyclin B2 Protein Expression Level in Lung Adenocarcinoma and Lung Squamous Cell Carcinoma
Immunohistochemistry (IHC) was conducted to detect the CCNB2 expression in 12 LUAD and LUSC tissues as well as three non-tumor lung tissues. The representative pictures of IHC staining are shown in Figure 5A. CCNB2 was found to be negatively expressed in the lung tissues. Significantly, CCNB2 was expressed in partial LUAD and LUSC cells ( Figure 5B). After analyzing 15 samples, we discovered that the expression intensity of CCNB2 in the LUAD and LUSC tissues was significantly higher than that in the non-tumor lung tissues ( Figure 5C, p < 0.05). Based on the cut-off value of the immunohistochemical score, we could divide the samples into two groups, with AUC = 0.92 ( Figure 5D).

Protein-Protein Interaction Network and Enrichment Analysis Suggested Key Status of Cyclin B2 in Cerebral Ischemic Stroke and Lung Cancer
Concerning 123 intersected DEGs screened from LUAD, LUSC, SCLC, and CIS, a PPI network was built to express the connectedness and generality among the DEGs. Ten molecules with strong degree bounds were identified as central genes, including CCNB2, cell division cycle associated 8 (CDCA8), PDZ-binding kinase (PBK), kinesin family member 4A (KIF4A), ribonucleotide reductase regulatory subunit M2 (RRM2), ubiquitin-conjugating enzyme E2 C (UBE2C), cell division cycle 45 (CDC45), ZW10 interacting kinetochore protein (ZWINT), cyclin-dependent kinase inhibitor 3 (CDKN3), and cell division cycle associated 3 (CDCA3) (Figure 6A). To further explore the signal pathways of these genes, a functional and pathway enrichment analysis of the 123 genes was conducted. According to the result, these genes were significantly enriched in the cell cycle, regulation of the cell cycle process, and the cell cycle checkpoints, among others (The top three terms, Supplementary  Figure 4), which indicates CCNB2 exerts a significant role in the gene set.

Biological Functions of the Cyclin B2 Module Genes in Cerebral Ischemic Stroke
The CCNB2-related genes obtained by correlation analysis were expressed in a co-expression network (Figures 6B,C), and a module comprising 76 genes was obtained ( Figure 6D). The hub gene included thymidylate synthetase (TYMS), DNA topoisomerase II alpha (TOP2A), Opa-interacting protein 5 (OIP5), and cell division cycle 45-like (CDC45L). Functional enrichment analysis indicates that these genes are mainly enriched in 5'-3' exonuclease activity, protein serine/threonine/tyrosine kinase activity, catalytic activity acting on DNA, and cyclin-dependent protein serine/threonine kinase regulator activity, and histone kinase activity ( Figure 6E). Pathway analyzes of module genes were enriched in the cell cycle, oocyte meiosis, and p53 signaling pathway ( Figure 6E).
TOP2A was found to have the highest mutation frequency, including missense mutation, frameshift insertion mutation, and multihit, in LUAD and LUSC. More importantly, the expression of CCNB2 was related to age, pathological stage, overall survival (OS), and progression-free survival (PFS) of Euro-American LUAD patients (Figures 7B,C).

Genotyping Based on Cyclin B2 Co-Expressed Genes Has a Significant Pathway Enrichment Difference
Based on the consensus clustering method, the CIS samples were divided into two subgroups (Figures 8A,B). The Student's t-test indicated a distinct discrepancy in the distribution of CCNB2 ( Figure 8D) and pathway enrichment ( Figure 8C) in two clusters. It was demonstrated that the significant specialties of clusters with higher CCNB2 expression are higher activity folate biosynthesis, basal transcription factors, and Parkinson's disease (PD) pathway ( Figure 8E).

Single-Cell Transcriptome Sequencing (scRNA-seq) Analysis Suggests That Cyclin B2 Is Up-Regulated on GABAergic Neurons in Brain Organoids and T Cells in Lung Adenocarcinoma
The harmony function was adopted to rectify the sample heterogeneity of the expression matrix. The cells were identified as different clusters using an unsupervised clustering method, and the up-regulated molecules of the cell subpopulation were gained through differential screening. The up-regulated genes compared markers collected in the CellMarker database 5 and marker molecules collected in the articles to select markers for the study. Overall, the cells subpopulation was marked as cells (Figures 9A,D) in brain organoids and LUAD (Brain organoids: GABAergic neurons, glutamatergic neuron, neural progenitor, upper cortical layer, deep cortical layer, neuronal, astrocyte, fibroblast, and oligodendrocyte; LUAD: T cells, B cells, tumor cells, dendritic cells, macrophages, cancerassociated fibroblasts (CAFs), endothelial cells, mast cells, and AT1 cells) indicating the involvement of different biological processes (Figure 10). CCNB2 and TOP2A were significantly over-expressed in GABAergic neurons (Figures 9B,C) which suggested that CCNB2 was related to repair after CIS damage. Simultaneously, CCNB2 was found up-regulated on T cells expressing proliferative molecules including TOP2A, marker of proliferation Ki-67 (MKI67), and proliferating cell nuclear antigen (PCNA) in LUAD (Figures 9E,F and Supplementary  Figure 5). DEGs among the cell groups were taken as the 5 http://biocc.hrbmu.edu.cn/CellMarker/index.jsp reference gene set, then a pseudo-temporal algorithm was employed to predict the dynamic development trajectory of the cells (Figures 11, 12).

DISCUSSION
This study indicated the up-regulation of the CCNB2 mRNA expression in CIS and LC through multilevel evidence. Simultaneously, by analyzing the CCNB2 co-expression genes and scRNA-seq data, the potential biological pathways CCNB2 might involve in were explored, which provides a new perspective for researching the mechanisms of the shared biological pathways between CIS and LC. Our study has the following advantages. First, the study included 103 high-throughput datasets, including 6946 experimental and 5592 control samples, which strongly confirmed the up-regulation of CCNB2 mRNA in CIS and LC through integrated analysis. Additionally, it is of great significance to study LC using tissue samples, because the oncogenesis of LC originates from cytology carcinogenesis. Second, it is the first systematic study of CCNB2 considering the shared biological function between CIS and LC, thus enriching the theoretical content of cyclin family genes. Last, there is still a lack of genome-based and high-throughput data analysis methods to mine the biological mechanism of CCNB2. More importantly, the expression pattern of CCNB2 at the singlecell level in brain organoids and LUAD was explored. Overall, this study comprehensively analyzed the potential biological mechanism of CCNB2 in CIS and LC, which indicated that CCNB2 promotes LC and CIS by regulating the cell cycle, protein kinase activity, and interacting with TOP2A.    Initially, through PPI analysis of crossed high-expression genes in CIS and LC, we found that CCNB2 may be an important gene in the common pathway. Furthermore, using MEGENA analysis, a gene module of CIS including 76 genes was obtained. The function and pathway enrichment analysis of the CCNB2 module genes showed that CCNB2 may participate in the formation of CIS and tissue damage caused by CIS, such as "cell cycle, " "protein kinase activity, " and "glycosphingolipid biosynthesis" through the module gene network. Importantly, previous studies suggested that the change in protein kinase activity was related to the occurrence of CIS (Riebeling et al., 2021). Receptor interacting serine/threonine kinase 3 (RIPK3) was an important member of the serine/threonine protein kinase family and a key regulator of programmed necrosis . It was found that in the advanced stage of atherosclerosis, about half of the death of macrophages in the plaque was caused by programmed necrosis regulated by RIPK . Meanwhile, the expression of RIPK3 in carotid atherosclerotic plaque with necrotic part was significantly higher than that in other tissues (Tian et al., 2016). Mitogen-activated protein kinase activated protein kinase 2 (MK2) was mainly activated by mitogen-activated protein kinase 14 phosphorylation (MPK14) . It was also found in vivo experiments that hypercholesterolemic low density lipoprotein receptor (IDLR) (-/-)/MK2(-/-) mice could resist atherosclerosis with reduced accumulation of lipids and macrophages in the aorta (Jagavelu et al., 2007). Meanwhile, the study of Ebrahimian T showed that the activation of MK2 mediates the production of inflammatory mediators and increased blood pressure regulated by angiotensin (Ebrahimian et al., 2011). Simultaneously, MAPKs family played an important role in the process of neuronal apoptosis by participating in the regulation of apoptosis relatedgenes . Whether it was transient or permanent ischemic injury, MAPKs expression level altered significantly, and inhibiting its cascade reaction could improve the cerebral ischemic injury in animal models . Therefore, the regulation of mitogen-activated protein kinase pathway may be an effective way to inhibit cerebral ischemic neuronal injury and protect neural function. The cell cycle is divided into the G1, S, G2, and M phases, and the development of cells depends on the normal operation of the cell cycle, which is an internal component of cell differentiation and proliferation (Lezaja and Altmeyer, 2021). The cell cycle is mainly regulated by cyclins and CDKs. Cyclins are abnormally expressed in CIS, while CDKs play an important role in CIS injury (Saraireh et al., 2012). Inhibiting the activity of the cell cycle has a significant neuroprotective effect. Trimsit et al. found that the expression of CyclinD1 and CDK4 in the penumbra, after cerebral ischemia, is a key marker to evaluate whether ischemic neurons enter the cell cycle again (Wood and Endicott, 2018). Additionally, Katchanov et al. found that in the rat model of cerebral ischemia-reperfusion, the deletion of endogenous CDKI p16INK4a can lead to the delayed death of striatal neurons (Katchanov et al., 2001), indicating that there is a relationship between the cell cycle and neuronal injury. Concurrently, many studies have shown that the composition and metabolism of glycosphingolipids on the surface of the tumor cell membrane will undergo changes, such as the emergence of new glycosphingolipid antigens and the transformation of glycosphingolipid structure, among others . The alternation correlates with the occurrence and development of tumors and can be distinguished as a tumor marker (Moghimi et al., 2021). Furthermore, inhibiting the high expression of related glycosphingolipids will have a clear palliative effect on the development of tumors and will help to block the metastatic pathway of tumors. For example, Salk found that the sialylation level of the glycosphingolipid sugar chain of cells increased with high malignancy or high metastatic potential, indicating that sialic acid glycosphingolipid may promote the metastasis and deterioration of tumor cells . These studies have suggested a potential association between cell cycle activity and the related metabolism between LC and CIS.
Subsequently, TOP2A as a hub gene of the CCNB2 coexpression module in CIS was found to have the highest mutation frequency, including missense mutation, frameshift insertion mutation, and multihit, in LUAD and LUSC. DNA topoisomerase is the general name of enzymes that can catalyze the conversion of DNA topoisomerase (Crewe and Madabhushi, 2021). TOP2A widely exists in eukaryotic and prokaryotic cells and mainly regulates DNA topology by participating in DNA division, repair, recombination, replication, and transcription (Di Felice and Camilloni, 2021). Among them, TOP2A is mainly distributed in the nucleus and participates in the cell cycle pathway, which is closely related to cell proliferation and apoptosis . By analyzing the effect of the expression level of TOP2A protein on the proliferation and invasion of tumor cells, the researchers found that the expression of TOP2A protein was closely related to the proliferation and invasion of LC cells (Kanne et al., 2021). However, there are few studies on TOP2A in the nervous system. Recently, Thakurela found that small molecule inhibitors of TOP2 can destroy the formation of embryonic stem cells and lead to the down-regulation of genes of embryonic pluripotency (Thakurela et al., 2013). Studies by Lauren F Harkin have shown that TOP2A is involved in regulating the acquisition and maintenance of activated astrocyte stem cell characteristics, which may be mediated by transcription factor Paired Box 6 (PAX6) (Harkin et al., 2016). More importantly, CCNB2 and TOP2A were found specifically over-expressed in GABAergic neurons by scRNA-seq analysis. GABA is an inhibitory neurotransmitter mainly existing in the central nervous system and widely distributed in the brain and spinal cord (Sanchez-Vives et al., 2021). The GABAergic signaling system comprises glutamate decarboxylase (GAD), GABA, GABA transporter, and GABA receptor (GABAR) . In the cerebral cortex, the excitatory signals provided by glutamatergic neurons and the inhibitory signals provided by GABAergic neurons always complement each other to achieve the balance of excitation and inhibition (Bell et al., 2021). Jiang conducted an in-depth study on how the activity of GABAergic neurons affected the rehabilitation of neurological function after stroke (Jiang et al., 2017). The results showed that inhibiting the activity of GABAergic neurons in the striatum could promote the improvement of neurological function and reduce the volume of brain atrophy in post-stroke mice. In addition, these protective effects may be related to the inhibition of GABAergic neuron activity and the promotion of the nutritional factor fibroblast growth factor (FGF) up-regulation in vascular endothelial cells. The above-mentioned studies will help us to further understand the pathological process of CNS injury. Additionally, CCNB2 was found up-regulated on T cells expressing proliferative molecules. Pathway enrichment analysis of upregulated genes on the cluster indicates that CCNB2 promotes LUAD via interleukin 18 (IL-18) signaling pathway and vascular endothelial growth factor A (VEGFA)-VEGFR2 signaling pathway. As a proinflammatory cytokine, IL-18 plays an important role in the occurrence and development of infectious diseases and tumors . Gu RH found that IL-18 exerts a fine prognosis effect on LC patients (Gu et al., 2020). VEGFA is a protein and specific mitotic agent for endothelial cells (Ling et al., 2021). Combined with receptors, VEGFA can promote the division and proliferation of vascular endothelial cells and increase microvascular permeability (Kajal et al., 2021;Lee et al., 2021). Many studies have confirmed that VEGFA is up-regulated in nonsmall cell lung cancer (NSCLC) tissues and is also significantly positively correlated with lymph node metastasis of NSCLC (Wei et al., 2021;. At the same time, it is also an effective prognostic factor for NSCLC Shi et al., 2021). In summary, the above results indicate that CCNB2 can be a potential target for treatment strategies for CNS injury repair as well as LC.
Although we attempted to analyze the relationship between LC and CIS using CCNB2 on multiple levels, the following limitations still existed in this study. First, the clinical specimens of patients were not collected, which hindered the collection of clinical parameters, such as tumor, node, and metastasis (TNM) staging. Second, incomplete information prevented us from performing adequate follow-up of patients. There were also insufficient data in the major gene databases and online analysis websites to analyze the influence of CCNB2 expression on the prognosis of CIS and LC. Additionally, the lack of comparison of CCNB2 expression data between LC-related CIS and LC combined with CIS leads to the deficiency in the persuasion of this study. Overall, our results still need to be further confirmed by experiments in vitro and vivo.
In summary, multiple inferences were provided for the upregulation of CCNB2 in LC and CIS. The effect of CCNB2 on the common pathway between LC and CIS may indicate a foreground of CCNB2 as a biomarker and therapeutic target for LC-related CIS.

ETHICS STATEMENT
This study was conducted in accordance with the declaration of Helsinki and approved by the Ethics Committee of The First Affiliated Hospital of Guangxi Medical University (protocol number #2015-KY-Guoji-259), China.

AUTHOR CONTRIBUTIONS
M-JL, Z-GH, and GC conceived and designed the present study and revised the manuscript. S-BY analyzed the data and drafted the manuscript. G-SL drafted the manuscript. TW, D-SH, Z-JL, LC, and ZY guided statistical analysis. YY analyzed the data. G-YC, JW, B-TY, R-XX, and L-YL filtered and extracted data. All authors contributed to the article and approved the submitted version.


Next-generation sequencing in advanced Chinese melanoma reveals therapeutic targets and prognostic biomarkers for immunotherapy

Limited studies have interrogated the genomic landscape of Chinese melanoma in which acral and mucosal melanoma are the mainstay. In this study, we carried out a retrospective analysis on 81 Chinese melanoma patients (15 acral, 25 mucosal and 41 cutaneous melanoma). With the identification of 1114 mutations spanning 248 genes, we summarized that the mutation spectrum varied significantly by subtypes. Acral melanoma and mucosal melanoma had significantly more CNVs. MYC amplification was one of the most commonly detected CNVs, other frequent CNVs in mucosal melanoma included NBN and KDR, which were associated with the poor survival of melanoma patients. A generally low TMB, with a median of only 5.1 mut/Mb, was observed in three groups including cutaneous melanoma. Additionally, over 50% variants in DNA damage repair pathway were detected in all three subtypes, most of which were HRD related genes. Patients with alterations of HRD related genes had a longer survival time after immunotherapy. This study revealed a molecular profiling of Chinese patients with advanced melanoma, and proposed the high variant rate in DDR pathway as a biomarker of immunotherapy, which might provide therapeutic targets and guidance in making clinical decision for different Chinese melanoma.

The Cancer Genome Atlas (TCGA) research defined molecular subtypes of cutaneous melanoma on the basis of the presence of specific "driver" gene (BRAF, RAS, and NF1) mutations 1 . In ultraviolet-shielded melanoma, mutations of BRAF, NRAS or NF1 are less frequent compared to cutaneous melanoma, but the existence of other cancer driver gene mutations are detected. Accurate profiling of the spectra of mutational changes in melanoma facilitates personalized management of the disease. Due to the subtype bias, the overall mutation frequencies of BRAF and KIT in Asian patients is approximately 16-25% and 6-10% 2-4 , respectively. However, the overall frequencies of BRAF and KIT mutations in Caucasians are about 10-60% and 0-28%, respectively 5,6 . More than 50% of Asian patients fail to gain benefits from BRAF and c-KIT targeted therapy because of the low frequency rates of BRAF and KIT mutations.
Several inhibitors targeting at poly (ADP-ribose) polymerase (PARP) have already been approved by FDA or undergoing clinical trials in various diseases and treatment settings. It has been reported that tumors displaying DNA repair dysfunction might exhibit a BRCA-like behavior, according to the concept of "BRCAness" 7,8 . It is noteworthy that the efficacy of PARP inhibitors has also been reported in non-BRCA related tumors. Tumors with "BRCAness" might therefore benefit from PARPi treatment [9][10][11][12] . DDR (DNA damage repair) gene alterations ubiquitously exist in many cancer types, and significant enrichment of somatic mutations in DDR genes is  13 . However, the frequency of DDR gene mutations has been reported to be rare in melanoma, of which most were genes related to UV-induced DNA Damage. Elucidating gene alteration signatures in all subtypes is crucial for strategic decision in melanoma management. However, the majority of melanoma genomic sequencing data regarding cutaneous melanoma and limited small-cohort studies focused on acral and mucosal melanoma were mainly reported in western countries. In addition, the whole exome sequencing (WES) and whole genome sequencing (WGS) reported in the previous studies are hardly feasible in wide clinical application in China. Therefore, in this study, we conducted a genetic comparison among three subtypes of melanoma with a selected genes panel of next-generation sequencing (NGS) and evaluated its feasibility in practical application.

Results
Clinicopathological features. In cutaneous melanoma, 36.6% of patients had congenital nevus and 29.3% of patients had exposure history to sun, which were both significantly higher than patients with acral (p = 0.001) and mucosal melanoma(p = 0.001). While in acral melanoma, 66.7% and 40% of patients had ulceration and trauma respectively, which were higher than that in mucosal and CSD/NCSD (chronic sun-induced damage / non-chronic sun-induced damage) group (p = 0.020 and p = 0.010). Data were shown in Table 1.
Alterations of driver genes of cutaneous, mucosal and acral melanoma. A total of 1114 mutations spanning 248 genes were found in 81 Chinese melanoma patients, in which 95 genes were detected with frequency over 5% (Fig. 1a). Mutations of 28 genes were observed in at least 5% patients (Fig. 1b). In cutaneous (CSD/NCSD) melanoma, the most frequently mutated gene was BRAF (56.1%, 23/41), which was similar to that in western population. Among patients with BRAF mutations, BRAF-V600E/K mutation was found in 22 patients, p.P402L mutation was detected in one patient, and BRAF amplification was observed in 4 patients. The 7 cases of NRAS mutations were targeted to hotspots on codon 61 and codon 12, which were hot spot mutations of cutaneous melanoma. KIT mutations on p.Y936 and p.A89T were found in 2 cases. NF1 mutations were observed in 3 patients (7.3%), in which 1 was missense mutation and 2 were splicing variants. Driver gene mutations were not found in 17 (41.5%,17/41) patients. LRP1B (17.1%, 7/45) was the most frequently mutated gene after BRAF in cutaneous melanoma, which was equal to TP53 (17.1%,7/45).  www.nature.com/scientificreports/ In mucosal melanoma, the mutation frequency of BRAF was 12.0% (3/25), which was significantly lower than that in cutaneous melanoma (p = 0.001). In addition, a case of BRAF amplification and a case of KCTD15-BRAF fusion were observed. Two cases (8.0%) had mutations in NRAS of p.Q61R and p.G12D, which were consistent with hot spots in cutaneous melanoma. All 4 cases of KIT amplifications were found in mucosal melanoma, among which 1 case was co-occurred with p.W557G mutation. NF1 mutations were observed in 4 patients (16.0%), including missense mutations in 2 patients and splicing variants in 2 patients. TP53 (20.0%, 5/25) was the gene with the highest mutation frequency in mucosal melanoma.

Comparison of variants in pathways.
According to the classification of KEGG Pathway database, we compared the gene alterations in different pathways among three groups. Since BRAF was found to possess the highest mutation frequency, MAPK pathway was consequentially considered as the most activated pathway, thus BRAF was removed from the comparison of each pathway. All of the three groups had more than 75.0% PI3K-AKT pathway mutations. Acral group showed 93.3% (14/15) mutations in Ras pathway, which was higher than that in cutaneous melanoma (61.0%, 25/41, p = 0.023). Acral melanoma also had more mutations in WNT pathway and TGF-beta pathway (80.0%, 12/15 and 66.7%, 10/15), compared with that in cutaneous melanoma patients (39.0%, 16/41, p = 0.014 and 26.8%, 11/41, p = 0.01) (Fig. 3e).
Association of TMB and gene alterations with prognosis of immunotherapy. Among 81 patients enrolled, 25 of them received palliative immunotherapy including anti-PD-1 antibody (pembrolizumab or nivolumab), anti-CTLA4 antibody (ipilimumab) monotherapy or combinations. The clinicopathological information of patients receiving first-line immunotherapy was shown in Table S2.
The median of TMB in three groups was 5.1mut/Mb (Fig. 4a). We derived a cutoff of 15mut/Mb, which can effectively distinguish clinical response. Patients with TMB > 15 mut/Mb had a significantly longer PFS than patients with TMB ≤ 15 mut/Mb (p = 0.049) (Fig. 4b). High TMB was associated with improved PFS.
Then mucosal and acral melanoma which are dominant in Chinese population were combined into one group for further comparison, namely Acral/Mucosal group and CSD/NCSD group. We found that PFS in CSD/NCSD group was similar to that in acral/mucosal group (p = 0.711), as shown in Fig. 4c. Because of the high alteration rate of DDR pathway in our cohort, especially in HRD related genes (CHEK1, CHEK2, BRCA1, BRCA2, ATM, PALB2, BAP1, IDH1), patients were divided into two groups, and we found that patients with alterations of HRD related genes showed a longer PFS (9.7 vs. 25.3 months, p = 0.040) and a prolonged survival time after immunotherapy (p = 0.117) (Fig. 4d,e). Then we validated our findings by other researches. This finding was then validated by gene mutation data of patients in other melanoma groups (UCLA, Cell 2016 + MSKCC, NEJM 2014 + DFCI, Science 2015, n = 202) [15][16][17] and patients with NSCLC (MSK, Cancer Cell 2018, n = 75) 18 and glioblastoma (Columbia, Nat Med. 2019, n = 42) 19 who accepted immunotherapy, taken from a publicly available database (CBioPortal) (Fig. 4f-h). DDR pathway alterations may serve as a biomarker of immunotherapy.

Discussion
In the current study, we performed NGS sequencing with a multiple-gene panel to investigate the comprehensive molecular characterization of 81 Chinese melanoma patients and evaluated the clinical correlations of gene status, we also evaluated the correlations between the response of immunotherapy and gene alterations. This study aimed to exploit tumor targeted NGS to compare different melanoma types in Chinese population. The mutation landscape of Chinese melanomas differed from that of western population, which was distinguished by melanoma  Fisher's exact test was performed, and the statistical significance was defined as *p < 0.05, **p < 0.01, ***p < 0.001.  20 , while Triple-WT was higher (37.8%, 17/45). Mucosal and acral melanoma are the main subtypes of Chinese melanoma. Compared with cutaneous melanoma, BRAF mutations were less observed in acral melanoma, while mutant NRAS (26.7%) and Triple-WT (7/15, 46.7%) were more common. In mucosal melanoma, mutations in BRAF (12.0%) and NRAS (8.0%) were both lower than that in western population, and lower than Chinese cutaneous melanoma. NF1 mutation (16.0%) of mucosal melanoma was similar to that in western population, but higher than that of Chinese cutaneous melanoma. Consistent with previous studies, our study further confirmed that gene mutations were more common in cutaneous melanoma, while more CNVs were observed in mucosal and acral melanoma 1,21,22 . In our study, the most common CNVs in mucosal and acral melanoma was MYC amplification, consistent with previous study in Asian patients 20 . In vitro experiments have confirmed that high expression of c-MYC is positively correlated with the aggressiveness of cutaneous melanoma, and the inhibitor can effectively inhibit tumor growth 23 . Although the role of MYC in the development and progression of mucosal and acral melanoma remains to be further studied, MYC may act as a new therapeutic target for the treatment of mucosal and acral melanoma with the application of effective inhibitors. NBN followed MYC also had a high CNVs rate in mucosal and acral melanoma. NBN is an important gene which is related to DNA damage repair, and plays an important role in protecting chromosome integrity. Currently, the variation and role of NBN in melanoma has not been reported. The variation frequency of KDR (VEGFR-2) in mucosal melanoma was only second to that of MYC and NBN. KDR is a main functional receptor of VEGF which plays a role in angiogenesis growth of tumor cells 24 . The cBioPortal database showed that MYC, NBN and KDR predicted poor survival of melanoma patients. Guo Jun et al. reported that Axitinib combined with PD-1 antibody achieved ORR at 48.3% in mucosal melanoma 25 , suggesting that Axitinib and other multi-target small molecule inhibitors www.nature.com/scientificreports/ targeting at angiogenesis may be used in mucosal melanoma and as a combination to improve the efficacy of immunotherapy. Guo Jun et al. reported that genetic aberrations in the cyclin-dependent kinase (CDK)4 pathway occur in 82% of patients with acral melanoma, showed CDK4 gain (39.5%) and CCND1 gain (26.7%). CCND1 amplification was found in 40.0% (6/15) patients in our acral group, similar to their report. Previous studies have reported that, TERT promoter mutations are more common in acral and mucosal melanoma patients 21 , and SF3B1 mutation are more common in mucosal melanoma, mainly in Caucasians 21 . However, because of the limitation of our panel, these genes are not detected in our study. Previous studies have shown that patients with higher TMB were more likely to gain a better efficacy of immunotherapy 16,[26][27][28] . Pembrolizumab has been approved by the FDA for the treatment of solid tumors with TMB over 10 mut/Mb 29 . TMB plays an increasingly important role in immunotherapy. However, in our study, a generally low TMB was observed in Chinese melanoma patients, and the median of TMB was only about 5.1 mut/Mb, with no difference found in cutaneous, mucosal and acral melanoma. In KEYNOTE-151, the overall response rate of second-line Pembrolizumab for Chinese patients was 16.7%, and first line anti-PD-1 antibody for Chinese patients was 15.6%, both of which were lower than that of western population [30][31][32] . The generally low TMB in Chinese population and poor immunogenicity might be the reason why the efficacy of immunotherapy in Chinese population, including cutaneous melanoma, is worse than that in western population.
PARP inhibitors are considered to be effective drugs for the treatment of BRCA germline mutations. It has been reported that PARP inhibitors also play a role in tumors with no BRCA germline mutations but BRCA somatic mutations or other HRD mutations 9,10,13 . In our study, we observed a high mutation frequency at 14.8% in HRD pathway, except for BRCA1 (1/81), BRCA2 (4/81), other genes like ATM (4/81), PALB2 (2/81), CHEK2 (2/81), BAP1 (4/81) and IDH1 (1/81) were included. It is reported that patients with homologous recombination deficiency (HRD) pathway mutations usually had an increased burden of neoantigens 33 . PARP inhibitors may be applied to patients with these variants to further improve the efficacy of immune checkpoint inhibitors.
To conclude, we compared the mutation profiling of three main subtypes of Chinese melanoma. We observed a generally low TMB of Chinese melanoma patients, but a high variants rate in DDR pathway, especially in HRD related genes, which may contribute to the exploration of new drugable targets. Additionally, this study highlighted the importance of implementing next-generation sequencing testing.

Methods
Ethics statement. This study was conducted in accordance with the 1964 Helsinki Declaration. All human studies were approved by the Ethics Committee of Sun Yat-sen University Cancer Center (GZR2017-207). Written informed consent was obtained from eligible patients.

Patient selection and sample collection.
Formalin-fixed paraffin-embedded tissues were obtained from 84 Chinese patients diagnosed with melanoma between September 2017 and September 2021 in the authors' clinic in Sun Yat-sen University Cancer Center. Standard histopathology was performed to confirm the diagnosis of malignancy and histologic subtype. Among 84 analyzed cases, 3 samples were excluded due to insufficient DNA quantity. A total of 81 melanoma cases were enrolled, including 25 (30.9%) mucosal melanoma, 15 (18.5%) acral melanoma and 41 (50.6%) cutaneous melanoma (namely CSD and NCSD melanoma). The median age of patients of each subtype at diagnosis was 44, 53 and 55 years old respectively.
Tissue DNA extraction, NGS detection and sequencing data analysis. DNA was extracted with the QIAamp DNA FFPE tissue Kit (Qiagen) according to the manufacturer's instructions. DNA concentration was measured by Qubit dsDNA assay.
Genetic profiles of all tissue samples were assessed by performing capture-based targeted deep sequencing with the OncoScreen panel (Burning Rock Biotech Ltd.) which covered 2.02 MB of human genomic regions, including all exons and critical introns of 295 genes, and genes included in the panel were listed in Table S1. Details of sequencing data analysis were described as previously reported 34 . The calculation of tumor mutation burden (TMB) was based on the ratio of the total number of mutations to the size of panel.
Validation data collection and analysis. The data of targeted gene mutations and survival data were obtained from cBioPortal database, and analyzed on cBioPortal, samples were divided into 'mut' or 'wt' . Gene mutation data and survival data of MYC, NBN and CCND1 were from the study of Melanoma (MSKCC, Clin Cancer Res 2021, n = 696) 14 . DDR related gene (CHEK1, CHEK2, BRCA1, BRCA2, ATM, PALB2, BAP1, IDH1) mutation data and survival data of patients treated by immunotherapy were also obtained from cBioPortal, melanoma (UCLA, Cell 2016 + MSKCC, NEJM 2014 + DFCI, Science 2015, n = 202) [15][16][17] , NSCLC (MSK, Cancer Cell 2018, n = 75) 18 and glioblastoma (Columbia, Nat Med. 2019, n = 42) 19 . P < 0.05 was considered to indicate a statistically significant difference. Statistical analysis. Patient follow-up data were acquired from medical records. The χ 2 test and Fisher's exact tests were applied to analyze the association. The significance of the association of the mutations between the three groups was analyzed using Fisher's exact test. The Mann-Whitney test was used to compare tumor mutation burden. The Kaplan-Meier method was utilized to conduct survival analysis. All statistical analyses were accomplished by SPSS V.20.0 software. P value < 0.05 was considered to be statistically significant.


The reference genome and organelle genomes of wasabi (Eutrema japoniacum)

Wasabi (Eutrema japonicum) is one of the most famous vegetable crops in the family Brassicaceae. However, a limited genomic resource is available, which hinders genomic breeding and understanding of the genetic basis of vital traits. Here, we generated the genome assembly of wasabi using the hybrid genome assembly strategy, which combined the Nanopore long reads and Illumina reads. The genome assembly contains 687M bp and 39,534 high-quality annotated gene models. Besides, we annotated 68.85% of the genomic sequences as repetitive elements, including 43.72% of retrotransposons and 18.99% of DNA transposons. Using the customized pipeline, we also generated the complete organelle genomes of wasabi. This reference genome could provide essential genomic resources for evolution, breeding, and exploring the unique biological traits of wasabi.

Introduction
Wasabi (Eutrema japonicum) is a perennial herb that plays a vital role in Japanese cuisine and culture (Yamane et al., 2016;Ramirez et al., 2020). A paste made from ground rhizomes is used as a pungent condiment, similar to hot mustard or horseradish (Ramirez et al., 2020;Yang et al., 2020). Wasabi is a plant of the family Brassicaceae (Yamane et al., 2016). Brassicaceae is a medium-sized and economically important family of flowering plants commonly known as the mustards, the crucifers, or the cabbage family (Guo et al., 2017;Nikolov et al., 2019). The family contains cruciferous vegetables, including species such as Brassica oleracea (cabbage), Brassica rapa (turnip), Brassica napus (rapeseed), Raphanus sativus (radish), Armoracia rusticana (horseradish), but also a cut-flower Matthiola (stock) and the model organism Arabidopsis thaliana (Guo et al., 2017;Nikolov et al., 2019). Despite its high economic value, only a few genomes from the family have been deciphered, e.g., Brassica oleracea, Brassica rapa, and Arabidopsis thaliana (Wang et al., 2011;Liu et al., 2014).
Eutrema is a genus in the Brassicaceae that contains more than 30 species, most of which are distributed in eastern Asia, e.g., the model plant E. salsugineum, which is used for studies of abiotic stress. Most of these species occur in the Qinghai-Tibet Plateau (QTP). Several studies have been conducted to understand the genomes of Eutrema. Yang et al. (2013) presented the reference genome sequence (241 Mb) of E. salsugineum at 8 × coverage sequenced using the traditional Sanger sequencing-based approach. Guo et al. (2018) provide de novo whole-genome assemblies for a pair of recently diverged perennials with contrasting altitude preferences, the high-altitude E. heterophyllum from the eastern Qinghai-Tibet Plateau and its lowland congener E. yunnanense.
As for the high commercial crop, the chloroplast genome sequence of wasabi and its relatives were decoded to understand the evolution and phylogeny (Haga et al., 2019). Besides, the high repeat sequence content in its relative species also suggests the potential high repeat in the wasabi genome (Guo et al., 2018). This study performed hybrid genome de no assembly towards wasabi, combining Nanopore long reads with illumine short reads. After detailed genome annotation, we obtained a highquality dataset of the wasabi genome, providing us with a powerful tool to conduct related genomic studies and breeding programs.

Plant materials
The plant material was obtained from the experimental plantation in Liangshan (27.88°N, 102.26°E), Sichuan, China. We selected a well-planted Eutrema japoniacum sample for all experiments. Three different tissues (leaf, root, and stem) were collected for RNA sequencing. Fresh and healthy samples were harvested and frozen in liquid nitrogen immediately after collection and stored at −80°C in the laboratory.

Whole genome sequencing and RNA sequencing
Genomic DNA was extracted with a QIAGEN ® Genomic Kit from leaves. The purified DNA was then prepared using the SQK-LSK109 genome sequencing kit protocol [Oxford Nanopore Technologies (ONT), Oxford, United Kingdom]. Single-molecule real-time sequencing of long reads was performed on a PromethION platform (ONT, Oxford, United Kingdom). Additionally, the Illumina sequencing libraries were constructed following the manufacturer's instructions. The libraries with an insert size of 300 bp were created using the TruSeq Sample Preparation kit and sequenced on the Illumina HiSeq X-ten platform (Illumina, San Diego, CA).
The samples of different tissues (leaf, root, and stem) were prepared for library construction. The quality and quantity of the RNA were evaluated using a NanoDropTM D-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE), a Qubit ® 3.0 Fluorometer (Thermo Fisher Scientific, United Stated), and an Agilent Bioanalyzer 2100 (Agilent Technologies, CA, United States). The constructed libraries were sequenced using the similar method above.
Genome survey, genome assembly, and genome annotation We counted the K-mers using Jellyfish software and calculated the characteristics of the genome using Genomescope 2.0 software (Vurture et al., 2017). We performed genome assembly using MaSuRCA software (Zimin et al., 2013), which adopted a hybrid strategy. Genome annotation was conducted using the Maker-P pipeline (Campbell et al., 2014). Using the RNA-seq reads from different libraries, we performed de novo assembly with SOAPdenovo-Trans (Xie et al., 2014). To identify other kinds of repeat sequences in the wasabi genome assembly, we first built a non-redundant repeat sequence library by searching for repetitive sequences using the EDTA pipeline (Ou et al., 2019). The tRNA genes were predicted using the tRNAscan-SE package (version 1.3.1) with default parameters (Chan et al., 2021). The rRNA genes (8S, 18S, and 28S) were predicted using RNAmmer algorithms with default parameters (Lagesen et al., 2007). The miRNA and snRNA were identified using INFERNAL (version 1.

Assessment of assembly quality
We assessed the completeness of wasabi genome assembly using BUSCO software (version 4.0.5) with the parameters "-l embryophyta_odb10 -g genome" (Simão et al., 2015). We performed a sequence identity assessment by aligning the non-redundant Illumina short reads to the wasabi genome assembly using BWA with default parameters (Li and Durbin, 2009). Finally, RNA-seq data of different tissues were aligned to the wasabi genome using Hisat2 with default parameters (Kim et al., 2019).

Gene functional annotation
To annotate the genes, we identified the homologous genes with the National Centre for Biotechnology Information (NCBI) Non-Redundant (NR) Protein Database, the Universal Protein Resource Knowledgebase (UniProtKB), Swiss-Prot Protein Database using BLASTP with an E-value threshold of 1e-5 (Altschul et al., 1990). Functional domains were identified using InterProScan (version 5.2-45.0) against several publicly available databases (Blum et al., 2021). Metabolic pathway annotations were performed by sequence De novo assembly of organelle genome Firstly, we used the conserved gene sequences from organelles as baits files to extract Nanopore organelle reads including mitochondria and chloroplast. Secondly, the long reads were assembled using Flye software (Kolmogorov et al., 2019) and polished by Racon software (Vaser et al., 2017). The assemblies were used as the new baits file to re-extract nanopore organelle reads and Illumina reads. Finally, the Unicycler software (Wick et al., 2017) was used to conduct the organelle genomes. Annotation was performed with GE-Seq using default parameters to predict protein-coding genes, tRNA genes, and ribosomal RNA (rRNA) genes (Tillich et al., 2017). Manual annotation was performed for genes with low sequence identity to determine the positions of start and stop codons depending on the translated amino acid sequence using the chloroplast/bacterial genetic code. OrganellarGenomeDRAW (OGDRAW) was optimized to create detailed, high-quality maps of organellar genomes (Greiner et al., 2019).

Genome survey and genome assembly
To obtain the wasabi genome assembly, we generated a set of 122fold-coverage Illumina paired-end short reads (75.62 Gb) and 30fold-coverage Nanopore long reads (26.83 Gb) (Supplementary Table  S1). We firstly conducted the genome survey using the K-mer analysis. The wasabi genome size and the genome heterozygosity rate were estimated to be 640 Mb and 0.68%, respectively (Supplementary Figure S1). We assembled the reads into 687 Mb of contig sequences using the hybrid genome assembly strategy (Supplementary Table S2). We noticed that the size of genome assembly is slight larger that the estimated genome size, which may be caused by the assembly duplications or the underestimation of the genome survey. The genome assembly contains 7,875 contigs with an N50 length of 356,067 bp, and the largest contigs were 5,563,810 bp (Supplementary Table S2). We identified 90% of BUSCO gene models in the wasabi genome assembly (Supplementary Table S3), suggesting the completeness of the genome assembly. We furthermore remapped the sequencing reads to the genome to investigate the completeness. Nearly 96% of the DNA Illumina reads, and 99% of the Nanopore reads could be mapped to the genome. Moreover, an average of 80% of RNA-seq reads were aligned to the genome (Supplementary Table  S4). All of the statistics reflected the relative completeness of the wasabi genome.

Genome annotation and gene family analysis
We first generate the high-quality library using the EDTA pipeline to produce a high-quality annotation of the repetitive sequences. We identified 9,042 intact repetitive sequence elements in the genome, including 6,241 LTR-RTs and 3,161 DNA transposons (Supplementary Table S5). Specifically, there were 1,302 Copia type-, 2583 Gypsy type-and unknown type-LTR-RTs (Supplementary  Table S5). In the wasabi genome, repetitive sequences occupied 473.7 Mb (68.85% of the 687.9 Mbp). The most abundant repeats in this genome are retrotransposons (43.72% of genome assembly) ( Figure 1B). The non-LTR retrotransposons (LINE and SINE) accounted for 0.27%. 18.99% of the genome sequences belong to the DNA transposons (14.99% for TIR and 4.00% for Helitrons) (Supplementary Table S6).
Using the RNA-seq data, we obtained a total of 120,945 transcripts. After the genome annotation, A total of 39,534 genes were predicted, 99.9% (39,495) of which could be functionally annotated referring to the currently available databases, and 84. 2% (33,279) were expressed in at least one tissue (Supplementary Table S7). These results would provide valuable genetic resources for future functional genomics and molecular breeding research.
We analyzed the gene families using the genomes of wasabi and other 11 other representative plant species, including Fragaria vesca, Vitis vinifera, Solanum tuberosum, Actinidia deliciosa, Brassica oleracea, Eutrema salsugineum, Arabidopsis halleri, Arabidopsis thaliana, Boechera stricta, Capsella grandiflora and Capsella rubella. Three hundred seventy-seven thousand three hundred one genes were clustered into 27,255 orthologous groups (Supplementary Table S8). Seven thousand one hundred ninety-six groups were present in all 12 species, with 381 species-specific orthologous groups (Supplementary Tables S8, S9). We identified a total of 1,030 wasabi specific orthologous groups containing genes with a range of biological processes and 70% of them were annotated (Supplementary Table S10). The phylogenetic tree constructed using the single-copy genes indicated the phylogenetic position of the wasabi ( Figure 1C).

Organelle genomes
We conducted the organelle genomes using the customized pipeline, which integrated the Nanopore sequencing data and the Illumina data. The Nanopore long reads from the organelle genomes were fist conducted assembly to obtain the raw assembly. The Nanopore reads and Illumina reads with high-accuracy were integrated to obtain the final assembly. The customized pipeline could take advantage of the continuity of long reads and the highaccuracy of short reads. As a result, the length of the chloroplast genome was 153,847 bp and was circular, suggesting the completeness of the chloroplast assembly (Supplementary Table S11 and Supplementary Figure S2). The chloroplast genome exhibited a typical quadripartite structure, consisting of a pair of inverted  Figure S3). We conducted the phylogenetic of both the organelle genomes with relative species indicated the phylogenetic position of the wasabi.

Author contributions
HL, XY, and QQ conceived the study. YG, BY, QX, and XZ performed bioinformatics analysis. XZ, HL, and XY extracted DNA and RNA. BL, and XZ wrote the manuscript. All authors approved the final submission.

Funding
This work was financially supported by the Scientific Research Foundation of Yulin Normal University for high-level talents (G2019ZK13, G2020ZK06, G2019ZK35); National Natural Frontiers in Genetics frontiersin.org