Super resolution microscopy and deep learning identify Zika virus reorganization of the endoplasmic reticulum

The endoplasmic reticulum (ER) is a complex subcellular organelle composed of diverse structures such as tubules, sheets and tubular matrices. Flaviviruses such as Zika virus (ZIKV) induce reorganization of ER membranes to facilitate viral replication. Here, using 3D super resolution microscopy, ZIKV infection is shown to induce the formation of dense tubular matrices associated with viral replication in the central ER. Viral non-structural proteins NS4B and NS2B associate with replication complexes within the ZIKV-induced tubular matrix and exhibit distinct ER distributions outside this central ER region. Deep neural networks trained to distinguish ZIKV-infected versus mock-infected cells successfully identified ZIKV-induced central ER tubular matrices as a determinant of viral infection. Super resolution microscopy and deep learning are therefore able to identify and localize morphological features of the ER and allow for better understanding of how ER morphology changes due to viral infection.


ZIKV infection induces the formation of a dense tubular matrix in the CER.
In order to study ER morphology following ZIKV infection, we generated stable U87 glioblastoma cells transfected with either ER monomeric oxidizing environment-optimized green fluorescent protein (ERmoxGFP), a lumenal ER reporter containing the bovine prolactin signal sequence and KDEL ER retention sequence linked to inert, monomeric, cysteine-less moxGFP 29 , or the membrane-associated ER reporter Sec61β-GFP. U87 cells stably expressing the ER reporters were then infected with ZIKV strain PRVABC59 (Puerto Rico 2015) at a multiplicity of infection (MOI) of 1 for 48 h. Cells were fixed with 3% paraformaldehyde/0.2% glutaraldehyde to preserve ER architecture 6,7,30,31 and labeled for dsRNA, a marker for ZIKV replication factories 19 . Maximum projections of 3D STED image stacks show high intensity ERmoxGFP and Sec61β-GFP labeling in a CER region and low intensity labeling in PER tubules in mock-infected cells (Fig. 1A), as reported previously by diffraction limited confocal microscopy 3 . Upon ZIKV infection, the CER reorganized to form an intensely labeled crescent-shaped region surrounding a lower intensity perinuclear region that overlapped extensively with dsRNA (Fig. 1A). Imaris Bitplane software fragments the ER into distinct segments that can then be analyzed for different features, including reporter density, segment Z-height and segment overlap with other labels, such as dsRNA. Densitybased segmentation of the ERmoxGFP-and Sec61β-GFP-labelled ER of ZIKV-infected cells showed that the higher density crescent-shaped CER region exhibited significant overlap with dsRNA-positive ER structures relative to the rest of the ER (Fig. 1B). This suggests that ZIKV dsRNA associates with an ER region of high density for both lumenal and membrane ER reporters.
We then investigated the relationship between the dense ZIKV-induced crescent-like region and the CER and PER. Examining the ER based on Z-height of ER segments showed that PER tubules were present in 3-5 serial sections of 210 nm each while CER regions were abundantly present across 7-10 serial sections ( Fig. 2A). After segmenting the ER in both ERmoxGFP and Sec61β-GFP labeled cells, average max Z-height of the PER tubular regions were 0.95 μm and for CER regions 1.8 μm ( Fig. 2A). Additionally, no PER segments reached a height greater than 1.26 μm (Fig. 2B). Based on this clear height separation of ER segments, ER segments with a height above 1.26 μm were therefore classified as being part of the CER mask and ER segments with a height below 1.26 μm as the PER mask. This height-based segmentation of ER segments effectively segmented the ER into CER-and PER-like regions, however it should be noted that the segmentation only approximated the CER and PER as the PER mask included short segments under the nucleus (Fig. 2B). Based on this ER segment heightbased classification of CER and PER, the CER was found to present a two-fold increase in volume and increased density relative to the PER (Fig. 2C). ER density in the height-based CER was elevated relative to PER and ZIKV infection showed minimal impact on the relative volume or density of CER and PER regions, except for Sec61β-GFP labeled cells where CER volume was reduced during ZIKV-infected compared to mock-infection (Fig. 2C).
Overlaying the CER and PER masks with the dsRNA-positive ER mask showed that the dsRNA-positive ER (> 80% volume/volume) is predominantly included within the CER mask (Fig. 3A). Indeed, only 10% of PER volume contains dsRNA while 35% of CER volume contains dsRNA for both ER reporters (Fig. 3A). The CER was composed of a convoluted network of tubules for both the ERmoxGFP-and Sec61β-labeled ER (Fig. 3B). 3D reconstructions confirmed that these regions were predominantly tubular with a few small sheet-like structures, morphologically similar to the tubular matrix of peripheral ER sheets 7 . 3D voxel-based visualization and quantification showed that the density of ER tubular structures in the dsRNA-positive ER is higher, for both the ERmoxGFP or Sec61β-GFP ER reporters, than in the dsRNA-negative CER regions of ZIKV-infected cells or the CER of mock-infected cells (Fig. 3B). The higher ER reporter density reflects reduced spacing between tubules in the dsRNA-positive ER, suggesting that ZIKV infection induces tubular matrix reorganization in a subdomain of the CER in U87 cells. To validate our finding of central ER tubular matrices in a more relevant model of ZIKV-induced microcephaly, we then infected cerebral brain organoids with ZIKV 32 . EM analysis Figure 1. 3D STED microscopy reveals ZIKV-induced ER reorganization in human U87 glioblastoma cells. (A) ERmoxGFP or Sec61β-GFP stably transfected U87 cells were mock-infected or infected for 48 h with the PRVABC59 2015 ZIKV strain (MOI = 1). ER reporter GFP and immunostained dsRNA-labeled ZIKV replication factories were imaged by 3D STED microscopy. (B) Fluorescence intensity of ERmoxGFP of infected cells using a spectrum heat map and a segmentation mask of the ER that colocalizes with dsRNA (grey), both generated on Imaris × 64 9.2.1 (Imaris), are depicted. Yellow squares in the panels indicate the magnified ROIs shown in the adjacent panels. Quantification of the mean normalized ER density ((Intensity sum of mask/total cell intensity sum)/ (volume sum of mask/ total cell volume sum)) was performed for both dsRNA-positive and dsRNA-negative ERmoxGFP and Sec61β-GFP in ZIKV-infected cells by Imaris segmentation. Scale bar: 10 microns. 5 cells per biological replicate (N = 3). Statistics were done using unpaired Student's T tests: ** = P < 0.01, Error bars represent SEM. and dsRNA (red) labeling in 3 × 3 μm ROIs of the ZIKV-infected dsRNA-positive (red) and dsRNA-negative (yellow) CER and mock-infected (green) CER are shown above Imaris 3D surface rendering of 1 × 1 μm regions of the above ROIs. Graph shows mean normalized ER density for each of the three CER zones by Imaris segmentation and masking. 5 cells per biological replicate were analyzed for a total of N = 3. ANOVA with post-hoc Tukey HSD: * = P < 0.05, ** = P < 0.01, and *** = P < 0.001. Error bars represent SEM. Scale bar: 10 microns (500 nm for zoomed ROIs). www.nature.com/scientificreports/ showed that ZIKV-induced ER reorganization from perinuclear stacked rough ER sheets to a perinuclear, circular region of convoluted smooth ER tubules (Fig. 4). These results further support the idea of ZIKV induction of a perinuclear tubular matrix.

