BCL6 is regulated by the MAPK/ELK1 axis and promotes KRAS-driven lung cancer

Mutational activation of KRAS is a common oncogenic event in lung cancer, yet effective therapies are still lacking. Here, we identify B cell lymphoma 6 (BCL6) as a lynchpin in KRAS-driven lung cancer. BCL6 expression was increased upon KRAS activation in lung tumor tissue in mice and was positively correlated with the expression of KRAS-GTP, the active form of KRAS, in various human cancer cell lines. Moreover, BCL6 was highly expressed in human KRAS-mutant lung adenocarcinomas and was associated with poor patient survival. Mechanistically, the MAPK/ERK/ELK1 signaling axis downstream of mutant KRAS directly regulated BCL6 expression. BCL6 maintained the global expression of prereplication complex components; therefore, BCL6 inhibition induced stalling of the replication fork, leading to DNA damage and growth arrest in KRAS-mutant lung cancer cells. Importantly, BCL6-specific knockout in lungs significantly reduced the tumor burden and mortality in the LSL-KrasG12D/+ lung cancer mouse model. Likewise, pharmacological inhibition of BCL6 significantly impeded the growth of KRAS-mutant lung cancer cells both in vitro and in vivo. In summary, our findings reveal a crucial role of BCL6 in promoting KRAS-addicted lung cancer and suggest BCL6 as a therapeutic target for the treatment of this intractable disease.

Introduction
Lung cancer is the leading cause of cancer incidence and mortality worldwide (1). The overall 5-year survival rate of lung cancer remains low (2). Previous molecular studies have identified several oncogenic drivers and promoted the development of satisfactory treatments for lung cancer driven by alterations in EGFR, ALK, RET, or ROS1 (3). However, few effective therapies are available for gain-of-function mutations in KRAS, which occur in approximately 25% of all lung cancer cases. Knowledge of KRAS protein structure, dynamics, and signal transduction remains unmet, largely preventing the development of specific inhibitors that target this oncogene directly or indirectly. The selective KRAS inhibitor (KRASi) sotorasib, which specifically targets KRAS(G12C) by creating a stable covalent bond with a mutant cysteine residue, has been approved by the US Food and Drug Administration as a second-line treatment for locally advanced or metastatic nonsmall cell lung cancer carrying the KRAS(G12C) mutation (4). However, drug resistance inevitably occurs during treatment (5,6). Importantly, other frequently mutated forms of KRAS in lung cancer, such as KRAS(G12D) and KRAS(G12V), remain undruggable (7). Targeting downstream effectors of KRAS, such as MEK, is an alternative strategy to inhibit oncogenic KRAS signaling. However, the use of MEKi as a single agent or in combination with other clinical drugs has failed to demonstrate significant survival outcomes in KRAS-driven lung cancer, partly because of reactivation of the MAPK pathway component or reduced dominant intratumoral T cell clones (8,9). Therefore, identification and further development of promising therapeutic strategies are urgently needed related to either novel KRAS pathway components or KRAS-related proteins that are synthetic lethal with mutant KRAS. B cell lymphoma 6 (BCL6) was initially discovered as a protooncogene and is required for the formation of germinal center and development of B cells (10,11). During clonal expansion and immunoglobulin affinity maturation, BCL6 enables B cells to tolerate DNA damage by regulating DNA damage-sensing genes (ATR, TP53, ARF, and others) (10). BCL6 also favors the survival and proliferation of germinal center B cells by modulating proliferation checkpoint genes (CDKN1A, CDKN2A, PTEN, and others) (12) and numerous target oncogenes (MYC, BCL2, BMI1, and others) (13,14). Therefore, BCL6 has been characterized as a critical transcriptional repressor in germinal center. However, constitutive expression of BCL6 by chromosomal translocations or aberrant somatic mutations causes diffuse large B cell lymphoma (12). Current inhibitors either block the interaction between the BTB domain and correpressors or promote BCL6 degradation. These have shown antitumor effects against BCL6-addicted tumors (15,16). Recent studies have implicated genomic amplification of the BCL6 locus in certain solid tumors such as glioma, breast cancer, and ovarian cancer (17)(18)(19). Moreover, the combined inhibition of BCL6 and STAT3 synergistically defeats intratumoral heteroge-Mutational activation of KRAS is a common oncogenic event in lung cancer, yet effective therapies are still lacking. Here, we identify B cell lymphoma 6 (BCL6) as a lynchpin in KRAS-driven lung cancer. BCL6 expression was increased upon KRAS activation in lung tumor tissue in mice and was positively correlated with the expression of KRAS-GTP, the active form of KRAS, in various human cancer cell lines. Moreover, BCL6 was highly expressed in human KRAS-mutant lung adenocarcinomas and was associated with poor patient survival. Mechanistically, the MAPK/ERK/ELK1 signaling axis downstream of mutant KRAS directly regulated BCL6 expression. BCL6 maintained the global expression of prereplication complex components; therefore, BCL6 inhibition induced stalling of the replication fork, leading to DNA damage and growth arrest in KRAS-mutant lung cancer cells. Importantly, BCL6-specific knockout in lungs significantly reduced the tumor burden and mortality in the LSL-Kras G12D/+ lung cancer mouse model. Likewise, pharmacological inhibition of BCL6 significantly impeded the growth of KRAS-mutant lung cancer cells both in vitro and in vivo. In summary, our findings reveal a crucial role of BCL6 in promoting KRAS-addicted lung cancer and suggest BCL6 as a therapeutic target for the treatment of this intractable disease. upregulated BCL6 expression in the context of KRAS mutational activation. Increased BCL6 expression enhanced the expression of prereplication complex (preRC) components, thereby maintaining cellular replication origins and cell cycle progression. Importantly, BCL6 was functionally required for KRAS-mutant lung tumorigenesis and tumor growth, as the genetic or pharmacological inhibition of BCL6 significantly impeded the growth of KRAS-mutant lung cancer in preclinical mouse models. Our findings provide a direct link between BCL6 and KRAS oncogenic signaling and highlight that BCL6 is an important target for the treatment of KRAS-driven lung cancer. neity in a subset of non-small cell lung cancers (20). We recently revealed that BCL6 enables solid tumor cells to evade genotoxic stress (21), and pharmacological inhibition of BCL6 enhances the sensitivity of KRAS-mutant lung cancer cells to clinical BETi (22). However, the biological function of BCL6 in KRAS-mutant lung tumorigenesis remains unclear. The fact that BCL6 functions as a protooncogene and is required for malignant transformation prompts us to speculate that BCL6 may be linked to oncogenic responses in lung cancer.
In this study, we demonstrated that the MAPK/ETS transcription factor ELK1 (MAPK/ELK1) signaling axis directly . LSL-Kras G12D mice were infected with intranasal Ad-Cre for 16 weeks. (C and D) mRNA (C) and protein (D) expression levels of BCL6 in MEFs. MEFs were isolated from LSL-Kras G12D mice and infected with Ad-Cre for 72 hours (n = 3). (E) Correlation analysis. Correlation between BCL6 and KRAS-GTP protein expression levels in different cancer cell lines (n = 20). KRAS-GTP levels were determined using GST-RBD, the GST-fusion of the RAS binding domain of C-RAF, to pulldown active GTP-bound RAS from cellular lysates by glutathione beads. R, Pearson's correlation coefficient. (F and G) KRAS mutant variants upregulated BCL6 expression at the mRNA (F) and protein (G) levels. 293T cells were transfected with different doxycycline-inducible G12 KRAS mutant variants: pLVX-TetOne-KRAS(G12C), pLVX-TetOne-KRAS(G12D), pLVX-TetOne-KRAS(G12S), pLVX-TetOne-KRAS(G12V). Transfected cells were then treated with or without 2 μM doxycycline (Dox) for 96 hours. KRAS activity was determined using GST-RAF-RBD to pulldown GTP-KRAS from cell lysates. (H) Heatmap showing differentially expressed BCL6 target genes (P < 0.05) in HPNE and HPNE/KRAS cells (n = 3). z score was calculated based on counts of exon model per million mapped reads. (I) The mRNA expression levels of BCL6 target gene in HPNE and HPNE/KRAS cells (n = 3). Data in A, C, F, and I are expressed as mean ± SEM of 3 technical replicates, representative of 3 independent experiments. *P < 0.05, **P < 0.01, ***P < 0.001, determined via unpaired 2-tailed Student's t test. The immunoblots in B, D, and G were contemporaneous and run in parallel from the same biological replicate. J Clin Invest. 2022;132(22):e161308 https://doi.org/10.1172/JCI161308 The MAPK/ERK signaling axis contributes to BCL6 expression. Next, we explored how the oncogenic KRAS protein promotes BCL6 expression using genetic and chemical approaches. We genetically silenced KRAS and BCL6 using siRNAs in H460 and H441 cells. Our results demonstrated that KRAS deficiency abolished both the mRNA and protein expression of BCL6 ( Figure 2, A and B), whereas BCL6 silencing did not affect KRAS expression ( Figure 2C), indicating that BCL6 was downstream of KRAS. To further delineate which KRAS effector pathway contributes to BCL6 expression, we employed a small library consisting of 48 small-molecule inhibitors targeting crucial KRAS effectors, including MAPK, phosphatidyl inositol 3-kinase, and Ral pathway components (Supplemental Table 1). We observed decreased BCL6 protein expression in treated H460 cells and demonstrated that the MEKi selumetinib was the most effective hit. Intriguingly, 20 top-ranked hits from this screen were all MAPK/ERK pathway inhibitors, including RAFi, MEKi, and ERKi ( Figure 2D), indicating the positive regulation of the MAPK/ERK axis on BCL6 expression. Moreover, BCL6 expression was reduced at the transcriptional level ( Figure 2E).
BCL6 has been reported to directly bind to 2 specific sites within the TP53 promoter region and further suppresses its transcription (24). Next, we performed a TP53 reporter assay (Supplemental Figure 2A) and examined the effects of the MAPK/ERK pathway inhibitors on BCL6-mediated TP53 transcription. Luciferase reporter assays showed that most MAPK/ERK pathway inhibitors led to an increase in reporter activity ( Figure 2F), reinforcing an essential role of the MAPK pathway in regulating BCL6 biological activity. In line with this, siRNAs targeting RAF1, MAP2K1/2, and MAPK1/3 dramatically inhibited BCL6 protein and mRNA levels in H460 (Figure 2, G and H) and H441 cells (Supplemental Figure  2, B and C), respectively. In contrast, constitutively active variants of HRAS(G12V), BRAF(V600E), MEK(C121S), or MEK(S222D) significantly augmented BCL6 expression ( Figure 2I). Analysis of human lung cancer data sets derived from The Cancer Genome Atlas (TCGA) showed that KRAS, RAF1, MAP2K1, MAP2K2, MAPK3, and MAPK1 expression levels were positively correlated with the BCL6 expression level (Supplemental Figure 3). Collectively, our results suggest that the MAPK/ERK signaling axis promotes BCL6 expression in KRAS-mutant lung cancer cells.
ELK1 directly binds to BCL6 promoter and promotes its expression. Our results showed that mutant KRAS activated BCL6 expression through the MAPK/ERK signaling axis at the transcriptional level. Next, we characterized the mechanisms underlying this regulation. ELK1, c-Fos, JunD, c-Jun, and c-Myc have been reported to be critical transcription factors involved in the MAPK/ERK signaling cascade (25). Therefore, we investigated their effects on BCL6 expression after efficient silencing (Supplemental Figure 4). Surprisingly, we found that ELK1 knockdown significantly decreased BCL6 expression at both the mRNA and protein levels in H460 cells ( Figure 3, A and B). Data from luciferase reporter assays further showed that the BCL6 promoter activity was specifically suppressed after ELK1 silencing ( Figure 3C). Notably, ectopic expression of ELK1 in BCL6-depleted cells restored BCL6 expression ( Figure 3D), suggesting that BCL6 is dependent on ELK1 regulation. Moreover, BCL6 depletion decreased the clonogenic growth of H460 and H441 cells, whereas this effect was apparently reversed by simul-

