Predicting Heart Cell Types by Using Transcriptome Profiles and a Machine Learning Method

The heart is an essential organ in the human body. It contains various types of cells, such as cardiomyocytes, mesothelial cells, endothelial cells, and fibroblasts. The interactions between these cells determine the vital functions of the heart. Therefore, identifying the different cell types and revealing the expression rules in these cell types are crucial. In this study, multiple machine learning methods were used to analyze the heart single-cell profiles with 11 different heart cell types. The single-cell profiles were first analyzed via light gradient boosting machine method to evaluate the importance of gene features on the profiling dataset, and a ranking feature list was produced. This feature list was then brought into the incremental feature selection method to identify the best features and build the optimal classifiers. The results suggested that the best decision tree (DT) and random forest classification models achieved the highest weighted F1 scores of 0.957 and 0.981, respectively. The selected features, such as NPPA, LAMA2, DLC1, and the classification rules extracted from the optimal DT classifier played a crucial role in cardiac structure and function in recent research and enrichment analysis. In particular, some lncRNAs (LINC02019, NEAT1) were found to be quite important for the recognition of different cardiac cell types. In summary, these findings provide a solid academic foundation for the development of molecular diagnostics and biomarker discovery for cardiac diseases.

Introduction
The heart is a complex organ containing various cardiac cell types, and the interaction between different heart cell types could realize the important functions of the heart. Previous pioneering studies have shown that the heart is composed of approximately 70% non-cardiomyocytes and 30% cardiomyocytes [1]. Cardiomyocytes could be divided into atrial myocytes and ventricular myocytes, while non-cardiomyocytes mainly include fibroblasts, smooth muscle cells, pericytes, and endothelial cells. These cells form four chambers with different morphologies and functions, and they complete the systemic blood circulation [2].
Cardiomyocytes are responsible for contractile function, and they are the most important part. However, they do not function in isolation. Fibroblasts account for more than 40% of the total cells in the ventricle. Their core function is to maintain the cardiac extracellular matrix homeostasis and provide structural and mechanical support for the cardiomyocytes [3]. The mural cells of the vessel wall are mainly composed of smooth muscle cells and pericytes, and these two cell types are important for vascular integrity and heart function [4]. Endothelial cells form the inner layer of blood and lymphatic vessels; they maintain blood circulation by regulating the permeability and caliber of blood vessels and play an important role in controlling and maintaining the growth, contractility, and rhythm of the heart [5,6]. Mesothelial cells are transitional mesodermal-derived cells with similar morphological and functional characteristics to endothelial cells. They can secrete angiogenic factors, which are important for angiogenesis [7]. Heart adipose tissue not only can supply energy locally but also has heart repair functions, such as new blood vessel formation and immune regulation [8,9]. Immune cells and neurons are also very important for the functional homeostasis [10,11].
Identifying cell components and cell types are important for understanding cell functions, especially in complex organs, where multiple cell types work together. There are two types of traditional methods for cell type annotation: (1) cell marker-based methods, such as CellAssign [12], which needs high quality cell specific expressed genes, but most cell types do not have very specific biomarkers. (2) reference dataset-based methods, such as SingleR [13], which compare the scRNA data with reference scRNA data with known cell types and make predictions. Previous studies have reported some markers for cardiac cells, such as in atria (NPPA and SLN), ventricles (MYL2 and MYL3), endothelial cells (FABP4 and AQP7), smooth muscle cells (ACTA2), fibroblasts (COL1A1), pericytes (PDGFRB) or immune cells (PTPRC), neurons (NEXN1), and adipocytes (GPAM and FASN) [14][15][16][17]. Although these genes are very important for each cell type, the maintenance of cell function depends on the interaction among different genes. Therefore, revealing the specific expression patterns of different cell types, especially the expression features that distinguish them from other cell types, is very important for an enhanced understanding of fate decisions and cell functions.
On the basis of existing single cell profiling datasets from the Human Cell Atlas study of adult human heart cells [17], machine learning algorithms were used in the present study to analyze data to extract gene expression characteristics and biomarkers to characterize different heart cell types. Machine learning algorithms can extract hidden biomarkers that cannot be found by traditional methods through in-depth research on a large number of cell gene expression data. By using the light gradient boosting machine (lightGBM) algorithm [18], a ranking feature list was generated on the basis of the importance of these features. Then, the incremental feature selection (IFS) method [19] with decision tree (DT) [20] and random forest (RF) [21] algorithms was applied to determine the best number of features and build the optimal classifier. As a result, the most relevant gene features and decision rules were identified. Through these rules, 11 cell types could be accurately classified. Meanwhile, the results of Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis suggested that the selected genes may have important significance for the phenotype of cells or function of the heart. Further research on these genes may help clarify the detailed mechanism of heart development. In short, this research identified a group of potential cardiac cell biomarkers and a precise classifier composed of many decision rules, thus providing insights for further research on cardiac development and function. The method we proposed to identify heart cell types was a type of reference dataset-based method. But unlike traditional reference dataset-based methods which use all the genes, we only use discriminative genes selected with feature selection Life 2022, 12, 228 3 of 17 methods. The method had the merits of both two types of traditional methods: (1) we only considered the discriminative genes; therefore, it was faster and more explainable than the refence dataset-based methods; (2) we quantitively considered the combinations of expression levels, i.e., expression rules, rather than only cell specific markers; therefore, it was more accurate than the cell marker-based methods. The following text was organized as follows: (1) Section 2 lists the dataset analyzed in this study and algorithms used in this study; (2) The results were presented in Section 3; (3) Extensive discussion on results was provided in Section 4; (4) Section 5 summarized this study.

Study Design
In this study, the interpretation of the model (i.e., classifier) was consisted of two parts, that is, interpretation on (1) single gene expression signatures and (2) combined gene expression rules. Single-gene interpretation explained the effect of single optimal genes on the classification, while combined-gene interpretation focused on exploring how multiple genes contribute to the classification together. The whole framework of this study is shown in Figure 1. dataset-based methods which use all the genes, we only use discriminative genes selected with feature selection methods. The method had the merits of both two types of traditiona methods: (1) we only considered the discriminative genes; therefore, it was faster and more explainable than the refence dataset-based methods; (2) we quantitively considered the combinations of expression levels, i.e., expression rules, rather than only cell specific markers; therefore, it was more accurate than the cell marker-based methods. The follow ing text was organized as follows: (1) Section 2 lists the dataset analyzed in this study and algorithms used in this study; (2) The results were presented in Section 3; (3) Extensive discussion on results was provided in Section 4; (4) Section 5 summarized this study.

Study Design
In this study, the interpretation of the model (i.e., classifier) was consisted of two parts, that is, interpretation on (1) single gene expression signatures and (2) combined gene expression rules. Single-gene interpretation explained the effect of single optima genes on the classification, while combined-gene interpretation focused on exploring how multiple genes contribute to the classification together. The whole framework of thi study is shown in Figure 1.
. Figure 1. Flow chart of the study design. First, lightGBM method was applied to rank the feature of single-cell gene expression profiles into a ranked list. Second, IFS method with machine learning Figure 1. Flow chart of the study design. First, lightGBM method was applied to rank the features of single-cell gene expression profiles into a ranked list. Second, IFS method with machine learning algorithms was used to detect the best number of features and build the optimal classifiers and decision rules. Finally, functional enrichment analysis was performed on the optimal gene feature set.