ER localization of ZIKV NS2B and NS4B structural proteins. 3D STED analysis of cells labelled for
the ZIKV NS proteins NS2B and NS4B showed a predominant distribution of both to the CER and more precisely to the dense ZIKV-induced crescent-shaped tubular matrix in ERmoxGFP transfected U87 cells (Fig. 5A). While NS2B is predominantly associated with the dsRNA-positive CER, NS4B labeling extended throughout the CER as well as to the PER (Fig. 5A). To quantify this, we identified NS2B-positive and NS4B-positive ER segments and determined their overlap with total ER, CER and PER (Fig. 5A). While NS4B was present at high levels on both PER and CER segments, NS2B was enriched in the CER relative to the PER and presented a similar ER distribution to dsRNA (Fig. 5A). The majority (> 55%) of dsRNA-labeled puncta were associated with NS2B or NS4B, consistent with the presence of both these NS proteins in the ZIKV-induced tubular matrix and a role in viral replication. An increased number of NS2B spots (~ 25% of total) overlapped with dsRNA spots relative to NS4B spots (~ 10% of total) (Fig. 5B). These two ZIKV NS proteins therefore exhibit distinct distributions within the ZIKV-induced tubular matrix when not associated with dsRNA replication complexes (Fig. 5B). Together with the differential distribution of NS2B and NS4B within the ER as a whole (Fig. 5A), these results highlight that these two ZIKV NS proteins do not associate exclusively with replication factories and argue that following synthesis of the ZIKV polyprotein, NS2B and NS4B undergo distinct biosynthetic pathways before reuniting in ER-associated replication complexes. www.nature.com/scientificreports/ Deep learning identifies ZIKV-induced ER reorganization. Deep learning has been successfully applied to the task of image recognition, distinguishing the category or class to which a given image belongs.
Deep learning architectures have outperformed shallow architectures when benchmarked on ImageNet, a dataset widely used by the computer vision community that contains roughly 14 million images belonging to 20,000 classes (e.g., cats, dogs, plant species, various modes of transportation) [33][34][35] . Deep convolutional neural networks (CNNs) are capable of learning local and global spatial patterns from raw training data (i.e. pairs of images and corresponding classes) enabling inference of the correct class of unseen images, and have achieved state of the art performance when benchmarked on ImageNet 36 . Importantly, compared to non-deep learning detection methods such as machine learning, deep learning approaches: (1) avoid the need to design and select features www.nature.com/scientificreports/ for highly complex ER structures that are required by non-deep learning methods; and (2) provide the ability to move beyond simple classification to inspect discriminative regions (i.e. subregions of the ER within each cell). We therefore applied deep neural networks to identify and distinguish the morphological features of the ER of ZIKV-infected cells. A pipeline outlining our approach is shown in Fig. 6A. We train a CNN using 2D slices (each representing a single Z-slice) from 3D STED volumes of ERmoxGFP and Sec61β-GFP labeled ZIKV-and mock-infected cells. Our CNN builds off of VGG16, a deep neural network architecture proposed by Simonyan and Zisserman 35 which achieved state of the art results on the 2014 ImageNet Large Scale Visual Recognition Challenge (ILSVRC) 37,38 .
To improve computation time, we initially performed the analysis on downsampled STED images. As the 3D STED data sets were relatively small to train the CNN from scratch, we leveraged a network already pretrained on natural images. This transfer learning technique speeds up the convergence rate of networks when dealing with small target datasets. Certain filters (combinations of weights) learned on the first dataset (i.e. ImageNet) may still be useful for classifying a second dataset (i.e. STED); as a result, less weight updates are needed before achieving good performance. Using a pretrained VGG16 as our base model we obtained a 20% boost in test accuracy, compared against a random weight initialization. Using ERmoxGFP labelled ER alone, the CNN was able to distinguish between ZIKV-and mock-infected cells with an 82% accuracy (Fig. 6B, top left). Accuracy using Sec61β-GFP labelled ER was only 66%, whereas combining cells labeled with both ER markers, increasing the number of training samples seen by the network, still resulted in only 71% accuracy. Based on the improved classification performance, further analysis was performed using ERmoxGFP labelled cells.
With a limited sample size of 56 ERmoxGFP labelled cells, we applied leave-one-out cross validation; 56 cells were split into a set of 55 cells used to train the CNN, and a test set consisting of the remaining cell. This process is outlined in Fig. 6A, and is repeatedly applied using each cell as the test set. The 55 cells used to train the CNN are further split into a training and validation set, composed of 44 and 11 cells, respectively. During training, all 2D slices from cells belonging to the training set are passed individually to the CNN while 2D slices of cells belonging to the validation set are used to evaluate the performance and update the network's parameters (weights) accordingly. This grouping of cells into training, validation, and test sets ensures that information is not leaked, i.e., the CNN is only presented with unseen 2D slices during testing. ZIKV-and mock-infected class activation maps (CAMs) 39 are generated for each 2D slice to help define regions used by the CNN when inferring a class label (ZIKV or Mock). We achieved 82% accuracy and 82% specificity to identify ZIKV-and mock-infected cells based on a majority of slice predictions per cell. STED Z stacks extend beyond the cell and include Z slices that contain minimal ERmoxGFP signal and show poor accuracy to predict class label (Supp Fig 1). When considering those slices containing ERmoxGFP signal intensity greater than the median, we achieved 84% accuracy and 86% sensitivity. On a per slice basis, accuracy for all slices was 78% and sensitivity 79% that increased to 81% and 84%, respectively, when considering slices expressing ERmoxGFP greater than median intensity (Fig. 6B, Supp. Fig. 1). Considering the confusion matrix for all slices, cell label (Fig. 6B, top left), the 82% of ZIKV-infected 2D slices correctly predicted to be infected represent true positives (TP) while the 18% of ZIKV-infected 2D slices not predicted to be infected represent false negatives (FN). The 82% of mock-infected 2D slices correctly predicted to be uninfected represent true negatives (TN) and the 18% of mock-infected 2D slices predicted to be ZIKV-infected represent false positives (FP). The CNN has therefore accurately predicted which cells are ZIKVinfected based on the reorganization of ERmoxGFP-labeled ER. Improved accuracy for slices expressing higher ERmoxGFP highlights that the network is using the ERmoxGFP label to identify ZIKV-infected slices and cells.
To determine the basis for CNN decision making we analyzed the respective CAMs for both classes, ZIKVand mock-infected cells. This method, first proposed by Zhou et al. 39 produces a heatmap (or CAM), corresponding to a given image input, to localize the discriminating regions that the CNN uses to infer a class label (Fig. 7A). For instance, identification of ZIKV-infected slices (true positive; TP; orange) is primarily based on regions highlighted by the ZIKV CAM while identification of mock-infected cells (true negative; TN; green) is based on the mock CAM (Fig. 7B). The range of values of the generated CAMs is between 0 and 1. Class activation maps have been used for weakly supervised segmentation, i.e. to delineate regions of interest (ROIs) 40,41 We considered ROIs to be areas of the CAM with values greater than 80% of the maximum. Representative examples of ZIKV (red) and mock (blue) CAMs are shown in Fig. 7B; as they transition to yellow, concentric rings correspond to smaller CAMs and increasing thresholds from 10-95%. In general, the mock CAM was centered over the cell for correctly labelled uninfected cells (TN) cells as was the ZIKV CAM for correctly labeled ZIKVinfected cells (TP) cells (Fig. 7B). Figure 7B, right, provides examples of slices mislabelled by the network; in the case of ZIKV-infected cells slices that are incorrectly labelled as Mock (FN). Representative downsampled ERmoxGFP-labeled images from true positive (TP) instances show the 80% ZIKV CAM in red covering the central ER region. Higher resolution TP patches (112 × 112 pixels) from these regions of high network attention (ZIKV CAM with a threshold of 80%) show ER structures used to correctly label ZIKV-infected cells (Fig. 7C). Fig. 8A, ROIs identified by the ZIKV CAM show consistently higher ER density at all thresholds when comparing correctly identified infected slices (TP, orange) to the mock-infected slices (TN, green). Similarly, for 2D slices correctly identified as mock-infected (TN), the regions found by the mock CAM show increased ER density at all thresholds compared to the TP slices. This suggests that the CAMs used to identify both ZIKV-and mockinfected cells correspond to high ER density regions localized over the cell (see Fig. 7B) and that VGG16 is using differences in the ER label (ERmoxGFP) to identify slices as either ZIKV-or mock-infected. Consistently, the density profile of the ZIKV CAM for cells falsely identified to be infected (false positive; FP; blue) matched that of TP (orange) ZIKV-infected cells while that of the mock CAM for cells falsely identified to be uninfected (false negative; FN; purple) matched the mock TN (green) profile. ER density for the ZIKV and mock TP and www.nature.com/scientificreports/ Figure 8. Normalized CAM ER density across subgroups. ROIs are defined by various CAM thresholds. For a given ROI, ER density is defined as total ERmoxGFP intensity within the ROI divided by ERMoxGFP area inside the ROI. ER density for each ROI defined by the CAM is then normalized by the ER density of the whole cell. (A) ER density of ROIs defined by CAM thresholds from 10-95% with increments of 5% is compared across 4 subgroups: ZIKV-infected 2D slices correctly predicted to be infected (true positives); mock-infected 2D slices correctly predicted to be uninfected (true negatives); ZIKV-infected 2D slices incorrectly predicted to be uninfected (false negatives); mock-infected 2D slices incorrectly predicted to be infected (false positives). (B) ER densities of 80% CAMs ROIs compared across subgroups. (C) Euclidean distances between center of mass of 80% CAMs ROIs and weighted center of mass of ERmoxGFP signal. Statistical analysis was performed using ANOVA with post-hoc Tukey HSD: * = P < 0.05, ** = P < 0.01, and *** = P < 0.001. Unless otherwise specified directly in line with each plot, we found all other pairwise comparisons to be statistically significant with a P < 0.01 (**). This helped to improve the visualization of our data points. www.nature.com/scientificreports/ TN profiles at a CAM threshold of 80% show clearly that regions identified by the thresholded ZIKV CAM have significantly increased ER density when comparing TP to TN slices (Fig. 8B, left). Conversely, regions identified by the thresholded mock CAM have increased ER density for TN compared to TP cells (Fig. 8B, right). We then calculated the average Euclidean distance between the weighted center of mass of the CAMs and the weighted center of mass for all ERmoxGFP labelled pixels (Fig. 8C). The center of mass of the ZIKV CAM for TP and for the mock CAM for TN were closest to the ERmoxGFP center of mass. As for the ER density profile of the CAMs, Euclidean distance of FPs for the ZIKV CAM matched TP, and FNs for the mock CAM matched TNs. Together, these data argue that the CAMs used to identify ZIKV-and mock-infected cells are located in dense ER regions near the center of mass of the ERmoxGFP label.

