- 1Laboratory of Computational Biology and Machine Intelligence, National Laboratory of Pattern Recognition, Institute of Automation, Chinese Academy of Sciences, Beijing, China
- 2School of Artificial Intelligence, University of Chinese Academy of Sciences, Beijing, China
- 3CAS Key Laboratory of Separation Science for Analytical Chemistry, Dalian Institute of Chemical Physics, Chinese Academy of Sciences, Dalian, China
- 4Tomas Lindahl Laboratory, The Seventh Affiliated Hospital, Sun Yat-Sen University, Shenzhen, China
Endoplasmic reticulum stress (ER stress) is a condition that is defined by abnormal accumulation of unfolded proteins. It plays an important role in maintaining cellular protein, lipid, and ion homeostasis. By triggering the unfolded protein response (UPR) under ER stress, cells restore homeostasis or undergo apoptosis. Chronic ER stress is implicated in many human diseases. Despite extensive studies on related signaling mechanisms, reliable image biomarkers for ER stress remain lacking. To address this deficiency, we have validated a morphological image biomarker for ER stress and have developed a deep learning-based assay to enable automated detection and analysis of this marker for screening studies. Specifically, ER under stress exhibits abnormal morphological patterns that feature ring-shaped structures called whorls (WHs). Using a highly specific chemical probe for unfolded and aggregated proteins, we find that formation of ER whorls is specifically associated with the accumulation of the unfolded and aggregated proteins. This confirms that ER whorls can be used as an image biomarker for ER stress. To this end, we have developed ER-WHs-Analyzer, a deep learning-based image analysis assay that automatically recognizes and localizes ER whorls similarly as human experts. It does not require laborious manual annotation of ER whorls for training of deep learning models. Importantly, it reliably classifies different patterns of ER whorls induced by different ER stress drugs. Overall, our study provides mechanistic insights into morphological patterns of ER under stress as well as an image biomarker assay for screening studies to dissect related disease mechanisms and to accelerate related drug discoveries. It demonstrates the effectiveness of deep learning in recognizing and understanding complex morphological phenotypes of ER.
Introduction
The endoplasmic reticulum (ER) is a continuous membrane-bound organelle network that plays key roles in protein synthesis and modification, lipid biogenesis, and ionic homeostasis (Friedman and Voeltz, 2011). Dysregulation of protein synthesis and modification caused by intra- and extra-cellular cues leads to excessive accumulation of unfolded proteins, which triggers ER stress (Lee et al., 2015; Schwarz and Blower, 2016). Unfolded protein response (UPR) often refers to the signal transduction pathway for cellular response to ER stress (Hetz et al., 2020). There are three main UPR transducers: inositol-requiring enzyme 1 (IRE1) (Tirasophon et al., 1998; Wang et al., 1998), protein kinase RNA-like ER kinase (PERK) (Harding et al., 1999), and activating transcription factor 6 (ATF6) (Harding et al., 2000; Ron and Walter, 2007). By orchestrating cellular processes such as mRNA splicing by endoribonuclease IRE1α (Tirasophon et al., 1998; Wang et al., 1998), translation attenuation by kinases PERK (Harding et al., 1999), and protein folding assistance by chaperone BiP (binding immunoglobulin protein) (Shimizu et al., 2017), the UPR engages different outputs to restore ER protein homeostasis under mild stress conditions or to activate apoptosis under chronic stress conditions (Yoon et al., 2012; Nie et al., 2017; Ramachandran et al., 2021).
ER stress is strongly implicated in the onset and progression of a wide range of human diseases, including neurodegenerative diseases, metabolic diseases, and cancer (Ozcan et al., 2004; Zhang et al., 2005; Wang and Kaufman, 2016). It can cause not only alterations of protein synthesis or folding but also deleterious cellular responses including accumulation of lipids and activation of autophagy. Modulating ER stress shows great potential in the treatment of these diseases. Several compounds, including IRE1
Changes in ER morphology correlate well with ER stress (Schuck et al., 2009; Mateus et al., 2018). The classical ER structure consists of a continuous envelope surrounding the nucleus in the perinuclear region and a polygonal network of interconnected tubules and sheets in the peripheral region (Friedman and Voeltz, 2011). In ER-stressed cells, ER membranes are compacted to form ER whorls (Nii et al., 1968; Snapp et al., 2003). Formation of ER whorls provides an effective structural response to prolonged ER stress (Xu et al., 2020). It accompanies the activation of PERK and works together with vesicle transport machinery such as ESCRT (endosomal complexes required for transport) and COPII (coat protein complex II) complex to counter-balance ER stress-induced protein translation and ER expansion (Bernales et al., 2006; Schuck et al., 2014; Schäfer et al., 2020; Xu et al., 2020). ER whorls have been observed in yeast and mammalian cells under ER stress activated by various stimuli, including drugs and herpes simplex virus infection (Nii et al., 1968; Schäfer et al., 2020). To use ER whorls as an image biomarker for phonotypic drug screening, reliable and automated detection is essential. However, detection of ER whorls within the complex and dense ER network morphology poses a substantial technical challenge (Kazama et al., 2018).
The past decade witnessed the rapid rise of deep learning as a transformative artificial intelligence technique that computes using deep neural networks (DNNs). It has achieved breakthrough performance in many challenging tasks of analyzing natural images (LeCun et al., 2015). It has also achieved breakthrough performance in analyzing cellular images that previously are considered intractable for traditional methods (Moen et al., 2019). Unlike traditional methods, which rely on manually designed features to represent phenotypes in images, deep learning models automatically learn phenotypic features through their supervised training (LeCun et al., 2015). Recently, for example, deep learning models such as ResNet (He et al., 2016) and DenseNet (Huang et al., 2017) have achieved great success in recognizing cell states (Godinez et al., 2017; Sommer et al., 2017) and protein subcellular localization patterns (Kraus et al., 2016; Pärnamaa and Parts, 2017).
In this study, by comparing ER morphology under normal conditions versus IRE1
Materials and Methods
Induction and Western Blot Analysis of ER Stress in HEK293T Cells
HEK293T cells were treated with DMSO or ER stress induction compounds Thapsigargin (Tg) or Dithiothreitol (DTT) for various lengths of time based on experimental designs, typically 6 h, then trypsinized, pelleted (800 × g, 4 min, room temperature) and lysed in RIPA buffer (Millipore) with protease inhibitor (Roche) and phosphatase inhibitor on ice. The supernatant was collected by centrifugation (13,000 × g, 15 min, 4°C) and the protein concentration was determined using BCA protein assay (Beyotime). Proteins were resolved by 12% SDS-PAGE, transferred to PVDF membranes, then blocked with 5% non-fat milk in TBST for 1 h at room temperature. The membranes were then washed with TBST (3 × 10 min) and incubated with anti-IRE1 (CST), anti-p-IRE1 (Abcam), or anti-Caspase 3 (CST) overnight at 4°C. Next, the membranes were washed with TBST (3 × 10 min) and incubated with goat-anti-rabbit-HRP or goat-anti-mouse-HRP in TBST for 1 h at room temperature. After washing with TBST for 3 times, the blots were developed using an ECL detection reagent.
Live Cell Imaging
High-resolution images of ER in live cells were acquired under two conditions, i.e., normal ER morphology in control cells and abnormal ER morphology (i.e., with ER whorls) in cells under induced ER stress. HEK293T cells labeled with ER marker GFP-sec61b were treated with DMSO, Tg or DTT for various lengths of time based on experimental designs, typically 6 h. To check abundance and location of ER proteins on whorls, BFP-KDEL was transiently expressed in HEK293T cells to label luminal proteins. And mCherry-Rtn4a, mCherry-ATL3, and GFP-REEP5 were transiently expressed in HEK293T cells to label ER morphology regulator proteins reticulon, atlastin, and REEP5, respectively. Treated cells were then imaged using conventional spinning disk confocal microscopy at ∼200 nm resolution (Nikon CSU-W1 under 100× and 1.45 NA, excitation wavelength: 488 nm, emission wavelength: 535 nm) or 3D-SIM at ∼70 nm resolution (Nikon N-SIM, 100× SR objective, excitation wavelength: 488 nm, emission wavelength: 535 nm).
Synthesis of the AIEgen Probe for Misfolded and Aggregated Proteins
A 250 ml round bottom flask was charged with isophorone (6.0 ml, 40 mmol), malononitrile (2.9 g, 44 mol) and piperidine (cat.). The mixture was heated to reflux and stirred for 30 h. After cooling to room temperature, the solution was slowly poured into water and the precipitated solid was filtered. Recrystallization from EtOH afforded S1 as a brown solid (2.1 g, 28%). To a 25 ml round bottom flask S1 (186.3 mg, 1.0 mmol), 4-(bis(2-hydroxyethyl)amino) benzaldehyde (313.8 mg, 1.0 mmol) and piperidine (cat.) were stirred in 10 ml EtOH at 85°C for 20 h. The solution was cooled to room temperature, and water was added to the solution. Then the mixture was extracted with DCM (3 × 20 ml). The organic phase was combined and dried over Na2SO4. After filtration and concentration in vacuo, the residue was purified via flash silica gel chromatography (10–40% EtOAc in hexane) to provide the compound AIEgen as violet solid (231.4 mg, 61.3%) and analyzed by 1H-NMR (400 MHz, CDCl3).
In vitro Thermal Shift Assay
AIEgen (25
ER-WHs-Analyzer
Image pre-processing and data augmentation – Each of the high-resolution ER images acquired contains multiple cells. Each single cell was first cropped from the acquired full-size images. The cropped images were resized into 256
Feature recognition module – It is the module used for recognizing and classifying ER morphological patterns, with or without ER whorls. Two types of representative DNNs, the ResNet (He et al., 2016) and DenseNet (Huang et al., 2017), were chosen as the backbone network of the feature recognition module. The recognition of ER whorls, present or absent, is a standard classification problem. The Softmax function was used to map output of the DNNs into a classification score:
In Eq. 1,
where
Training strategies - Two training strategies were used: training from scratch and finetuning from a pre-trained model, i.e., transfer learning (Tan et al., 2018), a widely used approach to stabilize the training of deep learning models. To this end, a large dataset was collected primarily from open-source microscopy images, named as CBMI-Extra, which includes ∼70k images from ∼120 classes. Some of the images in CBMI-Extra were acquired using imaging protocols similar as those for the ER images in this study. In this way, a deep learning model pre-trained using CBMI-Extra provided a sound starting point to stabilize its subsequent training using ER images.
Feature localization module – A deep learning model with good performance is expected to capture image features similarly as human experts. The feature visualization tool, Grad-CAM (Selvaraju et al., 2017), was used to check whether the DNN models can correctly recognize features of the ER whorls. Given an input image for a forward pass in a trained network, the feature visualization tool generates a class activation mapping (CAM) in the form of a heatmap, visualizing the importance of each location in the input image in terms of its contribution to the prediction of the network. In this way, the feature visualization tool Grad-CAM constitutes the first part of the feature localization module by identifying up to several hotspots (i.e., clusters with higher scores) in the heatmap. It has been observed that the identified hotspots in the ER images match well with the ER morphological features that a human expert would identify. Based on the hotspots in heatmaps, an image processing-based pipeline was used to detect locations of the regions of ER whorls. This pipeline, which consists of segmentation, instance labeling and bounding-box assignment, constitutes the second part of the feature localization module (Supplementary Figure S1).
Results
ER Stress Induces Morphological Deformation That Forms Whorls
We used ER stress activator Thapsigargin (Tg) to set up our experimental assay (Dibdiakova et al., 2019). Because Tg activation of ER stress increases the level of phosphorylated IRE1
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/767866/fcell-09-767866-HTML-r1/image_m/fcell-09-767866-g001.jpg)
FIGURE 1. Upregulated ER stress markers and abnormal ER shapes in HEK293T cells treated with Tg at different concentrations. (A) Immunoblot analysis of p-IRE1
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/767866/fcell-09-767866-HTML-r1/image_m/fcell-09-767866-g002.jpg)
FIGURE 2. Progressive changes of ER structure under ER stress over time. (A) Immunoblot analysis of p-IRE1
To investigate whether activated ER stress causes morphological changes, ER was labeled by stable expression of GFP-Sec61
To check whether ER whorls can be induced by drugs other than Tg, we treated HEK293T cells with several reagents reported in the literature, including MK-28, a PERK activator (Ganz et al., 2020); Palmitic Acid, a long-chain saturated fatty acid (Xu et al., 2020; Harada et al., 2002); Bufalin, a Na+/K+-ATPase inhibitor (Shen et al., 2014); and CB-5083, a p97 inhibitor (Bastola et al., 2016). These reagents are known to induce ER stress via different mechanisms (Ganz et al., 2020; Harada et al., 2002; Shen et al., 2014; Bastola et al., 2016). We generally found them to be less effective than DTT and Tg in inducing ER stress and performed all treatments at 1
Because prolonged ER stress may lead to cell death, we checked whether formation of whorls was caused by cell apoptosis. Using treatment of 10
Dynamic ER Whorl Formation and Structural Deformation
To examine the dynamic formation of ER whorls, we performed time-lapse live cell imaging. Along with the membrane expansion and aggregation under stress, initial formation of ER whorls started approximately 0.5–1 h after Tg treatment. Severe membrane deformation led to further local ER aggregation (Figure 2B). Over the next 1–5 h, the number of ER whorls continued to increase while existing whorls became more condensed. After 6 h, the number of whorls and their morphologies generally became stable. The whorls occupied most of the intracellular space, and ER network connections were mostly lost (Figure 2B). Quantitative analysis revealed that whorls appeared in ∼21.0% of cells after 1 h treatment (Figure 2C), and the percentage increased over time (Figure 2C). We also examined the three-dimensional structure of whorls (Figure 2D) and found that they were composed of warped ER membranes without ER tubules. Expansion and deformation of ER sheets could contribute to whorl formation. Importantly, canonical ER network connections were largely lost due to absence of ER tubules.
To check the abundance and location of ER proteins on the whorls, we examined fluorescently labeled Rtn4a, ATL3, and receptor accessory protein 5 (REEP5) under induced ER stress. We also checked ER luminal proteins by expressing fluorescently labeled lumen marker KDEL (Supplementary Figure S3A). Overall, Rtn4a rarely locates to ER whorls under Tg treatment. Luminal proteins labeled by KDEL locates to ER whorls in ∼10% of the treated cells, REEP5 locates to ER whorls in ∼20% of the treated cells but ATL3 locates to ER whorls in ∼80% of the treated cells (Supplementary Figure S3B). Overall, these results reveal differential abundance and location of ER proteins on ER whorls.
Unfolded and Aggregated Proteins are Attached to ER Whorls
The results so far have revealed tight connections between ER whorl formation and ER stress. Unfolded proteins are a key driver of ER stress (Braakman and Bulleid, 2011). To check whether ER whorls may be used as a reliable image biomarker for ER stress, we examined their relations with unfolded proteins using an AIEgen probe. It exhibited highly specific binding affinity to unfolded and aggregated proteins (Figures 3A–C and Supplementary Figure S4) such as E. coli dihydrofolate reductase (DHFR), a model protein for thermal shift assay that detects levels of protein aggregation (Figure 3C). The specific binding affinity of the AIEgen probe was further checked in mut-DHFR, sortase, human immunoglobulin, and superoxide dismutase (SOD1). Relative increase of fluorescence intensity ranging from ∼2 to ∼10 folds was detected, suggesting AIEgen can serve as a general probe for misfolded and aggregated proteins in vitro (Figure 3D). With illumination by AIEgen, misfolded and aggregated proteins were detected as puncta in vivo in live HEK293T cells and were found to be mobile or immobile (Figure 3E). The puncta were mostly rounded in shape under Structure Illumination Microscopy (SIM), with a diameter of ∼0.5
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/767866/fcell-09-767866-HTML-r1/image_m/fcell-09-767866-g003.jpg)
FIGURE 3. Detecting unfolded and aggregated protein in vitro and in vivo by AIEgen. (A) A cartoon illustration of the formation of misfolded oligomers and insoluble aggregated DHFR model proteins. (B) Structural features of the AIEgen. (C) Thermal shift assay using OD330 turbidity (black curve) and the fluorescence of the AIEgen (red curve). Fluorescence of AIEgen occurs slightly earlier than the formation of insoluble aggregates measured by OD330 turbidity assay, indicating that the fluorescence originates from misfolded oligomers. (D) Fluorescence of AIEgen detecting the folded verses misfolded and aggregated proteins. AIEgen (25
In addition to Tg treatment, ER stress can also be activated by Dithiothreitol (DTT) treatment, which causes protein misfolding and aggregation by blocking formation of disulfide bonds. Similar as under Tg treatment, DTT treatment at 10 mM for 6 h induced formation of ER whorls in 75.7% of treated cells (Figures 4A,B). In contrast, treatment of HEK293T cells with 1 mM of DTT for 6 h only induced formation of ER whorls in a small percentage (4.3%) of cells (Supplementary Figures S5A–C). When the concentration of DTT was increased to 3 mM, formation of ER whorls was detected in 45.3% of treated cells (Supplementary Figures S5A–C), indicating ER stress induced by DTT treatment drives whorl formation in a dose-dependent manner. The rate of whorl formation was significantly lower in control experiments, at 1.2 ± 0.2% (DTT 0 mM, 0 h), 1.2 ± 0.2% (DTT 0 mM, 6 h), and 0.8 ± 0.1% (DTT 10 mM, 0 h). In DTT treated cells, 81.7 ± 2.1% of them showed both puncta of unfolded proteins and ER whorls (Figures 4B,C), significantly higher than that of control cells (11.4 ± 1.6% in DTT 0 mM 6 h, 3.7 ± 0.3% in DTT 0 mM 0 h, and 4.1 ± 0.3% in DTT 10 mM 0 h; p < 0.0001; mean ± SEM; n = 235 cells from nine experiments). Interestingly, we noticed that the puncta of unfolded protein were tightly associated with deformed ER (Figure 4C). Specifically, 87.5 ± 0.2% of puncta of unfolded proteins were attached to ER whorls during their formation. They were either attached to outer surfaces of the whorls (68.3 ± 0.3%) or were wrapped inside the whorls (19.3 ± 0.2%) (Figures 4C,D) (average ± SEM, n = 355 whorls from 151 cells). Taken together, these results indicate that whorls can serve as a morphological image biomarker for ER stress. The tight attachment of unfolded proteins to ER whorls also suggests a direct role of whorls in isolating misfolded and aggregated proteins.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/767866/fcell-09-767866-HTML-r1/image_m/fcell-09-767866-g004.jpg)
FIGURE 4. ER whorls tether unfolded proteins under stress. (A) ER whorls and unfolded protein (UP) puncta appear only under stress. Scale bar = 10 μm. (B) Quantification of unfolded protein puncta and whorls under stress. A 6-h treatment of 10 mM DTT induced both unfolded protein and whorls formation at a final ratio of 81.7 ± 2.1% and 75.7 ± 2.0% (mean ± SD; from n = 235 cells), respectively. ****: p < 0.0001. (C) 3D view of colocalization of unfolded proteins and whorls. Unfolded proteins may be in contact with the outside surfaces of whorls (a) or exist inside the whorls (b). Scale bar = 1 μm. (D) Quantification of localization of unfold proteins with respect to whorls, 68.3 ± 0.3% (mean ± SD; from n = 355) of whorls have UPs attached to the membrane (ON), 19.3 ± 0.2% of whorls wrapped UPs inside (IN).
Development of ER-WHs-Analyzer for Automated Detection of ER Whorls
So far, we have shown that ER whorls can serve as an image biomarker for ER stress. To use it for screening studies, we developed a deep learning-based analysis assay, which we refer to as ER-WHs-Analyzer, for automated detection and analysis of the whorls. The overall workflow of ER-WHs-Analyzer is shown in Figure 5A. First, raw images were cropped from acquired full-size ER images. Then, cropped images were standardized in their sizes and enhanced in their quality through preprocessing. The feature recognition module (DNNs-based classification model) was trained using the preprocessed images along with their binary labels, i.e., WT (wildtype without whorls) or WHs (with whorls) (Figure 5C). The trained feature recognition module was then used to detect whether an ER image contains whorls. If ER whorls were detected, the feature localization module was used to generate a heatmap of features learned by the recognition module. Based on the heatmap, regions of whorls were localized by simple thresholding (Figure 5A; Supplementary Figure S1). Two types of representative DNNs, the ResNet (He et al., 2016) and DenseNet (Huang et al., 2017) (Supplementary Figures S6, S7), were used in this study. Their performance was compared using architectures with different depths including ResNet14, ResNet34, ResNet50, ResNet101, DenseNet121, DenseNet161, DenseNet169 and DenseNet201, where the numbers indicate the count of trainable model layers (Figure 5B; Supplementary Figure S6).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/767866/fcell-09-767866-HTML-r1/image_m/fcell-09-767866-g005.jpg)
FIGURE 5. Deep learning-based morphological classification of ER morphology. (A) Overall workflow of the proposed ER-WHs-Analyzer. (B) Performance evaluation of various DNNs-based classification models in the feature recognition module of ER-WHs-Analyzer. * indicates that the model was finetuned from a pre-trained model. Numbers on X axis indicate different depth of the deep learning models. Red arrow indicates the best performed model, the DenseNet-201 which consists of 201 weights-learnable layers. (C) Results of images from the ER-Stress-A dataset processed by ER-WHs-Analyzer-v1, including the raw images, their corresponding feature heat maps, their prediction scores, and their predicted labels. (D) Based on the phenotypic features learned, a t-SNE map was used to visualize the distribution of wildtype ER images (WT) and the ER images with whorls (WHs).
Training and Testing Deep Learning Models of ER-WHs-Analyzer
For training and testing of the feature recognition module in ER-WHs-Analyzer, we constructed a dataset that we referred to as ER-Stress-A. We collected 490 cell images of normal ER morphology (labeled as WT) or abnormal ER morphology with whorls (labeled as WHs). We split the images into a training set, a validation set, and a test set (Supplementary Table S1). The training and validation sets were used to finetune training hyper-parameters (Supplementary Table S2). The test set was used for standalone testing. Standard performance metrics for image classification were used, including F1 (F1 score), AUC (area under ROC curve), ACC (accuracy), Spc (specificity), Sen (sensitivity), Pre (precision). The validation set in ER-Stress-A was used to compare performance of different architectural configurations of backbone networks.
Validation results of different configurations of backbone networks are compared in Figure 5B. Several observations can be made. First, deeper networks generally provided better performance. For example, ResNet50 with 50 weight-learnable layers obtained a higher F1-score than ResNet34 and ResNet18. And DenseNet generally outperformed ResNet. Second, models pretrained with CBMI-Extra (indicated by *) generally outperformed models without pretraining. The benefit of pretraining was more pronounced for shallow models. For example, the ResNet34 pretrained with CBMI-Extra outperformed DenseNet121 without pretraining. Third, deeper networks such as ResNet101 and DenseNet201 were more prone to overfitting. This problem can be mitigated by finetuning. DenseNet 201 achieves overall the highest rate of recognizing ER whorls (98.78%) (Figure 5B). Other performance metrics are listed in Supplementary Table S5. We chose DenseNet201 with pretraining for the feature recognition module of ER-WHs-Analyzer. In standalone testing, the model achieved excellent performance with F1 = 98.27%, AUC = 99.65%, ACC = 98.54%, Spc = 99.16%, Sen = 97.70%, Pre = 98.84%. In model testing, the ratio of abnormal ER structures labeled by experts was 42.33%, while the ratio detected by ER-WHs-Analyzer is 42.16%, reaching a high level of agreement. After feature recognition, the feature localization module of ER-WHs-Analyzer determined positions of the regions of whorls based on a heatmap of learned feature and subsequent thresholding (Figure 5A).
Representative images from ER-Stress-A test set and their correct classification labels are shown in the second row and the first row of Figure 5C, respectively. The heatmaps of learning features are shown in the third row of Figure 5C, while their recognition results using the feature recognition module are shown in the fourth and fifth rows, respectively. Together, the results showed that ER-WHs-Analyzer can accurately recognize and localize individual ER whorls (Figure 5C; Supplementary Figure S8). We also used a t-SNE map (Maaten and Hinton, 2008) to visualize the representative features that our method learned to separate morphology of normal ER from ER with whorls (Figure 5D). We found that morphology of ER under normal condition and induced stress can be well differentiated.
Separating Different Sub-Phenotypes of ER Whorls Using ER-WHs-Analyzer
ER stress is induced by treatment of Tg and DTT via different mechanisms (Dibdiakova et al., 2019). Tg induces ER stress via interfering with calcium ion transport, while DTT induces ER stress by blocking the formation of disulfide bonds required for protein folding (Jiang et al., 2015). Differently from Tg or Tunicamycin (Tm), DTT is also considered as a robust pro-apoptotic ER stress inducer (Labunskyy et al., 2009). As expected, DTT treatment at a concentration of 3 mM caused formation of ER whorls. However, ER whorls induced by DTT (Figure 6A, second row, panels 5–6) showed morphological difference from those induced by Tg treatment (Figure 6A, second row, panels 3–4). Specifically, whorls induced by Tg treatment generally tend to be circular and small whereas whorls induced by DTT treatment tend to be elliptic and large. Cells treated with Tg generally tend to have a low number of whorls, typically one or two. Cells treated with DTT generally tend to have more whorls, typically two or more. We quantitatively analyzed the number and area of whorls in the cells treated by Tg and DTT using feature localization module (Supplementary Figure S1). Overall, DTT treatment induced an average of 2.43 ± 1.36 whorls per cell (mean ± SD; n = 112 cells) and an average area of whorls of 29.27 ± 16.75
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/767866/fcell-09-767866-HTML-r1/image_m/fcell-09-767866-g006.jpg)
FIGURE 6. Two sub-phenotypes of ER whorls were identified by ER-WHs-Analyzer-v2. (A) Results of images from the ER-Stress-B dataset processed by ER-WHs-Analyzer-v2, including the raw images, their corresponding feature heat maps, their prediction scores, and their predicted labels. (B) Based on the phenotypic features learned, a t-SNE map was used to visualize the distribution of wildtype ER images (WT), the ER images with WHs-Type I, and the ER images with WHs-Type II. (C) Double-blind testing results between human experts (H) and ER-WHs-Analyzer v2 (M). Purple bars show the distribution of WHs-Type I. Red bars show the distribution of WHs-Type II. Error bars indicate SEM, which was computed from 81 sampled images. **: p < 0.01; ***: p < 0.001; ****: p < 0.0001, ns: not significant.
To verify whether our ER-WHs-Analyzer can separate these morphological subphenotypes of ER whorls, we constructed another dataset that we refer to as ER-Stress-B. It contains three categories: Wildtype, WHs-Type I and WHs-Type II. WHs-Type I contains one or two whorls and WHs-Type II contains more than two whorls. ER-Stress-B contains a total of 1404 cell images. We partitioned all the images into a training set and a validation set (Supplementary Table S3). We used the same deep learning model and performance metrics as in previous experiments. For this multi-class recognition task, we replaced the previous two-class output layer with a three-class output layer and retrained the models. Again, we compared the performance of our models with and without pretraining with the CBMI-Extra dataset (Supplementary Table S4). From the results, we observed a similar trend as in previous two-class classification experiments. Overall, DenseNet outperforms ResNet, and models with pretraining (indicated by *) outperform those without. However, differently from previous experiments, DenseNet161 shows overall the best performance, whose F1 score reaches 98.86%. Detailed evaluation metrics are compared in Supplementary Table S6. We chose DenseNet161 with pretraining for the feature recognition module and referred to the overall assay as ER-WHs-Analyzer v2.
Next, we checked whether ER-WHs-Analyzer v2 can reliably separate different subphenotypes of ER whorls. Figure 6A shows the results, which confirms that ER morphologies of WHs-Type I and WHs-Type II can be reliably differentiated. The t-SNE map (Figure 6B) shows that ER-WHs-Analyzer v2 has learned representative features to distinguish the sub-phenotypes of ER whorls.
Finally, we performed a double-blind experiment to compare ER-WHs-Analyzer v2 against human experts in classifying ER images acquired under control, Tg treatment, and DTT treatment respectively. Based on the classification results, we counted the numbers of the cells belonging to different sub-phenotypes for the control group and the two experimental groups. From the results, we have observed clear distribution differences of these three groups of cells (Figure 6C). According to classification by human experts, the control group contains 96.1 ± 7.3% WT, 3.9 ± 1.2% WHs-Type I and 0.0 ± 0.0% WHs-Type II. The Tg treatment group contains more WHs-Type I, with 78.7 ± 4.8% WT, 14.6 ± 3.5% WHs-Type I, and 6.8 ± 2.1% WHs-Type II (Figure 6C). In contrast, the DTT treatment group contains more WHs-Type II, with 13.6 ± 6.8% WT, 17.0 ± 4.2% WHs-Type I, 69.3 ± 7.1% WHs-Type II (Figure 6C). Classification results by ER-WHs-Analyzer v2 generally matched those by human experts (p > 0.07, n > 81 from seven experiments) (Figure 6C). Together, these experimental results suggested that ER-WHs-Analyzer v2 can reliably detect morphological differences of ER whorls induced by different treatments.
Discussion
In this study, we examined morphological patterns of ER under stress and identified ER whorls as an image biomarker of ER stress for screening studies. ER whorls are found in both yeast and mammalian cells, and their formation is considered as an integral part of cellular response to ER stress (Bernales et al., 2006; Schäfer et al., 2020; Xu et al., 2020), an important target of drug development studies for treatment of cancer as well as metabolic and neurodegenerative diseases. Although a wide variety of chemical and genetic tools for assessing ER stress have been developed (Sicari et al., 2020), image biomarkers for efficient phenotypic screening have been lacking. Our study fills this gap.
Monitoring ER stress by detecting and analyzing ER whorls carry some important advantages. Canonical ER stress detection is carried out by Western blot analysis of the expression or phosphorylation of ER stress modulators or their transcription by quantitative reverse transcription polymerase chain reaction (RT-PCR) (Haynes and Wiseman, 2018). These are endpoint assays that are laborious. More importantly, they cannot measure ER stress in single cells. Another approach to detect ER stress is to monitor the transcription reporters of molecules key to UPR, such as XBP1 and ATF4, by fluorescence. But it is not a real-time approach because of the time delay between initiation of transcription and the illumination of reporter proteins. In comparison, our study shows that ER morphological changes correlate well ER stress, and ER whorls can be used as an image biomarker to detect ER stress. They appear as early as 1 h after stress induction, indicating that they respond quickly to ER stress. Furthermore, high-resolution live-cell imaging of whorls makes real-time and single-cell level monitoring of ER stress possible. Our study combines automated high-resolution ER microscopy with deep learning-based analysis using ER-WHs-Analyzer to achieve high-throughput observation and quantification. This automated monitoring and analysis assay can be used as an effective tool for screening or validation of ER stress-related targets or drugs. It can reliably separate different sub-phenotypes of ER morphology under stress induced by different drugs. It can also be trained for use with other specific image biomarkers for ER stress. However, experimental settings of our study such as concentration and duration of Tg and DTT treatment remain be further optimized.
The mechanism regulating ER whorl formation is complicated. Although the causal relation between formation of ER whorls and COPII, ESCRT has been established, the relation between formation of ER whorls and PERK signaling remains unclear (Schäfer et al., 2020; Xu et al., 2020). Our study shows that formation of ER whorls accompanies activation of IRE1
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Author Contributions
YG, WL and GY designed the study and wrote the manuscript. NL contributed to the design of the study. DS and YL synthesized AIEgen probe. YFZ, YTZ and YG developed ER-WHs-Analyzer and analyzed the data. JL, YY and WL performed wet-lab experiments.
Funding
This work was supported in part by research grants from the Natural Science Foundation of China (Nos. 91954201, 31971289, and 32101216), the Chinese Academy of Sciences (No. 292019000056), University of Chinese Academy of Sciences (No. 115200M001), and the Beijing Municipal Science & Technology Commission (No. 5202022).
Conflict of Interest
A patent application related to this study has been approved by the China National Intellectual Property Administration. The number of the approved patent is ZL 2021 1 0150528.7.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.767866/full#supplementary-material
References
Axten, J. M., Medina, J. R., Feng, Y., Shu, A., Romeril, S. P., Grant, S. W., et al. (2012). Discovery of 7-Methyl-5-(1-{[3-(trifluoromethyl)phenyl]acetyl}-2,3-Dihydro-1h-Indol-5-Yl)-7h-Pyrrolo[2,3-D]pyrimidin-4-Amine (GSK2606414), a Potent and Selective First-In-Class Inhibitor of Protein Kinase R (PKR)-like Endoplasmic Reticulum Kinase (PERK). J. Med. Chem. 55 (16), 7193–7207. doi:10.1021/jm300713s
Bastola, P., Neums, L., Schoenen, F. J., and Chien, J. (2016). VCP Inhibitors Induce Endoplasmic Reticulum Stress, Cause Cell Cycle Arrest, Trigger Caspase-Mediated Cell Death and Synergistically Kill Ovarian Cancer Cells in Combination with Salubrinal. Mol. Oncol. 10 (10), 1559–1574. doi:10.1016/j.molonc.2016.09.005
Bernales, S., McDonald, K. L., and Walter, P. (2006). Autophagy Counterbalances Endoplasmic Reticulum Expansion during the Unfolded Protein Response. Plos Biol. 4 (12), e423. doi:10.1371/journal.pbio.0040423
Braakman, I., and Bulleid, N. J. (2011). Protein Folding and Modification in the Mammalian Endoplasmic Reticulum. Annu. Rev. Biochem. 80, 71–99. doi:10.1146/annurev-biochem-062209-093836
Dibdiakova, K., Saksonova, S., Pilchova, I., Klacanova, K., Tatarkova, Z., and Racay, P. (2019). Both thapsigargin- and Tunicamycin-Induced Endoplasmic Reticulum Stress Increases Expression of Hrd1 in IRE1-dependent Fashion. Neurol. Res. 41 (2), 177–188. doi:10.1080/01616412.2018.1547856
Friedman, J. R., and Voeltz, G. K. (2011). The ER in 3D: a Multifunctional Dynamic Membrane Network. Trends Cell Biol. 21 (12), 709–717. doi:10.1016/j.tcb.2011.07.004
Ganz, J., Shacham, T., Kramer, M., Shenkman, M., Eiger, H., Weinberg, N., et al. (2020). A Novel Specific PERK Activator Reduces Toxicity and Extends Survival in Huntington's Disease Models. Sci. Rep. 10 (1), 6875. doi:10.1038/s41598-020-63899-4
Ghosh, R., Wang, L., Wang, E. S., Perera, B. G. K., Igbaria, A., Morita, S., et al. (2014). Allosteric Inhibition of the IRE1α RNase Preserves Cell Viability and Function during Endoplasmic Reticulum Stress. Cell 158 (3), 534–548. doi:10.1016/j.cell.2014.07.002
Godinez, W. J., Hossain, I., Lazic, S. E., Davies, J. W., and Zhang, X. (2017). A Multi-Scale Convolutional Neural Network for Phenotyping High-Content Cellular Images. Bioinformatics 33 (13), 2010–2019. doi:10.1093/bioinformatics/btx069
Han, D., Lerner, A. G., Vande Walle, L., Upton, J.-P., Xu, W., Hagen, A., et al. (2009). IRE1α Kinase Activation Modes Control Alternate Endoribonuclease Outputs to Determine Divergent Cell Fates. Cell 138 (3), 562–575. doi:10.1016/j.cell.2009.07.017
Harada, H., Yamashita, U., Kurihara, H., Fukushi, E., Kawabata, J., and Kamei, Y. (2002). Antitumor Activity of Palmitic Acid Found as a Selective Cytotoxic Substance in a marine Red Alga. Anticancer Res. 22 (5), 2587–2590.
Harding, H. P., Novoa, I., Zhang, Y., Zeng, H., Wek, R., Schapira, M., et al. (2000). Regulated Translation Initiation Controls Stress-Induced Gene Expression in Mammalian Cells. Mol. Cel 6 (5), 1099–1108. doi:10.1016/s1097-2765(00)00108-8
Harding, H. P., Zhang, Y., and Ron, D. (1999). Protein Translation and Folding Are Coupled by an Endoplasmic-Reticulum-Resident Kinase. Nature 397 (6716), 271–274. doi:10.1038/16729
Haynes, C. M., and Wiseman, R. L. (2018). “Coordinating Organismal Physiology through the Unfolded Protein Response Preface,” in Coordinating Organismal Physiology through the Unfolded Protein Response (Cham, Switzerland: Springer), 414, V–Vi .
He, K., Zhang, X., Ren, S., and Sun, J. (2016). “Deep Residual Learning for Image Recognition,” in Proceedings of the IEEE conference on computer vision and pattern recognition, Las Vegas, NV, USA, 27-30 June 2016 (IEEE). doi:10.1109/cvpr.2016.90
Hetz, C., Zhang, K., and Kaufman, R. J. (2020). Mechanisms, Regulation and Functions of the Unfolded Protein Response. Nat. Rev. Mol. Cell Biol. 21 (8), 421–438. doi:10.1038/s41580-020-0250-z
Huang, G., Liu, Z., Van Der Maaten, L., and Weinberger, K. Q. (2017). “Densely Connected Convolutional Networks,” in Proceedings of the IEEE conference on computer vision and pattern recognition, Honolulu, HI, USA, 21-26 July 2017 (IEEE). doi:10.1109/cvpr.2017.243
Jiang, H., Zou, J., Zhang, H., Fu, W., Zeng, T., Huang, H., et al. (2015). Unfolded Protein Response Inducers Tunicamycin and Dithiothreitol Promote Myeloma Cell Differentiation Mediated by XBP-1. Clin. Exp. Med. 15 (1), 85–96. doi:10.1007/s10238-013-0269-y
Kazama, H., Hiramoto, M., Miyahara, K., Takano, N., and Miyazawa, K. (2018). Designing an Effective Drug Combination for ER Stress Loading in Cancer Therapy Using a Real-Time Monitoring System. Biochem. Biophysical Res. Commun. 501 (1), 286–292. doi:10.1016/j.bbrc.2018.05.001
Kitakaze, K., Taniuchi, S., Kawano, E., Hamada, Y., Miyake, M., Oyadomari, M., et al. (2019). Cell-based HTS Identifies a Chemical Chaperone for Preventing ER Protein Aggregation and Proteotoxicity. Elife 8, e43302. doi:10.7554/elife.43302
Kraus, O. Z., Ba, J. L., and Frey, B. J. (2016). Classifying and Segmenting Microscopy Images with Deep Multiple Instance Learning. Bioinformatics 32 (12), i52–i59. doi:10.1093/bioinformatics/btw252
Labunskyy, V. M., Yoo, M.-H., Hatfield, D. L., and Gladyshev, V. N. (2009). Sep15, a Thioredoxin-like Selenoprotein, Is Involved in the Unfolded Protein Response and Differentially Regulated by Adaptive and Acute ER Stresses. Biochemistry 48 (35), 8458–8465. doi:10.1021/bi900717p
LeCun, Y., Bengio, Y., and Hinton, G. (2015). Deep Learning. Nature 521 (7553), 436–444. doi:10.1038/nature14539
Lee, W., Yoo, W., and Chae, H. (2015). ER Stress and Autophagy. Curr. Mol. Med. 15 (8), 735–745. doi:10.2174/1566524015666150921105453
Maaten, L. v. d., and Hinton, G. (2008). Visualizing Data Using T-SNE. J. Machine Learn. Res. 9 (Nov), 2579–2605.
Mateus, D., Marini, E. S., Progida, C., and Bakke, O. (2018). Rab7a Modulates ER Stress and ER Morphology. Biochim. Biophys. Acta (Bba) - Mol. Cell Res. 1865 (5), 781–793. doi:10.1016/j.bbamcr.2018.02.011
Moen, E., Bannon, D., Kudo, T., Graf, W., Covert, M., and Van Valen, D. (2019). Deep Learning for Cellular Image Analysis. Nat. Methods 16 (12), 1233–1246. doi:10.1038/s41592-019-0403-1
Moffat, J. G., Vincent, F., Lee, J. A., Eder, J., and Prunotto, M. (2017). Opportunities and Challenges in Phenotypic Drug Discovery: an Industry Perspective. Nat. Rev. Drug Discov. 16 (8), 531–543. doi:10.1038/nrd.2017.111
Nie, J., Liu, A., Tan, Q., Zhao, K., Hu, K., Li, Y., et al. (2017). AICAR Activates ER Stress-dependent Apoptosis in Gallbladder Cancer Cells. Biochem. Biophysical Res. Commun. 482 (2), 246–252. doi:10.1016/j.bbrc.2016.11.050
Nii, S., Morgan, C., and Rose, H. M. (1968). Electron Microscopy of Herpes Simplex Virus. J. Virol. 2 (5), 517–536. doi:10.1128/jvi.2.5.517-536.1968
Nunes, A. F., Amaral, J. D., Lo, A. C., Fonseca, M. B., Viana, R. J. S., Callaerts-Vegh, Z., et al. (2012). TUDCA, a Bile Acid, Attenuates Amyloid Precursor Protein Processing and Amyloid-β Deposition in APP/PS1 Mice. Mol. Neurobiol. 45 (3), 440–454. doi:10.1007/s12035-012-8256-y
Ozcan, U., Cao, Q., Yilmaz, E., Lee, A. H., Iwakoshi, N. N., Ozdelen, E., et al. (2004). Endoplasmic Reticulum Stress Links Obesity, Insulin Action, and Type 2 Diabetes. Science 306 (5695), 457–461. doi:10.1126/science.1103160
Pärnamaa, T., and Parts, L. (2017). Accurate Classification of Protein Subcellular Localization from High-Throughput Microscopy Images Using Deep Learning. G3 (Bethesda) 7 (5), 1385–1392. doi:10.1534/g3.116.033654
Ramachandran, G., Moharir, S. C., Raghunand, T. R., and Swarup, G. (2021). Optineurin Modulates ER Stress-Induced Signaling Pathways and Cell Death. Biochem. Biophysical Res. Commun. 534, 297–302. doi:10.1016/j.bbrc.2020.11.091
Ron, D., and Walter, P. (2007). Signal Integration in the Endoplasmic Reticulum Unfolded Protein Response. Nat. Rev. Mol. Cell Biol 8 (7), 519–529. doi:10.1038/nrm2199
Schäfer, J. A., Schessner, J. P., Bircham, P. W., Tsuji, T., Funaya, C., Pajonk, O., et al. (2020). ESCRT Machinery Mediates Selective Microautophagy of Endoplasmic Reticulum in Yeast. EMBO J. 39 (2), e102586. doi:10.15252/embj.2019102586
Schuck, S., Gallagher, C. M., and Walter, P. (2014). ER-phagy Mediates Selective Degradation of Endoplasmic Reticulum Independently of the Core Autophagy Machinery. J. Cell Sci 127 (Pt 18), 4078–4088. doi:10.1242/jcs.154716
Schuck, S., Prinz, W. A., Thorn, K. S., Voss, C., and Walter, P. (2009). Membrane Expansion Alleviates Endoplasmic Reticulum Stress Independently of the Unfolded Protein Response. J. Cell Biol 187 (4), 525–536. doi:10.1083/jcb.200907074
Schwarz, D. S., and Blower, M. D. (2016). The Endoplasmic Reticulum: Structure, Function and Response to Cellular Signaling. Cell. Mol. Life Sci. 73 (1), 79–94. doi:10.1007/s00018-015-2052-6
Selvaraju, R. R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., and Batra, D. (2017). “Grad-cam: Visual Explanations from Deep Networks via Gradient-Based Localization,” in Proceedings of the IEEE international conference on computer vision, Venice, Italy, 25 December 2017 (IEEE). doi:10.1109/iccv.2017.74
Shen, S., Zhang, Y., Wang, Z., Zhang, R., and Gong, X. (2014). Bufalin Induces the Interplay between Apoptosis and Autophagy in Glioma Cells through Endoplasmic Reticulum Stress. Int. J. Biol. Sci. 10 (2), 212–224. doi:10.7150/ijbs.8056
Shimizu, A., Kaira, K., Yasuda, M., Asao, T., and Ishikawa, O. (2017). Clinical and Pathological Significance of ER Stress Marker (BiP/GRP78 and PERK) Expression in Malignant Melanoma. Pathol. Oncol. Res. 23 (1), 111–116. doi:10.1007/s12253-016-0099-9
Sicari, D., Delaunay‐Moisan, A., Combettes, L., Chevet, E., and Igbaria, A. (2020). A Guide to Assessing Endoplasmic Reticulum Homeostasis and Stress in Mammalian Systems. FEBS J. 287 (1), 27–42. doi:10.1111/febs.15107
Snapp, E. L., Hegde, R. S., Francolini, M., Lombardo, F., Colombo, S., Pedrazzini, E., et al. (2003). Formation of Stacked ER Cisternae by Low Affinity Protein Interactions. J. Cell Biol 163 (2), 257–269. doi:10.1083/jcb.200306020
Sommer, C., Hoefler, R., Samwer, M., and Gerlich, D. W. (2017). A Deep Learning and novelty Detection Framework for Rapid Phenotyping in High-Content Screening. MBoC 28 (23), 3428–3436. doi:10.1091/mbc.e17-05-0333
Swinney, D. C., and Anthony, J. (2011). How Were New Medicines Discovered? Nat. Rev. Drug Discov. 10 (7), 507–519. doi:10.1038/nrd3480
Tan, C., Sun, F., Kong, T., Zhang, W., Yang, C., and Liu, C. (2018). “A Survey on Deep Transfer Learning,” in International conference on artificial neural networks, Alghero, Italy, 27 September 2018 (Springer). doi:10.1007/978-3-030-01424-7_27
Tirasophon, W., Welihinda, A. A., and Kaufman, R. J. (1998). A Stress Response Pathway from the Endoplasmic Reticulum to the Nucleus Requires a Novel Bifunctional Protein Kinase/endoribonuclease (Ire1p) in Mammalian Cells. Genes Dev. 12 (12), 1812–1824. doi:10.1101/gad.12.12.1812
Wang, L., Perera, B. G. K., Hari, S. B., Bhhatarai, B., Backes, B. J., Seeliger, M. A., et al. (2012). Divergent Allosteric Control of the IRE1α Endoribonuclease Using Kinase Inhibitors. Nat. Chem. Biol. 8 (12), 982–989. doi:10.1038/nchembio.1094
Wang, M., and Kaufman, R. J. (2016). Protein Misfolding in the Endoplasmic Reticulum as a Conduit to Human Disease. Nature 529 (7586), 326–335. doi:10.1038/nature17041
Wang, X.-Z., et al. (1998). Cloning of Mammalian Ire1 Reveals Diversity in the ER Stress Responses. EMBO j 17 (19), 5708–5717. doi:10.1093/emboj/17.19.5708
Xu, F., Du, W., Zou, Q., Wang, Y., Zhang, X., Xing, X., et al. (2020). COPII Mitigates ER Stress by Promoting Formation of ER Whorls. Cell Res. 42, 141–156. doi:10.1038/s41422-020-00416-2
Yoon, S. O., Park, D. J., Ryu, J. C., Ozer, H. G., Tep, C., Shin, Y. J., et al. (2012). JNK3 Perpetuates Metabolic Stress Induced by Aβ Peptides. Neuron 75 (5), 824–837. doi:10.1016/j.neuron.2012.06.024
Zhang, K., Wong, H. N., Song, B., Miller, C. N., Scheuner, D., and Kaufman, R. J. (2005). The Unfolded Protein Response Sensor IRE1α Is Required at 2 Distinct Steps in B Cell Lymphopoiesis. J. Clin. Invest. 115 (2), 268–281. doi:10.1172/jci200521848
Keywords: ER stress, morphological classification, image biomarker, deep learning, homeostasis
Citation: Guo Y, Shen D, Zhou Y, Yang Y, Liang J, Zhou Y, Li N, Liu Y, Yang G and Li W (2022) Deep Learning-Based Morphological Classification of Endoplasmic Reticulum Under Stress. Front. Cell Dev. Biol. 9:767866. doi: 10.3389/fcell.2021.767866
Received: 31 August 2021; Accepted: 31 December 2021;
Published: 21 January 2022.
Edited by:
Jianguang Ji, Lund University, SwedenReviewed by:
Stephen H. Howell, Iowa State University, United StatesRossella Venditti, Telethon Institute of Genetics and Medicine (TIGEM), Italy
Copyright © 2022 Guo, Shen, Zhou, Yang, Liang, Zhou, Li, Liu, Yang and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ge Yang, ge.yang@ia.ac.cn; Wenjing Li, wenjing.li@ia.ac.cn
†These authors have contributed equally to this work and share first authorship