Data Collection
Raw dataset was downloaded from the publicly available Human Cell Atlas Data Coordination Platform, with accession number: ERP123138 (https://www.ebi.ac.uk/ena/ browser/view/ERP123138, accessed on 29 January 2021) [17]. The processed 10X Genomics dataset included the expression profiles of 33,538 genes in 451,513 cells from 11 cell types. The 11 major cell types and the corresponding sample sizes are presented in Table 1.

Feature Ranked by LightGBM
LightGBM is a well-known boosting learning machine [18] that combines many weak classifiers to achieve a single strong one. It could be regarded as an improved version of gradient boosting DT (GBDT) [22], which recurrently fits a new DT by using the negative gradient of the loss function of the current DT as the approximate value of the residual. The main differences between lightGBM and gradient-based one-side sampling (GOSS) lie in the two new strategies adopted by lightGBM. These are GOSS and exclusive feature bundling (EFB), which both greatly improve the efficiency and ensure the accuracy of classification. In GBDT, the gradient of a sample in calculating the residual error of a DT reflects the contribution of the sample to subsequent classification. Therefore, GOSS down samples the training data by randomly screening out most of the samples with small gradient and keeping a small number of them to maintain the distribution of the data. EFB bundles the mutually exclusive features together to reduce the dimension of the data. The mutually exclusive features are those that rarely take nonzero values simultaneously, and no or very little information is lost by bundling them as a new feature. EFB is realized by solving a graph coloring problem with the use of a greedy algorithm with a constant approximation ratio. As described in LightGBM's documentation (https://lightgbm.readthedocs.io/en/ latest/, accessed on 10 May 2020), the advantages of lightGBM include faster training efficiency, low memory usage, higher accuracy, support for parallel learning, and being able to handle large-scale data. In addition to classification, lightGBM sorts features in accordance with their importance, which is quantified by the number of times the feature is selected to build DTs. The more times a feature is used, the higher the ranking. In this research, the features were sorted using lightGBM for further analysis. The lightGBM program was implemented through a Python module.

IFS Method
Once the feature list F was generated by sorting the gene features with the lightGBM method, the number of significant features could still not be determined. Here, the optimal number of features was discovered using the IFS method [19]. IFS first generates a series of feature subsets from the feature list F on the basis of the specific step size. For example, when the step size is 5, the first feature subset f 1 generated is the top five features in F, the second feature subset f 2 is the top 10 features in F, and so on. Next, each classifier trains on a training set whose samples are expressed by the features in a feature subset. The performance of the classifier is evaluated using 10-fold cross-validation [23] and synthetic minority oversampling technique (SMOTE) [24], and the classifier with the best performance is considered the optimal classifier, and its trained subset is regarded as the optimal subset.

Classifier Building with DT and RF
In this study, DT and RF were adopted to build classifiers. Their descriptions are as follows: RF [21,[25][26][27][28] is a classification algorithm that integrates multiple DT classifiers. Through a bootstrap resampling technique, a new training set is composed, and the DT is constructed by randomly selecting samples and features from the original data set. The prediction labels of RF are obtained by combining the prediction results of multiple DT classifiers by using the principle of minority rule. RF has few parameters to tune, users can only choose a proper number of DTs so that it can yield good performance. Because RF contains several DTs, it has excellent noise tolerance and can avoid the overfitting problem. Importantly, although DT is a relative weak classifier, RF is much stronger. Thus, it is widely used in omics research.
Through RF, investigators can build an efficient classifier. However, it is a black-box classifier. Its classification principle is not easy to understand. Accordingly, we cannot extract essential difference of heart cells in different types. In view of this, DT [20,29,30] was also employed in this study, which is widely used in the field of biomedical research, because the decision rules generated by DT can effectively elucidate how decisions are made for classification or regression tasks. In another word, a DT is a white-box classifier that splits the data many times on the basis of certain thresholds in the features by using the IF-THEN format. These IF-THEN formats constitute decision rules, which can clearly exhibit special patterns for each class. In this study, these patterns indicated specific characteristic of heart cell types. Although DT provides relative low performance, it can give novel insights to study heart cells. In the output file of DT, the value of "passed counts" indicates the number of samples satisfying the condition of the rule. The above two algorithms were performed by the scikit-learn program with default parameters in Python [31].

SMOTE
The different number of samples from different cell types leads to the problem of data imbalance. Synthetic minority oversampling technique (SMOTE) was applied to minimize the effect of sample imbalance on the construction of classifiers [24]. It generates many synthetic samples for minority cell categories on the basis of the principle of k-nearest neighbors [32]. For each cell type, except for the cell type with the highest number of samples, new synthetic samples were added to this cell type via SMOTE until the number of samples of each type was almost the same. The SMOTE program was accessed from https://github.com/scikit-learn-contrib/imbalanced-learn (accessed on 24 March 2021), and the parameters were set to default.

Performance Measurement
As the number of samples in each category in the dataset could be strongly unbalanced, the weighted F 1 score [33][34][35] is an appropriate measurement of the classifier's performance. First, the F 1 score was calculated using the following formula: Next, weighted average was performed to the F 1 score of each category, and the weight was the proportion of the number of categories in the correct label. This measurement was called weighted F 1 score.

Functional Enrichment
For functional enrichment analysis, the Clusterprofile package was applied to the annotation of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, where p < 0.05 was taken as the screening criterion [36]. The GO terms are divided into three subgroups, namely, biological process (BP), cellular component (CC), and molecular function (MF).

Results of LightGBM Method on the Dataset
In this study, single cell expression profiles of 451,513 cell samples and 11 cell types for heart disease were obtained, and each sample of heart cell type was represented by the expression of 33,537 genes. LightGBM method was first applied to rank the genes into a feature list on the basis of feature importance to filter out the important set of features from these genes, and the results are provided in Table S1. The top 20 genes in the list are shown in Table 2.

Results of IFS Method with RF
A feature ranking list was obtained by the lightGBM method, but the optimal number of features was not determined. The IFS method was applied to optimize the selected gene features. We first adopted RF to execute IFS method. The propuse was to construct an efficient classifier for classifying heart cells.
Based on the feature list provided in Table S1, the IFS method produced a series of feature subsets with the step size of 5. A SVM classifier was then built based on each feature subset to predict the label of each sample. To save time, we only considered top 1000 gene features in the list. The evaluation results of all SVM classifiers are shown in Table S2. The IFS curve with the number of features as the X-axis and the performance of each classifier, measured by weighted F1score, as the Y-axis is drawn in Figure 2. RF obtained the optimal weighted F1 score of 0.981 with the top 470 features. Accordinly, the optimum RF classifier can be built using these 470 features. Table 3 provides the detailed evaluation metrics of the optimal RF classifier. It also provided good performance on other metrics. In addition, the F1 score for each cell type under the optimal RF classifier is presented in Figure 3. Such optimum RF classifier showed excellent performance on the prediction of each category, indicating the effectiveness of this classifier. RF classifier can be built using these 470 features. Table 3 provides the detailed evaluation metrics of the optimal RF classifier. It also provided good performance on other metrics. In addition, the F1 score for each cell type under the optimal RF classifier is presented in Figure 3. Such optimum RF classifier showed excellent performance on the prediction of each category, indicating the effectiveness of this classifier.       Figure 3. Performance of two optimal classifiers on each cell type. Figure 3. Performance of two optimal classifiers on each cell type.
As mentioned above, the optimal RF classifier provided quite good performance. It can be an efficient tool to classify heart cells. To further elaborate its robustness, the following test was conducted. First, some noise was added to the dataset, in which cells were represented by features used in the optimal RF classifier. In detail, we randomly selected 10% cells. Each feature of these cells randomly increased or decreased by a small number. On such dataset, the optimal RF classifier was evaluated by 10-fold cross-validation. Above procedures were executed for ten times, producing ten weighted F 1 scores. A box plot was drawn (Figure 4) to show these values. It can be observed that the performance on datasets with noise was almost same as that on the original dataset. This proved the robustness of the optimal RF classifier. plot was drawn (Figure 4) to show these values. It can on datasets with noise was almost same as that on the robustness of the optimal RF classifier.

Results of IFS method with DT
Although the optimal RF classifier exhibited quite box classifier that fails to explain the decisions. To extr heart cells, DT, a white-box algorithm, was employed. W on DT that had been done for RF. The performance of subsets is listed in Table S2. An IFS curve was also plott observed that the highest weighted F1score was 0.957 Thus, an optimum DT classifier can be built with thes classifier are provided in Table 3. The performance of shown in Figure 3. Evidently, the performance of the than that of the optimum RF classifier, which conform powerful than DT. However, DT has its own merit that procedures of DT were completely open, which make u sification principle, thereby giving new insights to unco different types.

Classification Rules Generated by the Optimal DT Class
The optimum DT classifier was built based on top 3 . Box plot to show the performance of the optimal RF classifier on datasets with noise. The performance is almost same as that on the original dataset, proving the robustness of the optimal RF classifier.

Results of IFS method with DT
Although the optimal RF classifier exhibited quite good performance, it is a black-box classifier that fails to explain the decisions. To extract more insights for the study of heart cells, DT, a white-box algorithm, was employed. We conducted the same procedures on DT that had been done for RF. The performance of DT classifiers on different feature subsets is listed in Table S2. An IFS curve was also plotted, as shown in Figure 2. It can be observed that the highest weighted F1score was 0.957 when top 380 features were used. Thus, an optimum DT classifier can be built with these features. Other metrics of such classifier are provided in Table 3. The performance of such classifier on all cell types is shown in Figure 3. Evidently, the performance of the optimum DT classifier was lower than that of the optimum RF classifier, which conform to the general fact that RF is more powerful than DT. However, DT has its own merit that RF cannot have. The classification procedures of DT were completely open, which make us possible to understand its classification principle, thereby giving new insights to uncover the difference of heart cells in different types.

Classification Rules Generated by the Optimal DT Classifier
The optimum DT classifier was built based on top 380 features. Accordingly, we used these features to represent all heart cells. Such representation of all heart cells was learnt by DT. A large tree was obtained, from which we constructed 11,139 interpretable rules. The detailed rules are listed in Table S3. The number of rules for each category is shown in Figure 5. Endothelial cells obtained the largest number of rules, with 2588, followed by atrial cardiomyocytes and pericytes. A detailed description of these rules could be found in the Section 4.2.

Functional Enrichment Analysis with Optimal Gene Set
The best gene set was obtained, including the top 470 features, by using the IFS method. These genes were analyzed by GO and KEGG pathway functional enrichment, and the results are presented in Table S4 and Figure 6. Many genes were enriched in the KEGG pathway of hypertrophic cardiomyopathy and dilated cardiomyopathy, and some genes were enriched in the BP of heart process, indicating that these genes may be associated with the development of heart disease and further demonstrating the effectiveness of the method.

Functional Enrichment Analysis with Optimal Gene Set
The best gene set was obtained, including the top 470 features, by using the IFS method. These genes were analyzed by GO and KEGG pathway functional enrichment, and the results are presented in Table S4 and Figure 6. Many genes were enriched in the KEGG pathway of hypertrophic cardiomyopathy and dilated cardiomyopathy, and some genes were enriched in the BP of heart process, indicating that these genes may be associated with the development of heart disease and further demonstrating the effectiveness of the method.

Discussion
In this study, several machine learning algorithms were applied to classify different cell types from single-cell and single-nucleus transcriptome profiles of six adult heart regions. First, the lightGBM method was performed to obtain a ranking feature list in accordance with feature importance. Second, an RF algorithm was applied to construct a

Discussion
In this study, several machine learning algorithms were applied to classify different cell types from single-cell and single-nucleus transcriptome profiles of six adult heart regions. First, the lightGBM method was performed to obtain a ranking feature list in accordance with feature importance. Second, an RF algorithm was applied to construct a precise classifier with a high classification accuracy of 0.981. However, as a black-box classifier, this classifier cannot reveal the different expression patterns of different heart cell types. Therefore, the DT algorithm was further used to obtain a group of decision rules. By using the top 380 features, different cell types could be distinguished, with a high classification accuracy of 0.957. Furthermore, novel biomarkers or expression patterns may be identified by analyzing the expression pattern of 310 selected genes within these decision rules.

Candidate Gene Expression Features Discriminating Different Heart Cells
A total of 470 selected features (genes) were used to classify 452,136 heart cells into 11 cell types. Among them, the top ranked genes are usually more decisive for distinguishing different cell types. Some relevant experimental evidence that supported the results were presented. NPPA (ENSG00000175206) encodes natriuretic peptide A (ANP), which is highly expressed in the heart muscle and related to the control of extracellular fluid volume and electrolyte homeostasis. Studies have found that NPPA is expressed primarily in the heart, and atria have higher expression than ventricles. NPPA can regulate vasodilation, natriuresis, diuresis, and aldosterone synthesis and further influence blood pressure. Moreover, it is involved in inhibiting cardiac hypertrophy, cardiac fibrosis, and cardiac remodeling by inducing cardiomyocyte apoptosis and attenuating the growth of cardiomyocytes and fibroblasts [37]. In adipocytes, ANP can promote white adipocyte browning to increase energy expenditure via a PKG-p38 mitogen-activated protein kinase mediated pathway and make the heart as a central regulator of adipose tissue biology [38]. These findings are consistent with the results of the present study, which showed that NPPA has a crucial role in different functional heart cells.
Similarly, gene LAMA2 (ENSG00000196569), which encodes the laminin subunit, plays an important role in normal heart function. Studies have demonstrated that homozygous mutation of LAMA2 can cause unstable myotube formation in various cardiac muscle, and abnormal LAMA2 expression may lead to heart dieases, such as cardiomyopathy, heart failure, and dilated cardiomyopathy [39,40]. A previous adult human heart research also showed that LAMA2 has different expression levels in fibroblasts, cardiac adipocytes, and other cell types [17].
Other key features in our results are also important for cardiac function. For example, the gene DLC1 (ENSG00000164741) is highly expressed in endothelial cells and a small number of ventricular cardiomyocytes [17,41]; RYR2 (ENSG00000198626) encodes a calcium channel component associated with cardiomyocyte and smooth muscle cell contraction and thermogenesis in beige adipocytes [42][43][44]; TTN (ENSG00000155657) encodes titin, a large abundant protein of striated muscle, mainly found in cardiac and skeletal human muscle. Mutations in these genes may cause a variety of cardiac diseases [45][46][47].
More importantly, we found that some long non-coding RNAs (lncRNA) are important for differentiating cardiac cell types. For example, LINC02019 (ENSG00000273356), it is the top gene in our feature ranking. The product of this gene is a long intergenic nonprotein coding RNA (lincRNA). We used the starBase tool to study its related RNA-binding proteins [48]. The most related proteins include EIF4A3, ELAVL1, LIN28A. Studies have shown that EIF4A3 is associated with acute myocardial infarction, and knockout of EIF4A3 can lead to failure of heart looping [49,50]. ELAVL1 plays an important role in inhibiting hyperglycemia-induced cardiomyocyte pyroptosis and regulating ferroptosis in myocardial injury [51,52]. LIN28A is also implicated in various cardiac injuries or diseases [53,54]. NEAT1 (ENSG00000245532) produces a lncRNA that may act as a transcriptional regulator for numerous genes. NEAT1 was markedly downregulated in cardiomyocytes following ischemia-reperfusion-induced injury. Moreover, by interacting with microRNA-125a-5p, NEAT1 could modulate the concentration of BCL2L12, which in turn regulates cardiomyocyte apoptosis [55]. Other studies also found that NEAT1 may influence myocardial injury repair through the MAPK and TLR2/NF-κB signaling pathways [56,57].
In summary, these genes are all related to cardiac structure and function, and they show various expression levels in the ventricle, atrium, and other cell types. Therefore, these genes could be used as decisive feature for distinguishing different cardiac cells, and we also confirmed that some lncRNAs may have more specific roles in the maintenance of normal cardiac function.

Candidate Gene Expression Rules Discriminating Different Heart Cells
Through DT method, a classifier consisting of 11,139 decision rules involving 380 selected features was built. Each cell type was assigned some rules, as shown in Figure 4. According to the value of "passed counts", top three rules for each cell type were extracted and listed in Table 4. The genes involved in these 33 rules were analyzed in combination with the existing literature to prove the reliability of the results. Studying other rules may also help find some new characteristics of cardiac cell subtypes, which may provide new insights into the in-depth understanding of cardiac development and function. Endothelial 8820 a : "passed counts" indicates the number of samples satisfying the condition of the rule.

Cardiomyocytes
Cardiomyocytes generate contractile force; thus, they normally show high expression of sarcomere proteins and calcium-mediated processes [17]. The present study showed that atrial cardiomyocytes and ventricular cardiomyocytes highly expressed TTN, which was the most important factor for distinguishing cardiomyocytes and non-cardiomyocytes. As mentioned above, atria showed higher expression of NPPA than ventricles [58]. The decision rules in the present study also showed the same expression pattern. In addition, atrial cardiomyocytes required a higher expression of KCNJ3 (ENSG00000162989) and MYL7 (ENSG00000106631) than ventricular cardiomyocytes. KCNJ3 encodes an integral membrane protein and an inward-rectifier type potassium channel. Studies have found that KCNJ3 plays an important role in governing cardiac electrical activity, and atrial cardiomyocytes have high KCNJ3 levels than ventricular cardiomyocytes [59]. Other studies also demonstrated that the protein encoded by MYL7, which is a part of myosin, is highly expressed in atrial cardiomyocytes [60].

Fibroblasts and Vascular, Stromal, and Mesothelial Cells
In the decision rules, smooth muscle cells showed higher expression of MYH11 (ENSG00000133392) than other cell types. MYH11 encodes myosin heavy chain 11, and high MYH11 level is a marker of mature contractile phenotype of smooth muscle cell, while downregulation or mutation of MYH11 is associated with vascular disease [61]. ACTA2 (ENSG00000107796) is also a marker of smooth muscle cells; it is known as smooth muscle α actin, which is usually highly expressed in smooth muscle cells, pericytes, and myofibroblasts [62]. Studies found that mutation of ACTA2 could cause coronary artery disease and thoracic aortic disease. An experimental study on smooth muscle cells and myofibroblasts harboring ACTA2 mutations indicated that occlusive disease is associated with increased proliferation of smooth muscle cells [63,64]. ABCC9 (ENSG00000069431) encodes a protein that is a member of the ATP-binding cassette transporter superfamily. In this research, the decision rules of pericytes showed high ABCC9 expression requirement. This finding is in accordance with various published studies, which showed that ABCC9 is highly expressed in pericytes and could be a biomarker for it [17,65,66].
VWF (ENSG00000110799) was required to be highly expressed in endothelial cells in the decision rules. As a protein coding gene, VWF encodes a glycoprotein involved in hemostasis, and it is reported to be highly expressed in endothelial cells in human cells [67], thus confirming the results of the present study.
Unlike endothelial cells, smooth muscle cells, pericytes, and mesothelial cells, fibroblasts are not involved in the construction of basement membranes. As known, collagen IV, laminin-entactin/nidogen complexes, and proteoglycans are the major molecular constituents of basement membranes [68]. This information was confirmed in the decision rules because LAMA2 (ENSG00000196569) and CD36 (ENSG00000135218) showed a relatively low expression in fibroblasts. As CD36 encodes collagen IV and LAMA2 encodes a subunit of laminin, they both required to be highly expressed in other four cell types.
Although existing studies could not confirm why mesothelial cells need to relatively highly express PLA2G2A in the rules, this may be a coincidence caused by the limited number of mesothelial cells or the unknown effect of PLA2G2A on mesothelial cells. Similarly, other types of cells have some very meaningful genes that have not been studied in depth. However, the results of this study showed that they may be important.

Adipocytes and Immune and Neuronal Cells
Adipocytes, immune cells, and neuronal cells showed low expression of sarcomere proteins or basement membranes components. NRXN1 (ENSG00000179915) encodes neurexin 1, which is a cell surface receptor involved in the formation of synaptic contacts, and efficient neurotransmission depend on NRXN1 [69]. The same expression pattern could be observed from the decision rules of neuronal cells. Another highly expressed gene in the rules was neuronal growth regulator 1 (NEGR1, ENSG00000172260). NEGR1 mediates neural cell communication and synapse formation, and its downregulation is related to obesity, learning difficulties, intellectual disability, and psychiatric disorders [70]. Thus, appropriate NEGR1 expression is necessary to maintain neuronal cell function.
Hematopoietic cells are commonly classified into myeloid and lymphoid cells. The expression of CD163 (ENSG00000177575) is the main criterion used to distinguish between myeloid and lymphoid cells in the rules. As CD163 encodes a member of the scavenger receptor cysteine-rich superfamily, it is exclusively expressed in monocytes and macrophages [71]. Thus, CD163 could serve as a marker for myeloid cells, and this finding also confirmed the results of the present study. Adipocytes showed high PLIN1 (ENSG00000166819) and ACACB (ENSG00000076555) levels than other cells. PLIN1 and ACACB participate in the inhibition of lipolysis and the regulation of fatty acid oxidation, respectively. According to the recent publications, these two genes are essential for the maintenance of adipocyte functions [72,73].
In conclusion, the genes that are highly expressed in the decision rules are often markers of the cell types or essential for maintaining cell function, which reflects the reliability of the research results. Some genes that are meaningful for cell classification, which have not been investigated in detail before, may have important implications for the function of the corresponding cell type.

Functional Analysis of the Optimal Gene Set
We also performed GO and KEGG pathway enrichments in the 470 decisive features identified by the RF method, and here we presented the relevant enrichment results of some cell types.
The enrichment terms related to adipocytes include GO:0019216 (regulation of lipid metabolic process), GO:0036041 (long-chain fatty acid binding) and hsa04923 (Regulation of lipolysis in adipocytes), which reflect the energy supply and regulation of adipocytes. Similarly, there are also some terms related to lymphoid cells and neurons in the enrichment results, such as GO:0140058 (neuron projection arborization) and GO:0035325 (Toll-like receptor binding).
The primary function of the heart is to effectively pump blood to the body tissues through the contraction-relaxation cycle of myocytes, and the heart is mainly composed of cardiomyocytes and fibroblasts. In our GO and KEGG enrichement results, many of the results are related to these two types of cells. Fibroblast-related terms are mainly related to fibers, such as GO:0043292 (contractile fiber), GO:0030016 (myofibril), and GO:0030017 (sarcomere). Sarcomere is composed of a segment of myofibril between two adjacent Z discs, and it is the contractile unit of myofibrils. Studies have found that various genetic mutations in the cardiac sarcomere could lead to defects in sarcomere production and further lead to ventricular dilatation and cardiac dysfunction [74]. The enrichment terms of cardiomyocytes mainly include GO: 0048738 (cardiac muscle tissue development), GO: 0003779 (actin binding) and GO: 0008307 (structural constituent of muscle). Actin plays an essential role in the assembly of cardiac myofibrils, and it is strongly associated with muscle contraction. In the process of myofibrillogenesis, actin is assembled into highly ordered mature state, and the abnormal expression of actin dynamic-related genes usually leads to myofibril abnormalities and heart defects [75]. Meanwhile, the KEGG pathway annotation showed that the features are related to muscle contraction and various cardiomyopathies (KEGG: hsa05410, hypertrophic cardiomyopathy; KEGG: hsa04260, cardiac muscle contraction). The main cause of cardiomyopathy is muscle cell dysfunction caused by genetic variation, especially the dysfunction of genes related to the cytoskeleton-sarcomere connection [76,77]. This further illustrates the important role of these genes in cardiac muscle cell function.
In this part, an extended description about the decisive features was discussed. These terms have been shown to be associated with multiple cardiac cell types, further demonstrating the importance of our decisive features to the function of cardiac cell populations and supporting the reliability and effectiveness of their usefulness in constructing classifiers.

Limitations of This Study
There are several limitations of this study: (1) The data we analyzed was from Human Cell Atlas and it was not known how representative the Human Cell Atlas data is. When more data from various heart locations becomes available, we need to test the model on more data. (2) The tissue was only from adult human hearts. Can the model be applied to children or infants? Do they have different cell types? We need to compare the adult heart data with child and infant data. were measured with smart-seq, will other low expressed biomarkers and rare cell types be discovered? There may be more cell markers and cell types than the identified ones.

Conclusions
The feature selection method and machine learning algorithms were applied to build a workflow to determine the essential gene features and specific gene expression rules for classifying different heart cell types. The positive results received from this study indicated that the developed classification models achieved an excellent classification performance. The selected key genes and decision rules were demonstrated to be associated with cardiac structure and function in recently published literature and enrichment analysis. In the future, we will collect data from different heart locations, different age groups and different sequencing platforms to get robust heart cell type annotation model. With the knowledge of gene expression patterns at single cell resolution, we can decipher the cell composition changes for diagnosis and target the dysfunctional cell without harming other cells for precision treatment.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/life12020228/s1, Table S1: Ranking results of features, as obtained using lightGBM method; Table S2: Evaluation of the performance of each classifier with different numbers of features by using IFS method with 10-fold cross validation; Table S3: Decision rules generated by DT classifier when using the top 380 features; Table S4

Data Availability Statement:
The data presented in this study are openly available in Human Cell Atlas Data Coordination Platform at https://www.ebi.ac.uk/ena/browser/view/ERP123138 (accessed on 29 January 2021), reference number [15].

Conflicts of Interest:
The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.


Characterization and subcellular localization of Alongshan virus proteins

Alongshan virus (ALSV) in the Jingmenvirus group within the family Flaviviridae is a newly discovered tick-borne virus associated with human disease, whose genome includes four segments and encodes four structural proteins (VP1a, VP1b, VP2, VP3, and VP4) and two non-structural proteins (NSP1 and NSP2). Here, we characterized the subcellular distribution and potential function of ALSV proteins in host cells. We found that viral proteins exhibited diverse subcellular distribution in multiple tissue-deriving cells and induced various morphological changes in the endoplasmic reticulum (ER), and NSP2, VP1b, VP2, and VP4 were all co-localized in the ER. The nuclear transfer and co-localization of VP4 and calnexin (a marker protein of ER), which were independent of their interaction, were unique to HepG2 cells. Expression of NSP1 could significantly reduce mitochondria quantity by inducing mitophagy. These findings would contribute to better understanding of the pathogenesis of emerging segmented flaviviruses.

Introduction
In recent years, various pathogenic tick-borne viruses have been identified to be associated with human febrile illness (Dong et al., 2021;Ma et al., 2021;Wang et al., 2021). Alongshan virus (ALSV), a segmented flavivirus in the Jingmenvirus group of the family Flaviviridae, has been detected in tick-bitted patients and Ixodes persulcatus in Northeast China . Subsequently, ALSV is also found in Ixodes persulcatus in Russia and in Ixodes ricinus in Finland (Kuivanen et al., 2019;Kholodilov et al., 2020), whose emergence may be associated with its codon usage (Rahman et al., 2022). The genome of ALSV is divided into four segments (S1-S4): S1 encodes flavivirus NS5-like non-structural protein 1 (NSP1), S2 encodes structural glycoproteins VP1a, VP1b and VP4, S3 encodes flavivirus NS2b3-like non-structural protein NSP2, and S4 encodes Frontiers in Microbiology 02 frontiersin.org nucleoprotein VP2 and membrane protein VP3 Zhang et al., 2020). The structural features of ALSV NSP2 is conserved in comparison with flavivirus NS3 helicases (Gao et al., 2020). Patients with ALSV infection mainly present with fever and headache, and some patients are accompanied by clinical symptoms of fatigue, coma, depression, nausea, myalgia/ arthralgia, and skin rash . There is no effective vaccine and drug available for ALSV. A better understanding of the pathogenic mechanism of ALSV would contribute to the development of vaccines and antiviral drugs. Flavivirus particles enter cells via endocytosis, exposing them to an acidic endosome that triggers the release of viral RNA into the cytosol (Modis et al., 2004;Zhang et al., 2015). The released positivesense RNA is recognized by ribosomes and translated at the rough endoplasmic reticulum (ER) into a single polyprotein, which is cleaved into structural and non-structural proteins by the viral and cellular proteases. The non-structural proteins orchestrate ER membrane rearrangement to form the single-membrane invaginations that house the replication machinery (Neufeldt et al., 2018). Then, the viral gRNA molecules generated by the viral replicase complex are incorporated into viral particles which involve RNA encapsidation and budding into the lumen of the ER (Welsch et al., 2009). Finally, the newly synthesized virus particles are transported from ER and Golgi apparatus to the cell surface, and released into extracellular by exocytosis (Xin et al., 2017). During the replication process, viruses manipulate the cellular organelles of the host to propagate and cause disease. Various flaviviruses, such as Zika virus (ZIKV), Dengue virus (DNEV), and hepatitis C virus (HCV), have evolved multiple strategies to drive cellular reprogramming by inducing a series of organelles (such as endoplasmic reticulum, Golgi apparatus, and mitochondria) structural and functional changes that promote viral replication (Neufeldt et al., 2018).
The subcellular trafficking with host cells and biological functions of ALSV proteins remain largely unknown. In the present study, we cloned all the genes of ALSV proteins into eukaryotic vectors for expression in various mammalian cells and examined their subcellular location using immunofluorescent assay. We found that in multiple tissue-deriving cells, viral proteins exhibited diverse subcellular distribution and induced various ER morphological changes, and NSP2, VP1b, VP2, and VP4 proteins were all co-localized with ER. The nuclear transfer and co-localization of VP4 and calnexin, which independent of their interaction was unique to HepG2 cells, and NSP1 reduced mitochondria quantity by inducing mitophagy. These findings would contribute to the understanding of the pathogenesis of emerging segmented flaviviruses.

Plasmids construction and transfection
Flag-tagged ALSV NSP1, NSP2, VP1a, VP1b, VP2, VP3, and VP4 encoding genes were cloned into VR1012-based expression vector and confirmed by sequencing. Transfection of the plasmids into the indicated cells using Lipofectamine 2000 (Invitrogen) when the cells had grown to approximately 80% confluence.

Immunofluorescence
Immunofluorescence was performed as previously described . HEK293T or HepG2 cells cultured on 12-mm coverslips were transfected with the indicated plasmids. After 24 h, cells were fixed with 4% paraformaldehyde and permeated with 0.5% Triton X-100. After cells were washed with PBST, they were blocked in 1% BSA and stained with primary antibodies, followed by staining with CoraLite 594 or 488-conjugated IgG secondary antibodies. Nuclei were stained with DAPI (Yesen Biotechnology, Shanghai, China). Fluorescence images were obtained and analyzed using a confocal microscope (FV3000, OLYMPUS).

Co-immunoprecipitation and immunoblot analysis
Co-immunoprecipitation was performed as previously described (Liu et al., 2022). In brief, cells were lysed in the lysis buffer containing protease and phosphatase inhibitor cocktail (Selleck, Houston, Texas, United States). For co-immunoprecipitation, lysates were incubated overnight with ANTI-FLAG ® M2 Affinity Gel (Sigma-Aldrich). Then, proteins were separated by SDS-PAGE and electro-transferred onto the PVDF membrane. The membrane-containing proteins were blocked for 1 h at room temperature with 5% BSA in PBST, followed by treatment with the indicated primary antibodies overnight at 4°C. Subsequently, blots were incubated with secondary antibody for 1 h at room temperature and visualized via enhanced chemiluminescence reagents (Millipore, Billerica, MA, United States) with ChemiDoc XRS + Molecular Imager software (Bio-Rad, Philadelphia, PA, United States).

Nuclear and cytoplasmic extraction
Nuclear and cytoplasmic extraction was performed as previously described . Briefly, cells were treated using the nuclear and cytoplasmic protein extraction kit (Beyotime, China) to prepare the nuclear and cytoplasmic fraction, according to the manufacturer's instructions. The purified cytoplasmic and nuclear fractions were subjected to immunoblots assay according to the standard procedures with the relevant antibodies.

Mitochondrial isolation
Total cells were lysed in the lysis buffer containing the protease and phosphatase inhibitor cocktail (Selleck, Houston, Texas, United States). Mitochondrial and cytosolic proteins were isolated using the Cell Mitochondrial Isolation Kit (Beyotime, China) and subjected to immunoblots according to the standard procedure with the relevant antibodies.

Virus infection
ALSV, isolated from a tick-bitted patient in Northeast China , was used to infect the target cells. Briefly, the cultured cells were infected with ALSV at a multiplicity of infection (MOI) of 10 diluted in serum-free DMEM, and incubated for 2 h at 37°C. After washing three times with PBS, the infected cells were cultured with DMEM containing 2% FBS for the indicated times.

Construction of expression vectors for ALSV proteins
In order to clarify the function of these ALSV-encoded proteins ( Figure 1A), we constructed Flag-tagged NSP1, NSP2, VP1a, VP1b, VP2, VP3, and VP4 proteins expression plasmids and transfected them into HEK293T cells. At 24 h post transfection (hpt), the whole cell lysis was suffered to immunoblot assay with anti-Flag antibody to confirm the expression of viral proteins. The observed bands of viral protein were close to the calculated molecular weights, but VP3 showed an obvious band above 130 KD and VP4 appeared multiple diffusely distributed bands ( Figure 1B), indicating that they may undergo multimerization and/or post-translational modification.

The ER localization of ALSV proteins in HEK293T cells
ER is the major cellular organelle that viruses need to usurp because it is a factory for producing viral proteins.
First, we examined the locations of ALSV proteins with ER. Calnexin is an ER-residing chaperone protein and serves as a marker protein of ER, which is essential for controlling the production of glycosylated proteins in the process of viral infection and immune escape (Kohli et al., 2021). To this end, HEK293T cells were fixed at 24 hpt and stained with anti-Flag in red fluorescence and anti-calnexin to visualize the ER in green fluorescence. Consequently, we found that all ALSV viral proteins were distributed in the cytoplasm and associated with ER ( Figure 2). The colors merge showed that the red fluorescence of NSP2 and VP4 was overlapping with the green fluorescence of calnexin, and the ImageJ line scan analysis found that the change trend of fluorescence intensity profiles of NSP2 and VP4 proteins was consistent with that of calnexin, which demonstrated that NSP2 and VP4 proteins significantly colocalized with ER. The changing trend of fluorescence intensity profiles of VP1b and VP2 proteins is consistent with Calnexin to some extent, indicating that VP1b and VP2 partially colocalized with ER. In addition to serving as the virus replication site, the ER is also central to the cytopathic effects and death observed in the infected cells (Monel et al., 2017(Monel et al., , 2019. We also found that the cells expressing NSP1, VP1a, and VP3 proteins shrank, in which the ER aggregated and distributed in perinuclear, indicating that these viral proteins may be involved in inducing cell death. In addition, the VP1b protein induced the formation of ER-derived large cytoplasmic vacuoles (Figure 2 the white arrow), which may be followed by "implosive" cell death (Monel et al., 2017). The results indicated that ALSV proteins either colocalized with ER or led to ER condensation and cell death.

The ER colocalization of ALSV proteins in HepG2 and Vero cells
Considering that viruses have different sensitivities to different tissue-sourcing cells, we used the same method to characterize the co-localization of ALSV proteins with ER in HepG2 and Vero cells. In HepG2 cells, the co-localization analysis demonstrated that NSP1 and VP2 were significantly co-localized with ER; NSP2, VP1b, and VP4 were partially co-localized with ER ( Figure 3). Different from HEK293T cells, VP4 and calnexin were distributed in both cytoplasm and nucleus of HepG2 cells and appeared as obvious overlapping dots in the nucleus. A small amount of VP1a was also found in nucleus; none of the viral proteins caused obvious cell shrinkage (Figure 3). In addition to VP1b, VP2 also induced ER-derived large cytoplasmic vacuoles in HepG2 cells (Figure 3 the white arrow).
In Vero cells, NSP2, VP1b, VP2, and VP4 co-localized with ER, all viral proteins distributed in the cytoplasm, and NSP1 induced ER aggregated and cellular shrinking, all of which are the same as HEK293T cells (Supplementary Figure S1) Figure S1). Collectively, these results indicated that viral proteins exhibited diverse subcellular distribution and induced various ER morphological changes in different deriving cells. The NSP2, VP1b, VP2, and VP4 proteins were all co-localized with ER of these cells.

VP1a and VP4 induce calnexin transfer into the nucleus in HepG2 cells
We found that VP1a and VP4 proteins entered the nucleus in the liver-deriving HepG2 cells but not in human embryonic kidney epithelial HEK293T cells and African green monkey kidney-deriving Vero cells. We transfected VP1a and VP4 plasmids into HepG2 or HEK293T cells for immunofluorescence staining (IFA) assay to further verify this phenomenon. The results Construction of expression vectors for ALSV proteins. (A) Schematic genome organization and putative proteomic maps of ALSV, the open reading frame (ORF) is shown in the green box. The transmembrane area is displayed in red box, and the signal peptide is displayed in blue box. (B) Immunoblot analysis of ALSV protein expression. The ALSV proteins plasmids were separately transfected into HEK293T cells, and the Flag-vec was used as a control. At 24 h post-transfection (hpt), cellular lysis was measured by immunoblot with Flag antibody, GAPDH as a loading control. The red arrows, respectively, indicate VP3 and VP4 bands which consistent with the calculated molecular weights.
Frontiers in Microbiology 05 frontiersin.org The ER subcellular localization of ALSV proteins in HEK293T cells. HEK293T cells were transfected with the plasmids expressing Flag-NSP1, NSP2, VP1a, VP1b, VP2, VP3, VP4 or vector. After 24 h, the cells were immune-stained with anti-Flag (red) and anti-calnexin (green) antibodies and counterstained with DAPI to examine chromosomes (blue). (Right) ImageJ line scan analysis of the intensity profiles of the viral protein (red) and the calnexin (green) along the plotted lines. The distance in the x-axis represents the length of the plotted lines, and the consistency of the change trend of intensity profiles is proportional to the degree of co-localization of the two proteins.
Frontiers in Microbiology 06 frontiersin.org The ER colocalization of ALSV proteins in HepG2 cells. HepG2 cells were transfected with the plasmids expressing Flag-tagged viral proteins or empty vector. After 24 h, cells were immune-stained with anti-Flag antibody (red), anti-calnexin antibody (green), and DAPI (blue). The ImageJ line scan analyses are shown on the right. The distance in the x-axis represents the length of the plotted lines, and the consistency of the change trend of intensity profiles is proportional to the degree of co-localization of the two proteins.
Frontiers in Microbiology 07 frontiersin.org were consistent with the previous discovery. In HepG2 cells, VP4 and calnexin entered the nucleus and accumulated as scattered circular spots. VP1a was less in the nucleus than cytoplasm while mediated calnexin transferred into the nucleus ( Figure 4A). In HEK293T cells, both VP1a and VP4 were only distributed in the cytoplasm ( Figure 4B). We then confirmed the localization of VP1a and VP4 in HepG2 and HEK293T cells by subcellular fractionation assays, and the results were in line with our expectations. VP1a and VP4 entered the nucleus and also induced calnexin into the nucleus in HepG2 cells ( Figure 4C), whereas they did not exist in the nucleus but still mediated calnexin transferring into the nucleus in HEK293T ( Figure 4D). Previously, we found that VP4 protein induced calnexin into the nucleus and co-localized with calnexin in HepG2 nucleus ( Figure 4A). Therefore, we speculated VP4 might mediate calnexin nucleartransferring through protein-protein interaction and detected whether VP4 interacted with calnexin by co-immunoprecipitation assay. Unexpectedly, neither VP1a nor VP4 bound to calnexin in HepG2 cells, whereas VP4 significantly interacted with calnexin in HEK293T cells ( Figure 4E). These results indicated that VP1a and VP4 mediated calnexin nuclear-transferring independent of these interactions in HepG2 cells, and VP4 interacted with and induced calnexin into the nucleus in HEK293T cells. The nucleus is a double-membrane-bound organelle that contains many protein-composed dot-like domains with specialized functions (subnuclear organelles), such as nucleus, Cajal bodies (CBs), and speckles [35,36]. VP4 protein is morphologically consistent with the punctate distribution of the nucleolus, so we want to explore whether VP4 co-localize with nucleolus and plays key roles in viral protein trafficking and ribosome assembly. The HepG2 or HEK293T cells were fixed at 24 hpt and stained with anti-Flag in green and anti-B23 to visualize the nucleolus in red, and the results showed that VP1a and VP4 did not co-localize with nucleolus in both HepG2 and HEK293T cells ( Figure 4F; Supplementary Figure S2), which indicates that VP1a and VP4 may associate with other domains in the nucleus, such as splicing compartments which is essential for gene splicing (Fu and Maniatis, 1992).

NSP1 reduced mitochondria quantity by inducing mitophagy
Mitochondria is an essential functional organelle, supplying energy to cells and participating in various important cellular processes, including cellular oxidative stress, calcium homeostasis, antiviral immune signaling, autophagy, and apoptosis (Vandecasteele et al., 2001;West et al., 2011;Bertero and Maack, 2018;Pickles et al., 2018). The viral infection leads to many physiological changes in host cells, including the induction of mitophagy (Zhang et al., 2018). We screened the subcellular localization of ALSV proteins to mitochondria with the marker protein Cox IV and found that with the exception of VP1b, the co-localization of other viral proteins with Cox IV was not shown ( Figure 5A; Supplementary Figure S3A). Consistent with the results of ER staining, HEK293T cells shrunk when ectopic expression of NSP1, VP1a, and VP3. Additionally, we found that NSP1 could attenuate or disappear mitochondrial fluorescence signal, indicating that NSP1 specifically reduced mitochondria quantity. Mitochondria were also partially reduced or agglutinated in the cells transfected with VP1a, VP1b, and VP3 ( Figure 5A). To determine whether these proteins cause mitochondrial reduction and whether they are involved in mitophagy, we first examined the effect of viral proteins on the expression of mitochondrial marker Cox IV. The results showed that NSP1, VP1a, VP1b, and VP3 proteins reduced the expression of Cox IV (Supplementary Figure S3B). Then, we transfected NSP1 into HEK293T cells and performed an immunoblot assay, and found NSP1 also inhibited the mitochondrial marker proteins TIM23 and TOM20 expression ( Figures 5B,C), which was consistent with the results of the immunofluorescent assay. Simultaneously, NSP1 decreased the expression of autophagy-related protein P62 ( Figure 5B) and increased the expression of microtubuleassociated protein 1 light chain 3B (LC3B; Figure 5C), indicating that NSP1 can induce mitochondrial autophagy. Moreover, NSP1 dose-dependently inhibited the expression of mitochondrial marker proteins Cox IV, TOM20, and TIM23, but not the ER marker calnexin ( Figure 5D). These results indicate that NSP1 reduced mitochondria quantity by inducing mitophagy.

ALSV infection induces calnexin nuclear transfer and mitophagy
We further validated our findings on a real infection model, and found that ALSV infection induced calnexin transferring into the nucleus of HepG2 cells using the immunofluorescence and cell fractionation analysis (Figures 6A,B). Additionally, ALSV infection inhibited the expression of mitochondrial marker proteins TOM20, TIM23 and Cox IV in both Hela and Vero cells, but ALSV did not affect expression of the ER marker calnexin ( Figure 6C). ALSV also promoted the autophagy-related protein LC3B-II accumulation and p62 degradation ( Figure 6D). These results indicated that ALSV infection induces mitochondrial autophagy.

Discussion
ALSV usurps and manipulates the host-cellular organelles to generate a cell environment that is permissive to viral replication. Clarifying the subcellular distributions of the viral proteins will benefit us in understanding the infection and pathogenic mechanisms. In this study, we characterized the subcellular localization of ALSV proteins and found that viral proteins exhibited diverse subcellular distribution and induced various ER morphological changes in HEK293T, HepG2, and Vero cells. The nuclear transfer and co-localization of VP4 and calnexin were Frontiers in Microbiology 09 frontiersin.org unique to HepG2 cells; NSP1 protein reduced mitochondria quantity by inducing mitophagy, and all the ALSV proteins did not co-localize with the Golgi apparatus (Supplementary Figure S4).
To justify why the viral proteins were located in the corresponding subcellular compartments, we predicted the subcellular localizations of viral proteins using the Frontiers in Microbiology 10 frontiersin.org state-of-the-art computational methods FUEL-mLoc (Wan et al., 2017), Virus-mPLoc (Shen and Chou, 2010), and mGOASVM (Wan et al., 2012), which found that the results obtained by different methods vary widely (Supplementary Table S1). The Virus-mPLoc prediction results of all ALSV viral proteins distributed in host cytoplasm were consistent with our results, and the VP3 colocalized with host cell membrane is reasonable because it has nine transmembrane domains, which may explain the relationship between the viral protein amino acid sequences and their subcellular locations. The mGOASVM prediction results of the non-structural proteins NSP1 and NSP2 colocalized with host endoplasmic reticulum also further support our findings.
The endoplasmic reticulum, an important site for protein folding and post-translational modification, is an essential organelle for viral replication and maturation (Bagchi, 2020). After flavivirus particles enter cells, the released RNA is translated at the rough ER into a single polyprotein, which is cleaved into structural and non-structural proteins. Then, the non-structural proteins orchestrate ER membrane to form the viral replicase complex where the viral particles are produced and budded into the ER lumen (Tardif et al., 2004;Yu et al., 2006;Ambrose and Mackenzie, 2011). The functional viral proteins of ZIKV are synthesized in ER and exhibit a unique subcellular distribution in host cells, of which the E, PrM, NS2a, and NS4a proteins are co-localize with ER (Hou et al., 2017;Coyaud et al., 2018). In this study, we found that ALSV NSP1, NSP2, VP1b, VP2, and VP4 proteins were co-localized with ER in various cells which may be due to the fact that ER is the synthesis site of these viral proteins, and speculate that the non-structural proteins NSP1 and NSP2 have involved in ER membrane-associated replication complexes formation as with other flaviviruses. Virus and viral proteins are localized to ER and induce the reticular architecture remodeling, creating an ideal viral replication environment, simultaneously introducing additional burden and ER stress to host cells (Oakes and Papa, 2015). The accumulation of misfolded viral proteins in ER lumen after ZIKV infection overwhelms ER protein folding capacity, leading to ER stress and triggering the unfolded protein response (UPR) which attempts to restore the disturbed homeostasis, once the restoration failing will initiate the apoptosis or paraptosis-like cell death Frontiers in Microbiology 11 frontiersin.org (Monel et al., 2017(Monel et al., , 2019Turpin et al., 2021). Here, we found NSP1, VP1a, and VP3 induced HEK293T cells shrinking and ER aggregating in perinuclear, suggesting that ALSV proteins may activate apoptosis response; VP1b, VP2, and NSP2 induced the form of ER-derived cytoplasmic vacuoles in various cells, indicating that ALSV proteins also activate paraptosis-like cell death.
Human hepatoma cells exhibit cell-type dependent morphology changes after flaviviruses infection, such as ZIKV induced a unique replication organelle (named zippered ER (zER) membranes) and bigger vesicles in infected hepatoma cells (Cortese et al., 2017). Here, we found that VP1a and VP4 entered the nucleus and mediated calnexin nucleartransferring independent of their interaction, which was unique to HepG2 cells and argued for a cell type-specific role. Calnexin, a calcium ion-binding protein localized on the ER membrane is essential for flavivirus replication and assembly complexes formation. The hepatitis B virus (HBV) envelope protein requires calnexin for proper folding and transport (Prange et al., 1999), and the production of infectious dengue virions is significantly reduced when calnexin is knockdown (Limjindaporn et al., 2009). Additionally, calnexin is found in the nucleus of hepatocytes (Czubryt et al., 1996), which may explain why calnexin is distributed in the nucleus of HepG2 cells. However, the mechanisms and functions of calnexin nuclear-transferring mediated by VP1a and VP4 are not clear, and further research is needed. Although VP1a and VP4 cannot enter the nucleus in HEK293T cells, they still induce calnexin into the nucleus, and the mechanism remains to be further studied.
We also screened the co-localization of ALSV proteins with the Golgi apparatus with the marker protein GM130, and found that no viral proteins was co-localize with the Golgi apparatus (Supplementary Figure S4). However, it cannot be ruled out that ALSV proteins interact with the Golgi apparatus by forming the newly synthesized virus particles which are transported from ER and Golgi apparatus to the cell surface (Hansen et al., 2017).
In summary, we found ALSV proteins distributed in the cytoplasm, with the exception of VP1a and VP4 expressed in HepG2 cells. The nuclear transfer and co-localization of VP4 and calnexin independent of their interaction was unique to HepG2 cells, and the underlying mechanism still requires further study. Additionally, ALSV proteins were either co-localized with ER or induced various ER morphological changes, and NSP1 protein reduced mitochondria quantity by inducing mitophagy. Our results provide new insight into the biological functions of ALSV proteins because the subcellular distribution of viral protein hints at its possible biological function there. However, our study utilizes the surrogate expression of viral proteins, which may not reflect the real expression level upon ALSV infection and the complex environment achieved by simultaneously expressing all viral proteins.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material; further inquiries can be directed to the corresponding authors.

Ethics statement
The studies involving human participants were reviewed and approved by the ethics committee of the Inner Mongolia General Forestry Hospital. The patients/participants provided their written informed consent to participate in this study.

Author contributions
QL, YZ, and ZH designed the experiments. YZ, PW, LL, BM, YH, ND, and HY conducted the experiments. Z-DW, LS, and ZH analyzed the data and conducted statistical analysis. YZ, MP, and QL drafted and revised the paper. All authors contributed to the article and approved the submitted version.

Funding
This work was supported by grant from the National Natural Science Foundation of China (81972873), the Pearl River Talent Plan in Guangdong Province of China (2019CX01N111), and the Scientific and Technological Research Projects of Guangzhou (202103000008).