Deep learning identifies ZIKV-infected cells based on the high density CER. As seen in
To determine the relationship of the ZIKV CAM to the dense tubular matrix region that we identified by 3D STED super-resolution microscopy (Fig. 3), we assessed CAM overlap with dsRNA, NS2B, NS4B and ERmoxGFP labeling. To do this, we calculated the normalized CAM intensity, defined as the intensity sum of a given channel in the ZIKV CAM normalized by the total area of that channel. The 80% ZIKV CAM shows a higher degree of overlap with the NS4B-positive ER relative to the ERmoxGFP-, dsRNA-and NS2B-labeled ER (Fig. 9A, left). Further, based on Euclidean distance analysis, the center of mass of the ZIKV CAM was closer to the NS4B center of mass (Fig. 9A, right). Relative to NS2B and dsRNA, NS4B density is present throughout the CER and increased in the PER (Fig. 5). Therefore, to further define the ER zone contributing to CNN discrimination of ZIKV-infected cells, we assessed overlap of the ZIKV CAM with the dsRNA-positive and dsRNA-negative CER as well as the PER, as defined previously (Fig. 3). The PER, determined based on Z-height of ER segments, includes the subnuclear ER and to ensure we were assessing strictly the PER we segmented out the nuclear region using a semi-automated annotation approach (See Materials and Methods). As seen in Fig. 9B, the ZIKV CAM shows increased overlap and proximity with the dsRNA-positive CER relative to either the dsRNA-negative CER or PER. The region closest to and most overlapping with the ZIKV CAM is the nuclear region, cropped from the PER channel. As the center of mass of the ZIKV CAM presents the shortest Euclidean distances to the dsRNApositive CER and the nucleus, this locates the ZIKV CAM center of mass to the perinuclear CER region between the dense dsRNA-positive tubular matrix and the nucleus. This is visualized in Fig. 9C, where increased grey scale density highlights the weighted overlap of the ZIKV CAM with these four ER regions. Increased ZIKV CAM overlap with the dsRNA-positive ER, the perinuclear zone of the dsRNA-negative CER as well as the nuclear zone adjacent to this region can be clearly observed. The ZIKV CAM encompasses the dsRNA, NS2B and NS4B-positive dense tubular matrix region that houses ZIKV replication factories, but also includes the adjacent perinuclear region enriched for NS4B but not NS2B. Individual patches do not encompass the totality of the CER (see Fig. 7C) and the precise nature of the features that the CNN uses to discriminate between ZIKV-and mock-infected cells remains to be determined. CAM localization analysis shows that the neural network uses the same CER region that we have observed to be modified upon ZIKV infection. Deep learning therefore has the ability to identify the ER morphological changes associated with ZIKV infection.