Results
Mutant KRAS promotes BCL6 expression. The deregulated expression of BCL6, a master regulator of germinal center formation, is often linked to lymphomagenesis. Sustained BCL6 expression is also required to maintain malignant phenotypes in certain solid tumors (10). However, the role of BCL6 in lung cancer remains to be elucidated. Given that BCL6 elicits transcriptional activity and directly regulates multiple target genes, we hypothesized that BCL6 plays a role in lung cancer. Therefore, we determined BCL6 expression in an LSL-Kras G12D/+ mouse model of lung cancer. Strikingly, we observed a significant tumor burden accompanied by high BCL6 expression in mouse lung lesions 4 months after intranasal administration of adenovirus-Cre (Ad-Cre) (Supplemental Figure  1, A and B; supplemental material available online with this article; https://doi.org/10.1172/JCI161308DS1). Quantitative PCR (qPCR) and immunoblot analysis further showed that the mRNA and protein expression levels of BCL6 in mouse lung tumors harboring the KRAS(G12D) mutation were significantly higher than those in adjacent normal lung tissue ( Figure 1, A and B). Similar results of BCL6 induction were observed in mouse embryonic fibroblasts (MEFs) isolated from LSL-Kras G12D/+ mice ( Figure 1, C and D) and 2 pairs of engineered KRAS isogenic cell lines (Supplemental Figure 1C) upon KRAS mutational activation. Next, we examined BCL6 and KRAS-GTP protein expression in various human lung, colorectal, and pancreatic cancer cell lines and observed a positive correlation ( Figure 1E and Supplemental Figure 1D). To investigate the regulatory link between mutant KRAS and BCL6, we overexpressed different doxycycline-inducible G12 KRAS mutant variants (G12C, G12D, G12S, and G12V) in 293T cells, all of which commonly occur in human lung cancer. The mRNA and protein expression levels of BCL6 were robustly increased in the presence of activating KRAS mutations (Figure 1, F and G), reinforcing the general role of BCL6 in response to different KRAS mutant variants. We further showed that increased BCL6 abundance was observed in the nuclei of KRAS-mutant cells and was not due to altered protein stability (Supplemental Figure 1, E and F).
BCL6 is a known transcriptional repressor that regulates hundreds of target genes by recruiting different chromatin-modifying corepressor complexes (14,23). We next explored whether the KRAS oncogene enhances the global transcription of BCL6 target genes by performing RNA sequencing (RNA-Seq) in KRAS WT and mutant cells. Although some functions of BCL6 appear to be conserved between normal and malignant cells, BCL6 transcriptional programs would be altered in cells with different genomic profiles or in specific biological contexts (14,23). Given that BCL6 target gene sets have not been fully defined in solid tumors, we thus used comparative BCL6 target gene selection to identify the genes that were differentially expressed between KRAS WT and mutant cells. Our data revealed that the core of BCL6 target genes, including BCL6, was affected by the inducible activation of KRAS in engineered isogenic cells ( Figure 1H) and MEFs (Supplemental Figure 1G), supporting the notion that oncogenic KRAS regulates BCL6 signaling. The altered expression of representative, differentially expressed target genes, such as BCL6, NCOR1, TP53, CD69, PIM1, and PRDM1, was further confirmed by qPCR analysis ( Figure 1I and Supplemental Figure 1H). These results indicate that KRAS mutational activation leads to sustained expression of BCL6 and its target genes. correlated with BCL6 expression levels in human lung cancers (Figure 3H). Collectively, these data suggest that the MAPK/ERK/ELK1 signaling axis promotes BCL6 expression.
A relatively short region in the first noncoding exon1A of BCL6 (+254 to +300) has been characterized as its transcriptional binding site and is associated with high levels of BCL6 transcriptional activity (26,27). Thus, we performed in silico motif analysis using taneous ELK1 overexpression ( Figure 3E). To further validate the modulation of BCL6 expression by the KRAS/MAPK/ELK1 cascade, we genetically silenced KRAS and observed a potent and coincident reduction in p-ERK1/2, p-ELK1, and BCL6 expression (Figure 3F); additionally, ectopic expression of ELK1 in KRAS-depleted cells apparently restored BCL6 abundance ( Figure 3G). Analysis of TCGA database showed that ELK1 expression levels were positively  shown. Data in A, C, E, and K are expressed as mean ± SEM of 3 technical replicates. Statistical analyses in A, C, and E were performed using 1-way ANOVA with Tukey's multiple comparison test, and in K using unpaired 2-tailed Student's t test, *P < 0.05, **P < 0.01, ***P < 0.001. The immunoblots in B, D, F, and G were contemporaneous and run in parallel from the same biological replicate.
the JASPAR database (28) and identified an ELK1-binding motif in BCL6 noncoding exon1A region ( Figure 3I). As the transcriptional repressor BACH2 is directly controlled by ELK1 to sustain the T follicular helper cell phenotype and function (29), we used the reported BACH2 binding motif as the control. To show the direct binding of ELK1 to the BCL6 exon1A region, we performed electrophoretic mobility shift assays, where a 47-bp probe for the BCL6 exon1A region containing the ELK1-binding motif (WT probe) was incubated with purified ELK1 protein. Notably, the DNA-protein complex-shifted band was clearly observed in the presence of the BCL6 WT probe ( Figure 3J). In contrast, BCL6 probes with mutated bases or base depletion in CCGGAA motif reduced the binding activity to the ELK1 protein, suggesting the CCGGAA motif as the critical binding site. Using chromatin immunoprecipitation followed by qPCR analysis, we further confirmed the occupancy of ELK1 on the BCL6 exon1A region in H460 cells ( Figure 3K). Altogether, our results reveal a direct binding of ELK1 to BCL6, which is essential for BCL6 transcription.

BCL6 inhibition impedes the growth of KRAS-mutant cancer cells in vitro.
Given that the KRAS/MAPK/ELK1 axis directly regulates BCL6 expression, we next investigated whether activated BCL6 promoted KRAS-mutant lung cancer. First, we investigated the effects of BCL6 knockdown on tumor cell growth in vitro. We genetically silenced BCL6 in the HPNE cell pair ( Figure 4A) and found that BCL6 silencing elicited cytotoxic effects selectively toward KRAS-mutant cells, whereas their WT counterparts were minimally affected ( Figure 4A). We used a panel of human lung, colorectal, and pancreatic cancer cell lines and examined the effects of BCL6 knockdown on their survival. Although BCL6 depletion led to growth arrest in all tested cell lines, KRAS-mutant cancer cells were more susceptible to BCL6 inhibition than WT cancer cells ( Figure 4B and Supplemental Figure 5A). Furthermore, BCL6 knockdown dramatically reduced the clonogenic growth of various KRAS-mutant cancer cell lines ( Figure 4C and Supplemental Figure 5B). The reduced growth of H460 and H441 cells following BCL6 depletion could be rescued by BCL6 overexpression, suggesting a BCL6-dependent effect ( Figure 4D).
To test whether BCL6 induction creates an actionable vulnerability, we treated KRAS-mutant cancer cells with BCL6 chemical inhibitors (15,16,30,31). Our results showed that the 4 reported BCL6i elicited potent cytotoxicity in multiple KRAS-mutant cancer cell lines ( Figure 4E and Supplemental Figure 5C). When noted, COMP7, which targets the BTB/POZ domain of BCL6 and prevents partner binding (15), achieved the best therapeutic efficacy. BI-3802 also displayed antiproliferative activity in our study; however, its poor bioavailability limits its application in animal studies (32). Therefore, we selected COMP7 for further investigation. In agreement with the results of BCL6 genetic silencing, treatment with COMP7 resulted in much lower IC 50 values in KRAS-mutant cells than in WT cells, suggesting that BCL6 inhibition creates de novo vulnerability specific to KRAS-mutant cells ( Figure 4F). Moreover, COMP7 treatment attenuated the clonogenic growth of various KRAS-mutant cancer cell lines, eliciting a general inhibitory effect ( Figure 4G and Supplemental Figure 5D).
Single targeting of MEK was insufficient to inhibit KRAS-driven lung cancer, partly because of sustained ERK reactivation (33). Given that the MAPK/ERK/ELK1 axis regulates BCL6 expression, we additionally explored whether BCL6 pharmacological inhibitors, such as COMP7, in combination with MEKi, could lead to more potent tumor inhibition than single-agent treatment. To this end, we combined COMP7 with the MEKi trametinib and examined the combined effects on the viability of H460 and H441 cells by analyzing the drug combination index (CI). Our results showed that most CI values at diluted drug concentrations were less than 1 (Supplemental Figure 6), indicating a synergistic effect of dual targeting of MEK and BCL6.
BCL6 inhibition results in replication fork stalling and DNA damage. Thereafter, we investigated the underlying mechanisms by which BCL6 inhibition blocked KRAS-mutant cancer growth. Genome-wide RNA-Seq was performed in HPNE/KRAS cells and MEF isolated from LSL-Kras G12D mice (MEF/KRAS) upon BCL6 efficient silencing. We analyzed the enriched pathways of the differentially expressed genes between BCL6 silencing and the control groups. Our results showed that sets associated with cell cycle-related biological processes, such as E2F targets and MYC targets, were enriched and were in the leading rank of altered pathways (Supplemental Figure 7). Notably, the cell-cycle-promoting E2F target pathway was significantly downregulated after BCL6 silencing in both of the cell lines ( Figure 5A). A heatmap of shared E2F target signature genes across HPNE/KRAS and MEF/KRAS cells showed significantly decreased expression of preRC components ( Figure 5B), including hexameric minichromosome maintenance (MCM) 2, 4, 5, and 7, which form the core of replicative helicase at replication forks (34). Our qPCR results confirmed that BCL6 depletion significantly reduced the expression of MCM2, MCM4, MCM5, and MCM7, whereas ectopic BCL6 expression in BCL6-knockdown cells largely reversed these effects ( Figure 5C). Immunoblotting analysis showed similar results ( Figure 5D).
PreRC proteins bind to DNA in the G1 phase and mark all potential replication origins, a subset of which is activated in well-defined steps during the S phase (35). Since BCL6 silencing suppressed pre-RC expression, we deduced that BCL6 knockdown might reduce the binding of preRC components to DNA, leading to replication fork stalling. To test this possibility, we analyzed the cell cycle progression in KRAS-mutant lung cancer cells after BCL6 depletion alone or in combination with BCL6 reconstitution. As expected, MCM2 and MCM5 mRNA expression were significantly suppressed by BCL6 silencing in both G1 and S phases, and an inverse effect was observed after BCL6 overexpression ( Figure 5E). Furthermore, both genetic and chemical inhibition of BCL6 caused a significant decrease in DNA replication, as evidenced by DNA fiber analysis (Figure 5, F and G and Supplemental Figure 8, A and B).
MCM2-7 complexes are loaded onto DNA to license replication origins for use in the upcoming S phase. Inhibition of MCM loading is sufficient to induce DNA damage in eukaryotic cells (36,37). Next, we performed alkaline comet assays to measure cellular DNA integrity after treatment. Our results demonstrated that BCL6 silencing significantly increased DNA strand breaks, whereas this could be significantly decreased by BCL6 overexpression ( Figure 5H and Supplemental Figure 8C). Consistently, the BCL6i COMP7 showed similar effects ( Figure 5I and Supplemental Figure 8D). Meanwhile, the expression of γ-H2AX, a sensitive marker of DNA damage, increased after BCL6 depletion in both H460 and H441 cells ( Figure 5J). Importantly, the molecular The therapeutic efficacy of COMP7 in vivo was also tested using a patient-derived xenograft model of lung adenocarcinoma (LUAD) harboring the KRAS(G12V) mutation (LUAD-PDX/ KRAS) (41). Administration of COMP7 to mice bearing established LUAD-PDX/KRAS xenograft tumors significantly decreased tumor volume and weight ( Figure 7D). In addition, we observed a decrease in preRC gene expression and an increase in γ-H2AX expression in xenograft tumors after treatment (Figure 7, E and F), indicating that treatment induced replication fork stalling and DNA damage. COMP7 treatment was well tolerated at doses up to 50 mg/kg, as it did not cause significant changes in body weight and blood biochemistry parameters (Supplemental Figure 11, B and C). Our in vitro and in vivo pharmacological studies collectively demonstrate that BCL6 serves as a potential therapeutic target for the treatment of KRAS-mutant lung cancer.
Aberrant BCL6 expression is associated with worse survival outcomes. To uncover the link between BCL6 expression and clinical outcomes in KRAS-mutant lung cancer, we further examined BCL6 expression in human LUAD tissue, including KRAS-mutant subgroup (n = 62) and KRAS WT subgroup (n = 67). Immunohistochemical results demonstrated a much higher protein-expression level of BCL6 in the KRAS-mutant subgroup than that in the KRAS WT subgroup and adjacent normal tissue (n = 113; Figure 8A). The majority (89%) of KRAS mutant lesions showed positive staining for BCL6, of which 47% exhibited moderate or strong staining signals. In contrast, nearly all of the KRAS WT subgroup and adjacent normal lung tissue showed weak or negative staining for BCL6 ( Figure 8B), indicating a specific BCL6 expression pattern in KRAS-mutant LUAD.
Next, we tested the association between BCL6 protein expression and clinical outcomes using a commercial tissue microarray representing 41 KRAS WT and 33 KRAS-mutant LUAD samples. Survival analysis using Kaplan-Meier estimates showed that among KRAS WT LUAD, the BCL6-high expression population did not display a significant survival difference compared with the BCL6low expression population (P = 0.2239; Figure 8C). In contrast, increased BCL6 expression was significantly correlated with poor survival in KRAS-mutant LUAD patients (P = 0.0481; Figure 8D). We further determined the association between BCL6 expression and survival outcomes of patients with KRAS-mutant cancers (n = 689) using TCGA database. Our results demonstrated that tumor patients with higher BCL6 expression had a shorter 5-year survival rate than BCL6-low expression patients ( Figure 8E). Collectively, these results demonstrate that BCL6 expression correlates with poor patient survival in KRAS-mutant cancers, including LUAD.

Discussion
Mutant KRAS is a well-defined oncogenic driver in lung cancer. Although KRAS(G12C) and immune checkpoint inhibitors have been approved for the treatment of lung cancer (4), other KRAS-mutant subtypes lack effective targeted therapies. In our study, we identified BCL6, a master regulator of germinal-center response, as a promoting factor for KRAS-mutant lung cancer. BCL6 expression was profoundly upregulated by the KRAS/ MAPK/ELK1 axis and associated with poor patient survival. BCL6 maintained the MCM2-7 complex, which forms the core of replicative helicase to initiate DNA synthesis. BCL6 inhibition significantly triggered stalled replication forks and growth arrest, thereby regulation of BCL6 on preRC genes, such as MCM5 was functional and linked to clonogenic growth of KRAS-mutant lung cancer cells ( Figure 5K). It has reported that BCL6 overexpression attenuated while BCL6 knockdown enhanced the production of ROS levels in human vascular smooth muscle cells (38). As KRAS-mutated tumors are known to have high ROS levels and are more vulnerable to oxidative stress (25), it is possible that, in addition to DNA-replication regulation, the controlling effects of BCL6 on ROS levels might also contribute to the tumor-promoting effects of BCL6 in KRAS-mutant cancer cells. To test this possibility, we examined ROS levels upon genetic or pharmacological interference with BCL6 in H460 cells. Our results showed that BCL6 inhibition exerted little effect on ROS production in KRAS-mutant lung cancer cells (Supplemental Figure 9), thus excluding the role of ROS in BCL6-mediated tumor-promoting action in the context of KRAS mutational activation. Taken together, these data suggest preRC proteins as downstream targets of BCL6. Targeted inhibition of BCL6, therefore, leads to replication fork stalling, DNA damage, and ultimately, to growth arrest.
Lung-specific ablation of BCL6 inhibits KRAS-driven lung tumorigenesis. To assess the role of BCL6 in tumorigenesis and tumor growth in KRAS-driven lung cancer, we generated an engineered mouse model harboring the following alleles: (a) WT alleles of Bcl6 with LSL-Kras G12D/+ ; (b) heterozygous floxed alleles for Bcl6 with LSL-Kras G12D/+ ; and (c) homozygous floxed alleles for Bcl6 with LSL-Kras G12D/+ . The development of Kras-driven lung cancer in mice was monitored by noninvasive microcomputed tomography, as indicated in the experimental scheme ( Figure 6A). According to previously published studies (39,40), LSL-Kras G12D/+ mice developed lung adenomas 12 weeks after infection with Ad-Cre. Strikingly, our histological analysis showed that Kras activation with heterozygous or homozygous Bcl6 loss in lung epithelium was devoid of any hyperplasia or adenomas and even retained normal histology 12 weeks after Ad-Cre administration (Supplemental Figure 10). LSL-Kras G12D/+ mice developed lung adenomas and highly variable tumor latencies 20 weeks after infection, with a median survival of 131 days, whereas Bcl6 loss via heterozygous and homozygous recombination of floxed alleles significantly impeded tumor progression and prolonged mouse survival, with a median survival of 204 days and 276 days, respectively ( Figure 6, B-D). Histopathological analysis of lung tumors from moribund autochthonous mice revealed that the mice harboring Bcl6 floxed alleles displayed a significant decrease in the proportion of higher-grade adenomas or carcinomas ( Figure 6, E and F). Altogether, our data suggest that BCL6 is functionally required for KRAS-driven lung tumorigenesis.
Chemical inhibition of BCL6 impedes KRAS-mutant lung tumor growth in vivo. After characterizing the suppressive role of BCL6 knockout in KRAS-mutant lung tumorigenesis, we further investigated whether the pharmacological inhibition of BCL6 could exert similar effects. We evaluated the therapeutic efficacy of COMP7 in a LSL-Kras G12D/+ mouse model. LSL-Kras G12D/+ mice were randomly grouped 12 weeks after intranasal administration of Ad-Cre and then treated with COMP7 or vehicle for an additional 4 weeks (Supplemental Figure 11A). Our results showed that COMP7 administration significantly slowed lung tumor growth beginning as early as 2 weeks after administration (Figure 7, A and B) and prolonged mouse survival ( Figure 7C).  C, E, F, G, H, I, and K were expressed as mean ± SEM Statistical analyses in C, E, F, H, and K were performed using 1-way ANOVA with Tukey's multiple comparison test, and in G and I using unpaired 2-tailed Student's t test, *P < 0.05, **P < 0.01, ***P < 0.001. The immunoblots in D, J, and K were contemporaneous and run in parallel from the same biological replicate, representative of at least 3 independent experiments. particularly in the context of oncogenic drivers. Our data demonstrated that mutational activation of KRAS robustly promoted BCL6 expression (Figure 1), and genetic depletion of BCL6 exhibited a profound inhibitory effect on KRAS-driven lung tumorigenesis ( Figure 6), suggesting that BCL6 is a lynchpin for oncogenic KRAS dependence in lung cancer. Although BCL6-deficient mice expressing oncogenic KRAS still developed lung hyperplasia or adenoma, a dramatically reduced tumor incidence was observed upon heterozygous or homozygous loss of BCL6 in the lungs. This suggests that BCL6 may be a functional requisite in the early stages of lung tumorigenesis. In addition to its role in tumorigenesis, BCL6 maintains cell survival under various stress conditions (43).
blocking the malignant phenotype of KRAS-mutant lung cancer cells. Importantly, BCL6 depletion impeded lung tumorigenesis in a genetically engineered Kras G12D -driven mouse model, and pharmacological inhibition of BCL6 consistently reduced the growth of KRAS-mutant lung cancer cells and prolonged mouse survival. Our study highlights the promise and merits of targeting BCL6 to treat KRAS-mutant lung cancer.
BCL6 has been most thought of as a hematopoietic tumor-specific oncogene (42). Recent studies have also implicated BCL6 overexpression in certain solid tumors such as glioblastoma, breast cancer, and ovarian cancer (10). However, the link between BCL6 and the initiation and progression of lung cancer remains unclear, Statistical analyses in C and F were performed using 1-way ANOVA with Tukey's multiple comparison test, and in D using log-rank test, *P < 0.05, **P < 0.01, ***P < 0.001. en ELK1 occupancy of the BCL6 promoter is necessary for BCL6 expression and its biological function in promoting cell survival. In addition, BCL6-mediated silencing of multiple target genes has been reported in primary B cells and lymphoma (14), where BCL6 is typically involved in maintaining cell cycle progression and controlling DNA damage-sensing and proliferation checkpoints (12). However, few BCL6 target genes have been identified in solid tumors. During the cell cycle of the G1/S transition, the conversion of the MCM2-7 complex to an active DNA helicase leads to DNA synthesis. Therefore, inappropriate expression or assembly of the prereplicative complex has been linked to replicative stress and DNA damage, ultimately leading to premalignant dysplasia and cancer (36). A large-scale loss-of-function screen revealed that cancer cells expressing oncogenic KRAS are highly dependent on the transcription factor GATA2 and the DNA replication initiation regulator CDC6 (49). The KRAS/LKB1 oncological genotype imposes metabolic vulnerability dependent on DNA synthesis (50). Moreover, KRAS-mutant cancer cells are sensitive to suppression of the DNA replication licensing factor MCM7 (51). In our study, we identified E2F targetsincluding several preRC genesas downstream components of BCL6 in a KRAS-mutant Treatment with the tyrosine kinase inhibitor imatinib upregulated BCL6 expression, leading to transcriptional inactivation of the p53 pathway in Ph + acute lymphoblastic cells (44). Recently, we demonstrated that clinical BETi also increase BCL6 expression and proposed a model in which BRD3 maintained the autoregulatory circuit of BCL6 (22). Based on these findings, we assume that increased BCL6 expression may be an adaptive surge in the cellular response to oncogene-or nononcogene-induced stress, highlighting the multifaceted role of BCL6 in cancer and the broad application of BCL6-based therapy.
BCL6 is required for maintaining B cells in germinal-center compartments and maintains malignant phenotypes (17,45). However, the tumor-specific regulatory network of BCL6 remains elusive. The signal transduction protein STAT5 serves as a negative regulator of BCL6 in lymphomas (46), and FoxO3a promotes BCL6 expression in leukemia cells exposed to BCR-ABLi (44,47). The dependence of BCL6 expression on STAT1 has been observed in imatinib-treated chronic myeloid leukemia cells (48). While a limited number of transcription factors have been reported to regulate BCL6 expression, we revealed that ELK1 downstream of the MAPK/ERK pathway is a regulator of BCL6 (Figure 3). ERK-driv- Statistical analyses in B and D were performed using 1-way ANOVA with Tukey's multiple comparison test, and in C using log-rank test, *P < 0.05, **P < 0.01, ***P < 0.001. The immunoblots were contemporaneous from 4 independent samples and run in parallel from the same biological replicate.
ty of KRAS-mutant cells toward BCL6 suppression. In addition, several lines of evidence have shown that the combination of BCL6 with other therapeutic agents leads to enhanced tumor regression. The combination of a BCL6i with an Hsp90i or an HDACi was highly synergistic in vitro and yielded more potent antilymphoma effects in vivo (56). The fact that BCL6 and EZH2 cooperate to repress transcription leads to concurrent inhibition of EZH2 and BCL6 in diffuse large B cell lymphoma (57). Currently, MEKi-based therapy has failed to demonstrate significant survival benefits in KRAS-driven lung cancer, partly due to ERK reactivation (33). Since we identified BCL6 as what we believe to be an effector downstream of ERK, BCL6i may provide a potential means to block tumor cells with access to ERK reactivation. This holds the promise of concurrent targeting of BCL6 and other KRAS pathway effectors, such as MEK, to treat KRAS-driven lung cancer. However, further experiments are required to support the development of BCL6-based combination therapies for the treatment of this intractable disease.
Although mutant-selective KRAS(G12C) therapy has been clinically approved (52,53), most KRAS oncoproteins remain undruggable. Our results demonstrated that BCL6 was a crucial downstream effector of the KRAS/MAPK pathway. Notably, KRAS-induced BCL6 expression was independent of KRAS mutation status, as different G12 KRAS mutant variants (G12C, G12D, G12S, and G12V) that commonly occur in lung tumorigenesis could promote BCL6 expression to similar extents. Unlike KRAS G12C-selective inhibitors, targeted inhibition of BCL6 might be effective in the majority of KRAS-mutant cancer subtypes and even in tumors harboring NRAS mutations. Considering the high prevalence of RAS point mutations in human tumors and the responsive role of BCL6 in KRAS oncogenic signaling, BCL6 inhibition may improve antitumor activity as a widely applicable approach in cancer therapy. To execute transcriptional activity, BCL6 requires homodimerization and forms a complex with cofactors including BCoR, NCoR, and SMRT (54). BCL6 pharmacological inhibitors that block the interaction between BCL6 and BCoR/NCoR/SMRT exert cytotoxicity toward BCL6-addicted cancer cells (55). Our in vivo studies showed that inhibition of BCL6 transcriptional activity by COMP7 significantly impeded the growth of KRAS-mutant lung cancer cells in preclinical mouse models (Figure 7). Our findings revealed the vulnerabili- Cell cycle fractionation. Cell cycle fractionation was performed as described previously (36). Cells were seeded onto 6-well plates overnight and then treated with BCL6 siRNAs alone or in combination with pcDNA3.1-BCL6 for 48 hours. Cells were trypsinized and resuspended at a concentration of 1 × 10 6 cells/mL in complete media containing 10 mg/mL Hoechst 33342 stain (Beyotime Biotechnology) prior to sorting. The cell populations were fractionated into G1 and S phases based on Hoechst staining. Total RNA was isolated from each population using a RNAiso Plus Kit (TaKaRa Biotechnology).
DNA fiber assay. The DNA fiber assay was performed as previously described (58). Cells were plated at 50% confluence onto 15 cm plates and allowed to adhere overnight. The cells were treated with 10 μM COMP7 and BCL6 siRNAs alone or in combination with pcDNA3.1-BCL6 for 48 hours. The treated cells were subsequently incubated with 5-iodo-2′-deoxyuridine (IdU; I7125, Sigma) and 5-chloro-2′-deoxyuridine (CIdU; C6891, Sigma). After labeling, 2.5 μL of the cell suspension (approximately 2,500 cells) was spotted onto glass slides, followed by the addition of 7.5 μL lysis buffer. CIdU and ldU were detected by confocal microscopy using red (ab6326, Abcam) and green (347580, BD Biosciences) fluorescent secondary antibodies, respectively. The ratio of green and red fiber track lengths (IdU/ CIdU) was quantified using the ImageJ software.
Single cell gel electrophoresis (comet assays). Alkaline comet assays were conducted using a single-cell gel electrophoresis assay kit (Trevigen). Briefly, cells were seeded onto 6-well plates overnight and treated with COMP7 and BCL6 siRNAs alone or in combination with pcDNA3.1-BCL6 for 48 hours. Cells were diluted to 1 × 10 6 /mL in a solution of PBS and premelted LMAgarose (at 37°C) at a ratio of 1:10 (v/v) and pipetted onto CometSlide. Slides were immersed in lysis solution at 4°C for 2 hours and in freshly prepared alkaline unwinding solution (200 mM NaOH, 1mM EDTA; pH > 13) for 20 minutes at room temperature. Gel electrophoresis was performed using alkaline electrophoresis solution at 21 V for 30 minutes. The gels were kept at 37°C until the gel circle was completely melted and stained with 100 μL of 1 × SYBR Gold (TaKaRa Biotechnology) (diluted in TE buffer, pH 7.5, containing 10 mM Tris-HCl pH 7.5, and 1 mM EDTA). Fluorescent images were captured using an Olympus IX70 microscope and analyzed using the Casp software (https://sourceforge.net/projects/ casp/). Data are expressed as the percentage of tail intensity versus total DNA intensity.
Tumor induction and inhibitor administration in mice. Mice carrying the Bcl6 flox allele were obtained from Biocytogen Pharmaceuticals Co. Ltd. Conditional knockout of exons 5 and 6 of Bcl6 led to a truncated BCL6 protein that is not functional in vivo. Bcl6 flox/flox mice were crossed with LSL-Kras G12D mice (Jackson Laboratory) to generate LSL-Kras G12D/+ ; Bcl6 +/+ , LSL-Kras G12D/+ ; Bcl6 +/flox and LSL-Kras G12D/+ ; Bcl6 flox/flox mice on a C57BL6/N background. at 8 weeks of age, mice were infected with intranasal adenoviral-Cre (Ad-Cre; HanBio) at a dose of 2.5 × 10 7 PFU as previously described (59). At weeks 12, 16, 20, and 30 after infection, the lungs were imaged using a Quantum GX microCT imaging system (PerkinElmer Inc.). For inhibitor treatment, 8-week-old LSL-Kras G12D/+ mice were infected intranasally with Ad-Cre. Mice were randomly grouped 12 weeks after induction and treated with vehicle (0.5% CMC-Na in sterile water) or COMP7 (25 or 50 mg/kg, orally dissolved in 0.5% CMC-Na in sterile water) for an additional 4 weeks. At the end of the treatment, formalin-fixed lung lobes were bisected, sectioned, and stained with H&E. In survival studies, the mice were sacrificed when into H522 cells using an adeno-associated virus system (Addgene). HPNE and HPNE/KRAS cells were cultured in medium containing 1 volume of M3 Base culture medium (InCell) and 3 volumes of glucose-free DMEM (Thermo Fisher Scientific) supplemented with 10% fetal bovine serum (FBS, Thermo Fisher Scientific), 5.5 nM glucose, 10 ng/mL epidermal growth factor (R&D System), and 50 μg/mL gentamycin. 293T cells were maintained in DMEM supplemented with 10% FBS. Other cell lines were cultured in RMPI-1640 medium containing 10% FBS. KRAS mutational status of cancer cells was characterized by sequencing or genotyping, as shown in Supplemental Table 2. The cell lines were authenticated using short tandem repeat analysis and tested for mycoplasma contamination.
Luciferase activity assay. Full-length BCL6 and TP53 proximal promoter regions were obtained from 293T cell genomic DNA by PCR assays. The PCR product was subcloned into a pGL3 luciferase reporter vector (Promega). Cells were transiently transfected with 10 ng Tp53 or BCL6 promoter reporter constructs using Lipofectamine TM 2000 reagent (Thermo Fisher Scientific). Transfected cells were treated with the indicated compounds, siRNAs, or pcDNA3.1-BCL6 plasmids for 48 hours. Cells were lysed in 1 × passive buffer for 20 minutes at room temperature, frozen, and thawed once. Luciferase activity was performed as previously described (22). Cells transfected with modified pGL3 containing firefly luciferase, but lacking the 3′UTR, were used as control.
Tool compound screen. BCL6 expression levels were determined using immunoblotting and luciferase reporter assays. For immunoblot analysis, H460 cells were treated with different compounds targeting KRAS downstream effectors at a final concentration of IC 50 s or double IC 50 s for 48 hours. BCL6 protein levels were quantified using ImageJ software (NIH). The BCL6 expression level (xi) for each compound was normalized to that of the DMSO-only control. The readout of each compound was transformed into a z score, as follows: z score = (xiμ)/σ, where μ represents the mean BCL6 expression level of each compound in a given plate, and σ is the SD of the BCL6 expression level. The absolute value of the z score represents the degree of BCL6 expression-level variation induced by a given compound. For luciferase reporter assays, H460 cells were treated with different compounds at a final concentration of IC 50 s for 48 hours. The readout of treated cells was normalized to that of DMSO-only control cells.
ChIP-qPCR analysis. The ChIP assays were performed using a SimpleChIP Plus Enzymatic Chromatin Immunoprecipitation Kit (Cell Signaling Technology) according to the manufacturer's instructions. Cells were collected and digested with micrococcal nuclease for 30 minutes, generating an average chromatin fragment size of 150-900 bp. Chromatin (4 × 10 6 cells) and 5 μg of ELK1 antibody (ab125085; Abcam) were used for ChIP assays. Isotype IgG was used as a negative control. The relative enrichment of BACH2 and BCL6 was quantified using qPCR assays. Primers used are listed in Supplemental Table 4.


Editorial: Advanced Strategies for the Recognition, Enrichment, and Detection of Low Abundance Target Bioanalytes

Biosensors and bioassays provide important analysis capabilities for biological, medical, and environmental applications. Owing to active and creative exploration of emerging materials, devices, and methods, there have been great advances in biosensors and bioassays, making detection more sensitive, rapid, accessible, and intelligent. For biodetection, especially for low abundance, strategies for target recognition, enrichment, and detection are of vital importance. The goal of this Research Topic is to provide a forum sharing recent research and novel ideas on related techniques. We have collected eight papers focused on biosensors for target or event detection: three research articles, two opinion articles, and three reviews. The techniques include microfluidics, electrochemistry, CRISPR, emerging materials, and so on. In terms of research papers, Zhu et al. report a microfluidic device integrated with an electrical impedance spectroscopy (EIS) biosensor to perform in-situ impedance measurement of yeast proliferation in single-cell resolution. Hydraulic shear force is utilized to detach a daughter cell from its mother cell, and time-lapse impedance measurement is performed to monitor the cellular process. The main highlight of this work is a combined sensor achieving in-situ and real-time monitoring of a dynamic event for single-cell reproduction, which uses simple EIS instead of image observation, saving hardware resources for image recording. Meanwhile, the capacity to monitor for longer than several hours still needs to be improved. Salvador et al. present the development of a fluorescent microfluidic device based on an antibody microarray for therapeutic drugmonitoring of acenocoumarol (ACL). A fully integratedmicrofluidic system is realized using a specific antibody for ACL on the glass slide in a microarray chip, with a fluorescent label for target binding. A limit of detection (LOD) of 1.23 nM is achieved in human plasma. As a proof-of-concept point-of-care device, this system is automatic and sensitive, although the exact turnaround time is not mentioned. The selectivity also needs to be demonstrated for situations such as the presence of other molecules similar to ACL existing in the blood. Yuan et al. have developed a one-step electropolymerized biomimetic polypyrrole membranebased electrochemical sensor for selective detection of valproate (VPA). The critical point of this strategy is to fabricate molecularly imprinted polymer (MIP) with simple one-step electropolymerization and employ the MIP/gate effect to realize sensitive concentration-electrical signal response. This sensor has a LOD of 17.48 μM and good selectivity against five drugs in Edited and reviewed by: Guozhen Liu, The Chinese University of Hong Kong, Shenzhen, China

Edited and reviewed by:
Guozhen Liu, The Chinese University of Hong Kong, Shenzhen, China combination therapy with VPA. This work furthers a step toward conveniently controlling VPA concentrations for patients, and the main shortfall is the lack of confirmatory tests in practical blood or serum, which is a much more complex matrix than deionized water. For opinions, Jiang et al. introduce the working modes and applications of biosensors for point mutation detection, including the biosensors for RNA mutation rate detection of SARS-CoV-2. As an alternative to traditional methods, biosensors for gene mutation detection are fast, low-cost, and suitable for large-scale applications. Using ssDNA probes for hybridizing with the point mutations, gene sequencing, and PCR-based methods become unnecessary. Meanwhile, versatile signal conversion approaches enable various sensitive sensors. As a result, the LOD for point mutation can be as low as 0.5 fM. ssDNA is commonly designed as an oligonucleotide for biomolecule recognition, while the basic function of hybridization promotes the detection of point mutation.
Jiang et al. describe a possibility for the application of artificial intelligence-based triboelectric nanogenerators (AI-TENGs) in biosensor development. The mechanism of triboelectric nanogenerators and their combination with artificial intelligence are introduced. Representative research on AI-TENG for biomolecules sensing and other biotechnology is then highlighted. TENG as an energy collection and conversion technique is of great value for portable biosensors needing a power supply, while more investigations might be needed on the integration of AI and TENG for biosensor development.
In terms of the reviews, Li et al. review recent advances in metal-organic framework (MOF)-based electrochemical biosensing applications. The main content includes the synthesis and modification of MOFs for electroactive materials and emerging biosensing applications of these MOF materials. The emphatically introduced applications are small-biomolecule sensing, biomacromolecule sensing, and pathogenic cell sensing. The specific improvement to figures of merit by MOFs, such as LOD, is not presented in this paper, and oxidative degradation in certain solutions might be considered when MOFs are applied in electrochemical sensors.
Liu et al. discuss nanomaterial-based immunocapture platforms for recognition, isolation, and detection of circulating tumor cells (CTCs). Three parts are outlined: design principles for immunoaffinity nanomaterials, nanomaterial-based platforms for CTC immunocapture and release, and recent advances in single-cell release and analysis for CTCs. Nanomaterials of different shapes, sizes, and structures are introduced, showing their good performance for agent linking and CTC isolation and detection. In this article, capture efficiency is highlighted as an important figure of merit. Although the overall detection time of most mentioned techniques is not introduced, this review remains a comprehensive review of CTC separation and detection.
Chen et al. have created an overview of the recent development of clustered regularly interspaced short palindromic repeats (CRISPR)-based biosensing techniques and their integration with microfluidic platforms. A representative microfluidic CRISPR sensor achieves a LOD of 20 pfu/ml of purified Ebola RNA within 5 min, and most similar sensors are all accompanied by high sensitivity and fast responses.
Working with the two powerful techniques, these biosensors, being applied in nucleic acid-based diagnostics, protein tests, metal ion monitoring, and protein/small molecule interactions screening, have great advantages for point-of-care testing of various bio-targets with low abundance.
In summary, the articles collected in this Research Topic include various emerging materials, methods, strategies, and platforms promoting technological innovation for recognition, enrichment, and detection of low abundance target bioanalytes. In these articles, microfluidic devices are most utilized for target collection and enrichment, while electrochemical biosensors are widely adopted transducers for signal conversion. Advanced and interdisciplinary fields such as CRISPR, MOF, and AI-TENG are represented in these articles.
For the goal of trace bio-target analysis, the decisive techniques are target enrichment and sensing mechanism, possibly including signal amplification. Efficient target enrichment and sensitive signal transformation remain critical challenges to be solved, and efforts are being made to seek advanced strategies besides the methods mentioned in this Research Topic. As a promising strategy, AC electrokinetics is demonstrated to be a simple and versatile approach for rapid bio-target enrichment (Salari and Thompson, 2018). Meanwhile, researchers have developed an ultrasensitive sensing mode by monitoring the solid-liquid interfacial capacitance of an electrode array. Combining the two mechanisms, the LOD level of 0.1 fM (Zhang et al., 2021), fg/mL (Qi et al., 2022), or 10 CFU/ml (Zhang et al., 2020) is achieved when ions, proteins, or bacteria are detected with the turnaround time less than 1 min. With technological advances being continuously on the horizon, more sensitive, convenient, rapid, and intelligent detection is expected to emerge in the future.

AUTHOR CONTRIBUTIONS
JZ and JJW wrote the editorial, which was proofread and accepted by all the authors.


Comparative Analysis of Codon Bias in the Chloroplast Genomes of Theaceae Species

Theaceae species are dicotyledonous angiosperms with extremely high ornamental and economic value. The chloroplast genome is traditionally used to study species evolution, expression of chloroplast genes and chloroplast transformation. Codon usage bias (CUB) analysis is beneficial for investigations of evolutionary relationships and can be used to improve gene expression efficiency in genetic transformation research. However, there are relatively few systematic studies of the CUB in the chloroplast genomes of Theaceae species. In this study, CUB and nucleotide compositions parameters were determined by the scripts written in the Perl language, CodonW 1.4.2, CU.Win2000, RStudio and SPSS 23.0. The chloroplast genome data of 40 Theaceae species were obtained to analyse the codon usage (CU) characteristics of the coding regions and the influence of the source of variation on CUB. To explore the relationship between the CUB and gene expression levels in these 40 Theaceae plastomes, the synonymous codon usage order (SCUO) and measure independent of length and composition (MILC) values were determined. Finally, phylogenetic analysis revealed the genetic evolutionary relationships among these Theaceae species. Our results showed that based on the chloroplast genomes of these 40 Theaceae species, the CUB was for codons containing A/T bases and those that ended with A/T bases. Moreover, there was great commonality in the CUB of the Theaceae species according to comparative analysis of relative synonymous codon usage (RSCU) and relative frequency of synonymous codon (RFSC): these species had 29 identical codons with bias (RSCU > 1), and there were 19 identical high-frequency codons. The CUB of Theaceae species is mainly affected by natural selection. The SCUO value of the 40 Theaceae species was 0.23 or 0.24, and the chloroplast gene expression level was moderate, according to MILC values. Additionally, we observed a positive correlation between the SCUO and MILC values, which indicated that CUB might affect gene expression. Furthermore, the phylogenetic analysis showed that the evolutionary relationships in these 40 Theaceae species were relatively conserved. A systematic study on the CUB and expression of Theaceae species provides further evidence for their evolution and phylogeny.

INTRODUCTION
Theaceae, a characteristic family of subtropical evergreen broadleaved forests, is widely distributed in subtropical regions . Different Camellia (Camellia sinensis L. Kuntze.) varieties are popular in many countries around the world, including India, Vietnam and Japan (Teixeira and Sousa, 2021;Wang et al., 2021). Camellia originated in Southeast Asia, and 80% of the species of this genus come from China. Theaceae species have high ornamental value and high economic value, such as medicine, food and drink (Yuan et al., 2021). Tea is water soaked with Camellia leaves, known as "oriental magic water", which is one of three most famous drinks in the world (Balentine et al., 1997). In addition, the seeds of Camellia are rich in oil that can be extracted and used for various purposes, such as edible and industrial oils.
The chloroplast genome usually consists of double-stranded closed circular DNA molecules, but the chloroplast DNA of a few plants is linear. The chloroplast genome is relatively small, generally ranging from 107 to 218 kb in length, accounting for approximately 10%-20% of total plant DNA (Yao et al., 2015). The chloroplast not only has a relatively stable evolutionary rate but also has a relatively slow rate of nucleotide substitution. Thus, the investigation of the chloroplast genome is of great significance for revealing the evolutionary rules, genetic relationships, expression of chloroplast genes and chloroplast transformation among plants (Olmstead and Palmer, 1994).
Codons are the link between amino acids, proteins, and genetic materials in organisms. Meanwhile, codons play an irreplaceable role in transmitting genetic information. There are 64 codons that collectively code for 20 amino acids. Except for methionine (Met, coded by AUG) and tryptophan (Trp, coded by UGG), which have unique codons, the remaining amino acids generally have 2-6 corresponding codons; this phenomenon is known as codon degeneracy (Ikemura, 1985). Codons coding the same amino acid are called synonymous codons, and their frequencies of use in the translation process are different. This phenomenon of unbalanced codon use is called codon usage bias (CUB; Sharp and Li., 1987;Wang et al., 2020). The study of CUB is beneficial for exploring genetic evolution, understanding gene expression characteristics and providing guidance for molecular breeding.
In recent years, the rapid development of high-throughput sequencing technology has led to increasingly enriched chloroplast genome databases, laying a solid data foundation for plant genome structure and phylogenetic analyses. Thus, the chloroplast genome has been studied more extensively and indepth. Moreover, related research on CUB in the chloroplast genome has been carried out in a variety of plants, including Amaranthus hybridus (Chaney et al., 2016), Camellia oleifera (Zhang et al., 2019), Medicago truncatula (Guo et al., 2020), Nicotiana tabacum (Arévalo-Gallegos et al., 2020); meanwhile, Magnoliaceae (Chang et al., 2006), Bambusoideae (Zhang et al., 2011), Euphorbiaceae plants (Wang et al., 2020) and Oryza species (Chakraborty et al., 2020). There are numerous reports on the chloroplast genome of the Theaceae species, but relatively rare on its CUB.
From the perspective of the coding sequence, we systematically analysed the codon usage (CU) characteristics of the chloroplast genomes of Theaceae species, the source of variation, the highfrequency (HF) codons, and the relationship between the CUB and gene expression levels based on measures independent of length and composition (MILC) and synonymous codon usage order (SCUO) parameters, and phylogenetic analysis was used to investigate the genetic relationships between the plants. This study revealed the relationships between the CUB and gene expression frequency, and provided a reference direction for further genetic evolution research and transgenic research in these 40 Theaceae species.

Source of Material
The coding regions (coding sequences, CDSs) of the chloroplast genomes of 40 Theaceae species, namely, 33 species of Camellia, one species of Hartia (Hartia laotica) and six species of Stewartia, were downloaded and collected from the National Center for Biotechnology Information (NCBI) database (https://www.ncbi. nlm.nih.gov/) on May 17, 2021. Detailed information on the chloroplast genomes of the 40 Theaceae species is provided in Table 1. To ensure the accuracy of the experimental results, the scripts written in the Perl language (Wang et al., 2020) were used to process the original CDSs of the 40 chloroplast genomes, and the selected sequences were required to meet the following conditions: 1) the CDS is composed of only A, T, G, and C bases; 2) the sequence length is a multiple of three; 3) the sequence starts with the start codon ATG and ends with the stop codon TAA, TAG or TGA; 4) there is no stop codon in the middle of the sequence (means to exclude prematurely terminated coding sequences); 5) the length of the CDS is ≥ 300 bp (He et al., 2016;Li et al., 2016;Wang et al., 2020).
Coding Usage Index GC1, GC2, and GC3 respectively refer to the GC content of the first, second, and third base positions, respectively, of the chloroplast codons in the 40 Theaceae species. The GC1, GC2, and GC3 contents and the average GC content among the three sites of the 40 Theaceae species were calculated by using the scripts written in the Perl language (Wang et al., 2020). The CU index of the chloroplast CDSs was analysed by the CodonW 1.4.2 program, including T3s, C3s, A3s, and G3s, representing the frequency of T/C/A/G usage at the third base of the codon. The GC content was calculated as a percentage of the total base content, and the RSCU was determined. positive bias, as they are used more frequently than any other synonymous codons; while an RSCU value < 1 represents negative bias, indicating that these codons are used less frequently than any other synonymous codons; an RSCU value = 1 indicates no bias, which means that synonymous codons are used equally and that CU is random (Sharp and Li, 1986). For the calculation formula of RSCU, please refer to the following formula (Eq. 1) by Sharp and Li. (1986): A Toolkit for Biologists integrating various biological datahandling tools (TBtools) software was used to create a heatmap of the average RSCU values (Chen et al., 2020).
RFSC refers to the ratio of the number of codons observed in the plants to the total number of synonymous codons. The calculation formula for the frequency of use under the influence of gene length and amino acid abundance is detailed in the study of Zhou et al. (2007a) and Zhou et al. (2007b). The RFSC value can be calculated according to the actual number of codons using the CodonW 1.4.2 program (Sharp and Li, 1986). RFSC can be calculated using the following Eq. 2 (Sharp and Li, 1986): When the RFSC value of a codon is greater than 60% or more than 1.5 times the average frequency of its synonymous codons, this codon is called an HF codon (Zhou et al., 2007a;Zhou et al., 2007b).

Source Analysis of Codon Variation
Correspondence analysis (COA) multivariate statistical methods are often used to study the variance and gene distribution of RSCU, and each CDS was represented as a 59-dimensional vector (excluding the codon ATG encoding Met, the codon TGG encoding Trp, and the three stop codons TAA, TAG, and TGA) (Sharp and Li, 1986). The COA in CodonW 1.4.2 is generally used to study the CU variation in the genome. To reduce the influence of amino acid composition on codon utilization, each dimension corresponds to the RSCU value of each relevant codon, eliminating the influence of the stop codons TAA, TAG, and TGA. The effective number of codons (ENc) evaluates codon use bias at the genomic level; the lower the ENc value is, the more significant the codon use preference is (Wright, 1990;Zhang et al., 2020).
In the ENc-GC3s plot, GC3s and ENc values were used as horizontal and vertical coordinates, respectively, to analyse the relationship between the base composition and CUB (Wang and Hickey, 2007).
Parity rule 2 (PR2) plot analysis is widely used to explain whether there is an impact on nucleotide composition due to mutation pressure and selection pressure in DNA double strands, which includes G3/(G3+C3) and A3/(A3+T3) as the horizontal and vertical coordinates, respectively. The positions of A = T and G = C (PR2) are the centres, and the coordinates are (0.5, 0.5; Meyer, 2021). The centre position of the graph indicates that there is no deviation of mutation and selectivity (replacement rate) between two complementary DNA strands. When the proportions of G and C (or A and T) are similar, the gene CUB is affected only by mutation pressure; if the composition ratio is too large, it indicates that the CUB originates from natural selection and other factors (Sueoka, 1999).
Neutral analysis is a method for quantitatively analyzing the influence of directed mutation pressure and natural selection on CUB. We examined the neutral plot and calculated the GC12 and GC3 values of the genome using the scripts written in the Perl language (Wang et al., 2020). When the slope of the regression curve goes to 0 and there is no significant correlation between GC12 and GC3, it is fully influenced by natural selection. In contrast, when the slope is close to or equal to 1 and the correlation is significant, mutation pressure is speculated to have a major effect on genes (Sueoka, 1988).

SCUO Analysis
To better illustrate the relationship between the CUB and gene expression, we used the CU.Win2000 program to calculate the SCUO of the chloroplast sequences of the 40 Theaceae species (Wan et al., 2006). The SCUO value represents the CUB of synonymous codons relative to the whole sequence. SCUO values range between 0 and 1. The larger the value is, the stronger the CUB is. The calculation formula is explained in detail in the study by Wan et al. (2006).

MILC Analysis
MILC is a measure of the gene expression level used to quantify the difference in CU between a gene and some expected distribution of codons. The codon distribution can either be calculated from the background nucleotide composition or derived from a single gene or a gene group (Supek and Vlahovicek, 2005). For each chloroplast genome sequence, CU analysis in R uses the package coRdon 1.13.0. First, we need to count occurrences of each codon in each sequence. The default parameter values are set during the calculation, and the MILC value of each sequence in the set is calculated for the average CU deviation of the entire sample (https://bioconductor.org/packages/devel/bioc/vignettes/coRdon/). A high MILC value indicates a high gene expression standard, and a low MILC value indicates a low gene expression standard (Dorić and Bilela, 2014).

Phylogenetic Analysis of the 40 Theaceae Species
Two codons (AUG and UGG) and three stop codons (UAA, UAG, and UGA) encoding only one amino acid were removed from the chloroplast genomes of the 40 Theaceae species. The RSCU values of the remaining 59 synonymous codons were used as variables to calculate the distance between genes, namely, the squared Euclidian distance (Hou et al., 2016). Then, to explore the relationships between the CUB of the Theaceae species and their phylogenetic evolution in association with species relatedness, the intergroup linkage method in the SPSS 23.0 software (Pallant, 2013) was used to draw a cluster diagram for cluster analysis. The chloroplast genome sequences of 40 Theaceae species were aligned using the Muitiple Sequence Alignment Program (MAFFT v.7.4.0) (Katoh and Standley, 2013). The Maximum Likelihood phylogenetic tree of 40 Theaceae species was constructed with RAxML v.8.2.12 (Stamatakis, 2014) with 1000 bootstrap replicates, referring to the method of Liu et al. (2021). Other parameters were set by default in the software.

Nucleotide Composition Analysis
The chloroplast coding genes of the 40 Theaceae species had very similar base compositions in the whole coding region ( Table 2). The GC1, GC2 and GC3 contents were 35.84%-39.61%, 35.95%-39.88% and 35.61%-39.65%, respectively, and the mean values of the three were 37.14%-38.29%, which indicates an overall preference for codons containing A and T and ending in A/T. Therefore, from the perspective of base composition, the CUB of the chloroplast genomes of these 40 Theaceae species is extremely similar.

RSCU, RFSC Analysis and HF Codon Screening Reveal Biased Use of A/T
As shown in Figure 1, the number of preferred codons (RSCU > 1) in the 40 chloroplast genomes was 29, among which, 28 (96.55%) ended in A or U (excluding UUG, 12 ended in A and 16 ended in U), the probability of the preferred codon ending in A/T, was much higher than that ending in G/C. In all the Frontiers in Genetics | www.frontiersin.org March 2022 | Volume 13 | Article 824610 Theaceae species, AGA, encoding arginine (Arg), was strongly preferred (RSCU > 1.7), followed by GAU, encoding aspartic acid (Asp). The RFSC value was calculated by using the actual number of codons used by the CodonW 1.4.2 program (Supplementary Table S1). In summary, the codons in the chloroplast genomes of the 40 Theaceae species all showed a preference for ending in A/U. Furthermore, in Figure 1, the colours of the 40 Theaceae species for the same codon were extremely similar, indicating that the RSCU values for a given codon in the chloroplast genomes of the Theaceae species were similar. The HF codons in the chloroplast genomes of the 40 Theaceae species are shown in

Analysis of Codon Variation Suggests Natural Selection Dominates
In the COA of the 40 Theaceae species (Figure 2), the points representing different genes were separated from each other in all plants of Camellia. In the chloroplast genomes of the 40 Theaceae species, axis 1 and axis 2 accounted for 9.32%-10.73% and 8.05%-9.95% of the synonymous CU variation, respectively. Genes with different GC contents are marked with different colours to reveal the influence of GC content on CUB. The GC content of all the genes in the 40 Theaceae species was less than 45%. This high consistency indicates that the sources of CU variation in different Theaceae species are extremely similar. In the ENc-GC3s map of the 40 Theaceae species (Figure 3), all the genes were scattered in small clusters on the left side. Based on observation and analysis, the chloroplast genomes of the 40 Theaceae species showed scattered small clusters overall. Most of the genes were distributed far away from the standard curve and others were distributed in the appendix of the standard curve. This shows that the CUB of the chloroplast genomes is affected not only by mutation pressure but also by natural selection and other factors.
Natural selection may affect the use of bases A and T together with G and C. In this study, we found that the chloroplast genes of the 40 Theaceae species were distributed unevenly in the four quadrants of the PR2 plot (Figure 4), and they were distributed mainly in the regions of G3/(G3+C3) > 0.5 and A3/(A3+T3) < 0.5, indicating that the use of the third codon was unbalanced. Figure 5 shows an important positive correlation between GC12 and GC3 of the chloroplast genomes (P ≤ 0.01), which indicates that all the Theaceae species have been subject to direct mutation pressure on chloroplast gene CU. In the neutral analysis chart, the distribution ranges of GC12 and GC3 were relatively narrow. GC12 was distributed between 0.32 and 0.57, and GC3 was  distributed between 0.14 and 0.36. Most points were distributed in small clusters, indicating that both natural selection and mutation pressure affect CUB.
The neutral analysis chart can further reveal the strength of mutation pressure and natural selection. When the slope of the regression line is close to 0, natural selection plays a leading role,  and when the slope of the regression line is close to 1, rapid mutation pressure is assumed to play a main role. Among six Camellia species, namely, Camellia amplexicaulis, Camellia brevistyla, Camellia fascicularis, Camellia gauchowensis, Camellia granthamiana, and Camellia leptophylla, the absolute value of the angle at which the regression line slopes was slightly     Table 1.
Frontiers in Genetics | www.frontiersin.org March 2022 | Volume 13 | Article 824610 13 most genes in the chloroplast genome of Theaceae species have a low CUB. In summary, from the perspective of SCUO, the CUB of the chloroplast genomes of the 40 Theaceae species is very similar.

CUB Affects Gene Expression
The MILC value reflects the level of gene expression. The greater the MILC number is, the higher the gene expression levels is. We calculated the MILC value of chloroplast genes of  Table S2), and the results showed that 31 Theaceae species had a MILC value of 0.55, and the other 9 species had a MILC value of 0.56, which revealed that the expression level of most of the chloroplast genes was moderate. Additionally, we used correlation analysis to explore the relationship between the gene expression level and the CUB. We observed a significant positive correlation between the SCUO and MILC values ( Table 4), indicating that the gene expression levels may be affected by the CUB.

Phylogenetic Analysis of the 40 Theaceae Species
The phylogenetic tree constructed based on the cluster analysis of the RSCU values and the ML analysis of the chloroplast genomes has a highly similar topology. The greater the squared Euclidean distance between plants is, the greater the difference in their CUB. According to the RSCU calculations, the squared Euclidean distance between the chloroplast genes of the 40 Theaceae species was extremely small (all less than 0.12). At the same time, the clustering results divided the plants into two major branches: the first branch contained ten species of Camellia, all belonging to the genus Camellia; the second branch contained 30 of Theaceae, including 23 species of Camellia, one species of Hartia (Hartia laotica) and six species of Stewartia species (Figure 6). On the first branch, Camellia brevistyla, Camellia gauchowensis and Camellia yuhsienensis formed a small group, while the remaining plants belonged to another group: Camellia azalea, Camellia petelotii, Camellia grandibracteata, Camellia pubicosta, Camellia leptophylla, Camellia Nitidissima, and Camellia reticulata. On the second branch, two species of Camellia (Camellia monadelpha and Camelliaserrata), one species of Hartia (Hartia laotica) and six species of Stewartia formed a small group, with the two species of Camellia (Camellia renshanxiae and Camellia sasanqua) forming their own small group. The other 21 species of Camellia belonged to another small group. The results showed that the 21 species of Camellia in the second clade were more closely related, while the 10 species in the first clade were more distantly related.
The results of the constructed phylogenetic tree based on the genomes are samely divided into two major branches: the first branch contains 33 species of Camellia and all belong to the genus Camellia, and the second branch contains the remaining seven species of Camellia, including one species of Hartia (Hartia laotica) and six species of Stewartia. However, the six species of Stewartia belong to a small branch in phylogenetic tree according to the RSCU values. Overall, the topological structure of the phylogenetic tree and the clustering results based on the RSCU values of the chloroplast genome are similar.

DISCUSSION
Codons are an important core element that connects amino acids, proteins and genetic materials in living organisms (Sundararajan et al., 2016). Studies of the preferred usage bias of codons provide reliable information for the study of protein expression and its corresponding functions (Wang et al., 2018). Many biological factors have been demonstrated to influence CUB, including gene length, tRNA quantity, gene expression level, chain-specific mutation bias, and GC composition (Wan et al., 2004;Meyer, 2021).
Therefore, this study systematically analysed the CUB characteristics, CUB variation sources and phylogenetic relationships of the chloroplast genomes of 40 Theaceae species. As claimed by the results of RSCU, RFSC analysis and HF codon screening, the codons tended to end with A/T, which is consistent with the assumption that "higher plants tend to use codons ending with A/T" (Campbell and Gowri, 1990;Wang et al., 2020). According to Sharp et al. (1988) proposed that "organisms with close genetic relationships have very similar codon usage bias", the CUB in the chloroplast genomes of the 40 Theaceae species are extremely similar based on the RSCU values and comparative analysis of CU. In addition, a CUB analysis of 18 rice chloroplast genomes revealed that their CUB was low and AT base preference was strong; different rice varieties have different CU patterns that are mainly subjected to mutation pressure and natural selection (Chakraborty et al., 2020). We also analysed the preferred codons and HF codons of the 40 Theaceae species   (Campbell and Gowri, 1990); specifically, there was very high similarity in the use of codons among the 40 Theaceae species. Eighteen codons were both preferred codons and HF codons: GCU, UGU, GAU, GAA, UUU, GGA, CAU, AUU, AAA, UUA, AAU, CCU, CAA, AGA, UCU, ACU, GUA, and UAU. Compared to the use of one type of codon (preferred codons or HF codons), the use of the two types (preferred codons and HF codons) may have a higher reference value for the optimization of codons in terms of heterologous expression.
Regarding the source of codon variation, our results further suggested that both mutation pressure and natural selection affect CUB in the chloroplast genomes of the 40 Theaceae species. ENc-GC3s and PR2 plot analysis showed that the CUB of the 40 Theaceae species was influenced by mutation pressure and natural selection, and the similar results were found in other tea plants (Rawal et al., 2020;Rawal et al., 2021). In addition, a natural selection study of the chloroplast genomes of 29 species of Magnoliaceae systematically compared their CU characteristics and elucidated their evolutionary relationships (Chang et al., 2006). Meanwhile, Zhang et al. (2011) get the similar conclusion on 3 species in the bamboo subfamily. Comparisons of plant chloroplast genomes have revealed the main reason for their CUB. Furthermore, the quantitative results of neutral analysis showed that the CUB of the 40 Theaceae species were mainly affected by natural selection, which was consistent with the influencing factors of CUB in the chloroplast genomes of Magnoliaceae (Chang et al., 2006) and Euphorbiaceae plant ( Wang et al., 2020). In summary, natural selection plays a leading role while mutation pressure plays a secondary role in determining the CU of chloroplast genes in 40 Theaceae species.
In our study, the calculated SCUO values were lower than 0.25, which means that most genes in the chloroplast-encoding genome of Theaceae plastomes have a low bias. This result is consistent with the result of chloroplast genes in Oryza species (Chakraborty et al., 2020). In Camellia amplexicaulis, Camellia fascicularis, Camellia granthamiana, Camellia japonica S288C, Camellia kissii, Camellia ptilophylla, Camellia reticulata, Hartia laotica, Stewartia micrantha, Stewartia monadelpha, Stewartia obovata, Stewartia serrata, Stewartia sinii, and Stewartia villosa, the average value of SCUO was 0.23, which reflects the low codon bias. In the results of 18 Oryza species studies, the SCUO values ranged from 0.24 to 0.27, indicating that the CUB was low and varied in most of the rice species. Likewise, the 40 Theaceae species had low chloroplast CUB. We regarded MILC as a measure of gene expression in our current study and observed high MILC values. The MILC value of 31 of the Theaceae species was 0.55, and the value of the other nine kinds of MILC was 0.56. In the related research on Oryza species, the expression level of each chloroplast gene of 18 rice species based on the MILC value was calculated to be in the range of 0.74 to 0.78, which means that most of the chloroplast genes are expressed at high levels in rice and that there are differences in gene expression levels (Chakraborty et al., 2020). The results showed that, from the perspective of MILC, there may be certain differences in gene expression between dicotyledonous plants and monocotyledonous plants of chloroplast genomes. Karl Pearson's correlation analysis between SCUO and MILC was used to explore the effect of CUB (Sareen et al., 2018). Based on the gene expression level, we observed a significant positive between SCUO and MILC values (P < 0.01; Table 4), indicating that the gene expression levels may be affected by the CUB. All in all, the higher the CUB is, the higher the gene expression level is (Li et al., 2016). In related reports by Zhao et al. (2016) and Liu et al. (2020), there is a close correlation between gene expression frequency and CUB, that is, the stronger the CUB, the higher the gene expression frequency. The efficiency of gene expression can be improved by codon optimization. For example, enhancing the expression of viral genes in mammalian cell lines by complete resynthesis of the gene (Gustafsson et al., 2004). Based on the common use characteristics of the codons identified in the above analysis, the expression efficiency of some exogenous genes in the Theaceae species can be improved, and the functional genes regulating the growth and development of the Theaceae species can be targeted (Peng et al., 2021).
There are many species of plants in Theaceae, and the interspecific relationship of this economically important plant is still a very controversial issue (Vijayan et al., 2009;Yang et al., 2013). The two classification approaches proposed by Chang (1998) and Ming (2000) diverge in many respects, especially with regard to the division of subgenera, sections and species. So far, it is uncertain which system has most accurately described the phylogenetic relationships within the Theaceae family. Therefore, it is necessary to look for other evidence to reconstruct the establishment of the Theaceae classification system. The squared Euclidean distances between plants within the genera Camellia and Stewartia were smaller (≤ 0.1), indicating that the CUB of the chloroplast genes were very similar (Wang, 2007). Based on the cluster analysis map constructed by RSCU values, most species of the same genus are clustered together, indicating that closely related species tend to have the same codon bias pattern. This is consistent with the conclusion in Liu et al. (2005) research that the closer the genetic relationship between species, the more similar the CU patterns. In addition, the clustering results based on RSCU values were similar to the topological structure of the phylogenetic tree constructed based on chloroplast genome sequences. In the topology of the phylogenetic tree, Camellia chekiangoleosa, Camellia crapnelliana, Camellia cuspidat, Camellia granthamiana, Camellia gymnogyna, Camellia impressinervis, Camellia pitardii, Camellia ptilophylla, Camellia renshanxiangiae, Camellia yunnanensis belongs to the first clade, which has the same opinion with . And Camellia pitardii was clustered with Camellia crapnelliana, Camellia cuspidate, Camellia impressinervis, and Camellia yunnanensis, indicating that Camellia pitardii might belong to subgenus Camellia, and it means they're very closely related, which is similar to . In addition, there are some inconsistencies between the phylogenetic tree of this study and others' Theaceae phylogeny.  is mainly focus on the phylogeny of Camellia, Camellia renshanxiangiae and Camellia sasanqua belong to the genus Camellia. However, in cluster analysis based on RSCU values, the two species of Camellia grouped into the Stewartia-Hartia clade. This provides a clearer direction for our subsequent research in the further. Therefore, the clustering method based on RSCU values can be used as a supplement to the molecular taxonomy and phylogenetic studies of Theaceae species. Furthermore, based on the genetic relationships revealed by cluster analysis, it is expected that in the future, transgenic technology can be used to conduct molecular experiments on genes underlying important traits in Theaceae species in order to facilitate in-depth studies on the cultivation of Theaceae species.

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

AUTHOR CONTRIBUTIONS
ZZ and ZJW designed the research. ZJW, QC, YW, ML, and CW analyzed the data and performed the experiments. ZJW, QC, YW, and ML wrote the paper. ZXW, CJ, CX, and HW conducted manuscript editing. ZZ revised the manuscript. All authors contributed to the article and approved to the published version of the article.


Rediscovery of Serangium montazerii Fürsch in Georgia and updated list of the Coccinellidae of Georgia

This article provides an updated list of the Coccinellidae of Georgia, recording 84 species. Serangium montazerii Fürsch nec S. parcesetosum Sicard is also rehabilitated as a species actually present in the country.

Here we propose an updated list of the Coccinellidae found in Georgia, taking into account the actual phylogeny and nomenclature. We also report the rediscovery of Serangium montazerii Fürsch, 1995nec Serangium parcesetosum Sicard, 1929 in an orange grove near Batumi, and make address nomenclatural confusion about the species in Georgia.

Site description
The sampled site is located at Chakvis Tskali, Kobuleti District, Adjara (41.71843°N, 41.73841°E) in a familial orange grove and garden. Further orange trees located 6 km eastward along the Chakvistskali River were also sampled. Observations were conducted on 25-VI-2019.
The climate of the region is characteristic of the Batumi area and the Black Sea Georgian coast, with high precipitations (avg. 2515 mm/year), mild winters (avg. 7 °C in January) and mild summers (avg. 23 °C in July). The region has warm temperate, almost subtropical, climate conditions (Bohn et al. 2000(Bohn et al. -2003. The vegetation is represented by hygrophilous thermophytic mixed deciduous broadleaved forest (Bohn et al. 2000(Bohn et al. -2003. The presence of many citrus orchards, with their specific pests, have made Batumi a focus of biocontrol attempts, with release of exotic Coccinellidae during the Soviet period (Timofeyeva and Hoang Duc Nhuan 1978).

Sampling method
Beating trays were used to collect ladybird beetles by beating and shaking the branches of orange trees. The ladybird beetles were then collected with an aspirator and subsequently identified to species level in the laboratory using available literature (Gourreau 1974, Iablokoff-Knzorian 1978, Poorani 1998, Bieńkowski 2018, 2020.

Serangium montazerii on Citrus sp.
This non-native species intentionally introduced in Georgia (Adjara) was initially identified as Serangium parcesetosum Sicard, 1929 (Timofeyeva andHoang Duc Nhuan, 1978). These authors also reassigned it to its original genus after the transfer to the genus Catana by Chapin (1940). They were followed by latter authors and by Escalona and Slipinski (2012).
However, Fürsch (1995) described a new species, Serangium montazerii Fürsch, 1995 from Iranian specimens (Mazandran Province). Subsequently Duverger (1998) as-signed all the introduced populations of Serangium in Western Europe, which were resulting from the first one introduced by Timofeyeva and Hoang Duc Nhuan (1978), to this latter species. However, after examination of specimens from Corsica, Coutanceau (2006), and Coutanceau and Malausa (2014) assigned the French introduced populations to S. parcesetosum. Their works, published in French in an amateurs' journal, remained largely unknown. Escalona and Slipinski (2012) in their revision of the Microweiseinae subfamily did not examine any members of the two taxa. The specimens introduced in France (continental and Corsica, Malausa et al. (1988)), and in Turkey (Yigit and Canhilal 2005)) came from this Georgian source (Timofeyeva and Hoang Duc Nhuan 1978). Likewise, neither Fürsch (1995), nor Duverger (1998), nor Coutanceau (2006) examined or compared the two taxa. However, Poorani (1998) noted that the two species can be distinguished by the shape of the right paramere and by the inter-ocular distance. Bieńkowski (2020) also underlined the same diagnostic characters and revised illustrations of the genitalia of both taxa. Moreover, the population introduced in Georgia in the 1970's (Timofeyeva and Hoang Duc Nhuan 1978) came from North India (Raniket, Uttarakhand) where S. montazerii occurs, while S. parcesetosum is distributed into Central and Southern India (Poorani 1998).
We found one specimen of S. montazerii on orange trees ( Figure 1). It was female and we could not perform a determination using genitalia. However, the shape of the head with an inter-ocular distance equal to twice the size of the eyes ( Figure 1B) allows us to assign this specimen to S. montazerii. The location is close to Batumi, and corresponds to the area of first historical introduction of the species in Georgia, the country of first use of the species for biological control of D. citri (Timofeyeva and Hoang Duc Nhuan 1978). According to Yigit and Canhilal (2005), the species was still present in Georgia in 1990. It has probably been present in the country since 1975, but was not listed by Merkviladze and Kvavadze (2002). The species is also present northward in the vicinity of Sochi (Orlova-Bienkowskaja and Bieńkowski 2017; Bieńkowski 2018, 2020; Bieńkowski and Orlova-Bienkowskaja 2020) in Russia.

Tribe Diomini Gordon, 1999
• One species is recorded from Georgia.

Conclusion
The Georgian fauna of Coccinellidae is containing 84 species. Two species (S. biflammulatus and S. montazerii) have been added to the Merkviladze and Kvavadze (2002) list, but one species is dubious for the country (T. conatus); a further two species are now considered as synonyms. Thirty species were not listed in the Catalogue of Palearctic Coleoptera (Kovář 2007).  (Kovář 2007).