Discussion
ZIKV infection is characterized by re-organization of the ER to create replication factories and convoluted ER membranes involved in viral replication, whose ultrastructure has been elegantly characterized by EM 19,25 . Here, 3D STED super resolution microscopy reveals the formation of a novel, perinuclear, crescent-shaped region of dense tubular ER, or convoluted membranes, within the CER of ZIKV-infected cells. This region is enriched for dsRNA and ZIKV non-structural proteins NS2B and NS4B and therefore corresponds to the site of ZIKV replication factories and genomic replication 19,20 . 3D STED microscopy therefore shows that ZIKV infection induces CER reorganization. Importantly, based largely on the formation of this dense CER region, deep learning image analysis is able to identify ZIKV-infected cells based solely on ER morphology. This highlights the potential use of deep learning to further define ER morphology in terms of biologically relevant computational features 15,42 .
The ER is a morphologically complex organelle, containing smooth ER tubules and ribosome-studded rough ER sheets identified ultrastructurally by EM since over 60 years 43 . In confocal images of cultured cells, these morphological structures correspond, respectively, to PER tubules and the dense perinuclear CER 3 . Here, we approximate CER and PER regions based upon Z bounding box height of segmented 3D STED ER volumes. By use of our segmented ER volumes, the dense ER zone enriched for dsRNA is determined to be within the CER. The CER of both ZIKV-and mock-infected U87 cells is composed of tubular networks that correspond to previously described PER tubular matrices 7 . While we were unable to detect ER sheets by super-resolution analysis of cultured U87 cells, EM of brain organoids shows the transformation of ER sheets to convoluted membranes upon ZIKV infection. Overall, the fixed cell 3D STED analysis applied here suggests that structures associated with ZIKV replication, previously denoted as convoluted membranes, derive from tubular matrix reorganization in the CER.
The ZIKV-induced CER-localized, high ER density tubular matrices are enriched for dsRNA localizing active replication factories to this CER domain 19 . Functionally, the formation of tubular matrix is suggested to stabilize high curvature regions of the ER with increased membrane surface area that act as a reservoir for membrane proteins and lipids 7 . ZIKV-induced dense tubular matrix may serve to provide material required for formation of replication factories and viral replication. Previously, tubular matrices were shown to shaped by both RTN and ATL proteins 7 . Consistently, depletion of RTN3.1A, ATL2 and ATL3 reduce ER reorganization and replication of multiple flaviviruses, including ZIKV 25,45,46 . Here, NS2B and NS4B colocalization with dsRNA localizes these proteins to replication factories. Interestingly, dsRNA does not show complete co-localization to NS4B or NS2B. Future studies investigating the localization of NS proteins to replication factories using live cell imaging could further define the partial co-localization observed here. Furthermore, NS2B distribution is restricted to Scientific Reports | (2020) 10:20937 | https://doi.org/10.1038/s41598-020-77170-3 www.nature.com/scientificreports/ www.nature.com/scientificreports/ the ZIKV-induced tubular matrix, consistent with its recruitment of NS3 to the ER and role as a cofactor for NS3 protease function, necessary for ZIKV genomic replication and polyprotein cleavage 25,28,47 . In contrast, NS4B, is a multifunctional, membrane protein involved in many processes of ZIKV infection and pathogenesis, including promoting ER membrane proliferation, curvature of ER membranes to produce replication factories and disruption of ER-mitochondria contacts allowing the virus to successfully subvert the host innate immune response 26,48,49 . The broad distribution of NS4B throughout the ER suggests a role for this non-structural protein in ZIKV-induced ER reorganization. We show that deep convolutional neural networks can accurately classify ZIKV-and mock-infected cells based on ER labeling. Using CAMs, we focus on regions which are more likely to be used by the network to infer a class label. Importantly, our analysis shows that the VGG16 neural network uses dense ERmoxGFP in the CER as an indicator of ZIKV infection. This demonstrates that the neural network is detecting the ZIKV-induced reorganization of the CER that leads to the formation of the tubular matrix associated with ZIKV replication. Accuracy obtained for VGG16 classification of ZIKV-infected cells from 3D STED image stacks is comparable to prior classification using VGG16 37,38,[50][51][52] . CNN deep learning analysis was better able to distinguish ZIKVinfected from mock-infected cells based on ERmoxGFP compared to Sec61β-GFP labeling. ERmoxGFP was developed as an inert ER lumenal marker to only minimally impact that native structure and physiology of the ER 29 . Fluorescently-tagged ER markers such as the ER membrane marker Sec61greek beta-GFP can impact ER organization, including ER nanodomain localization 6,44 . Indeed, we did detect a significant change in CER/PER volumes using Sec61β as a reporter (Fig. 2C). The lumenal reporter, ERmoxGFP, therefore better detects ER morphology differences induced by ZIKV.
The interpretability of artificial intelligence is an evolving field and we believe that interpretable methods, such as Grad-CAM 53 , are important tools for the understanding of deep neural networks applied to exploratory data sets. This approach has now allowed us to identify features of discriminatory regions, and has not, to our knowledge, been applied to subcellular morphology, nor to 3D super resolution images of the ER. This highlights the sensitivity of deep learning image analysis and augurs well for future identification not only of tubular matrix but also of the peripheral tubules, 3-way junctions, sheets and organelle contacts that comprise the ER 54 . However, defining the morphological characteristics associated with deep learning-based decision making from isolated 2D patches of ER remains a challenge. Future application of similar deep learning approaches to superresolution 3D analysis of both fixed and live cells may lead to novel understanding of the dynamic, local changes to ER morphology associated with the diverse functions of this cellular organelle. To the best of our knowledge, interpreting the network's decisions is a novel approach to study ER morphology using deep learning.
Application of a deep learning framework to host-cell ER labeling correctly identified ZIKV-infected cells from 2D sections of downsampled STED super-resolution images; this indicates that VGG16 and other deep CNN architectures are capable of finding discriminatory features in lower resolution images for correct identification of infected cells. Importantly, our deep learning approach is translatable to ER reorganization induced by other pathogenic human enveloped viruses 15,42,55 . www.nature.com/scientificreports/ serum (ThermoFisher Scientific) and 1% bovine serum albumin (Sigma) in PBS-CM for 1 h; (5) incubated with primary antibodies in the same blocking solution as described above for 1 h at room temperature (RT) then washed quickly with PBC-CM once followed by three 5  3D STED microscopy. 3D STED imaging was performed with the 100X/1.4 Oil HC PL APO CS2 objective of a Leica TCS SP8 3X STED microscope (Leica, Germany) using white light laser excitation with HyD detectors, and Leica Application Suite X (LAS X) software. Sample acquisition was done at a scan speed of 600 Hz with a line average of 6. Time-gated fluorescence detection was used to further improve lateral resolution. GFP was excited at 488 nm and depleted using the 592 nm depletion laser. Alexa Fluor 532 was excited at 528 nm, and Alexa Fluor 568 was excited at 577 nm. Both were depleted using the 660 nm depletion laser. To avoid crosstalk, image stacks for each channel were acquired sequentially, in the order of Alexa Fluor 568, Alexa Fluor 532 and then GFP. A step size of 210 nm was utilized for 3D STED. Images were then deconvoluted using Huygens Professional software (Scientific Volume Imaging). To determine predicted 3D STED resolution, full width at half maximum (FWHM) values were obtained from theoretical point spread functions with Huygens Professional. Predicted resolution is 126 nm for XY and 347 nm for Z. 3D STED ER segmentation and ER region mask generation were done using Imaris × 64 9.2.1. dsRNA associated ER and NS4B/NS2B associated ER regions were generated by partial overlap analysis after surface generation. All data produced from Imaris were parsed and analyzed with a custom made Jupyter notebook. Two-tailed unpaired Student's t-tests or ANOVA with post-hoc Tukey HSD were done using Graph Pad Prism 6.0.

Antibodies
Transmission electron microscopy. 63-day old cerebral organoids were infected with the PRVABC59 ZIKV strain, MOI of 1, for 48 h. Organoids were then processed by replacing media with the primary fixative [0.1 M sodium cacodylate, 1.5% paraformaldehyde, 1.5% glutaraldehyde, pH 7.3, room temperature] 56 and allowed to sit for 1 h before placing them at 4 °C overnight. The organoids were washed three times (10 min each wash) with buffer (0.1 M sodium cacodylate) and then post-fixed for 1 h on ice in 0.1 M osmium tetroxide in 0.1 M sodium cacodylate (pH7.3). They were then washed three times with distilled water (10 min each wash, room temperature) and stained en bloc with 1% aqueous uranyl acetate for 1 h at room temperature. They were washed again three times with distilled water and then dehydrated through an ascending concentration series of ethyl alcohol (30%, 50%, 70%, 95%, 100%), passed through two changes of propylene oxide, and then infiltrated overnight in a 1:1 solution of propylene oxide:EMBED 812 resin. The organoids were then passed through two changes of 100% EMBED resin, embedded and the resin polymerized at 60 °C for 48 h. Thin sections were cut using a diamond knife and sections were collected on copper grids. The sections were then stained with uranyl acetate and lead citrate, and images collected using a FEI Tecnai G2 Spirit transmission electron microscope operated at 120 kV.
Image data processing for CNN analysis. A custom pipeline was created using MATLAB R2019a (The Mathworks, Inc., Natick, MA) and Python 3 (Python Software Foundation, Scotts Valley, CA) scripts to convert data stored in Leica proprietary format (.lif) to an accessible and lightweight format, which could easily feed as input to state of the art CNNs. Once an imaging session is complete, the biologist exports Leica STED data without need for further processing or file conversion. The pipeline takes Leica .lif files of 3D STED volumes of multiple cells as input and outputs a single NumPy array per cell. The pipeline relies on open-source Matlab package Bio-Formats, modified to handle our needs. This pipeline has been validated through manual inspection of cells with Imaris image rendering software.

Deep learning.
The Keras open-source neural network library was used to modify and train neural networks using Tensorflow backend. VGG16 was used as a base model to predict infected from non-infected 2D slices (removing the original fully-connected output layers). This architecture and pretrained weights (learned using the ImageNet dataset) are easily accessible using Keras. The pretrained weights can be loaded when instantiating the base model. On top of the VGG16 max pooling layer, we add (in order): a global average pooling layer, a dropout layer (with a dropout rate of 0.5), dense layer (using ReLu activation, output dimension 1024), a second dropout layer (0.5 dropout rate) followed by a final dense layer (with softmax activation) for making predictions. These additional layers helped prevent overfitting of the CNN to the training data. 56 individual CNNs (with identical architecture as outlined above) were trained using a softmax loss function and RMSprop optimizer. Each model was then used to predict the class labels (i.e. assign either 'Mock' or 'ZIKV' to each 2D slice) of a single test cell. The CNNs were fine-tuned for a total of 24 epochs, after which only the best weights (resulting in highest classification performance on the validation set) were then reloaded during test time. CAMs were generated following the procedure outlined by Zhou et al. 39 , where global average pooling is applied to the output of the final convolutional layer of VGG16 followed by a weighted sum, according to the specified class (e.g., ZIKV). Models were trained on 2 Nvidia GPUs (GeForce GTX TITAN X, each with 12 GB of VRAM; the host machine had 32 GB of RAM).
Regions of interest generated through the deep learning framework were analyzed using Python programming. Statistical analysis across channels and ER regions relied on several open-source Python packages. SciPy's statistical functions module, was used to perform Kruskal-Wallis non-parametric tests, statsmodules Python