Subnanometer-resolution Structural Determination of Hemagglutinin from Cryo-electron Tomography of Influenza Viruses

Cryo-electron tomography (cryoET) has been applied over the past decades to capture snapshots of protein complexes, viruses, cells, and organisms. A modality of cryo-electron microscopy (cryoEM), cryoET is a structural biology method where a biological sample is flash frozen, then imaged through a variety of orientations through tilting^1,2,3. Images taken at each orientation are then computationally aligned to their common tilt axis and reconstructed into a tomogram to provide a three-dimensional view⁴.

Whereas X-ray crystallography and single particle cryoEM demand purified, structurally homogenous molecules, cryoET can image a molecule directly within its native context⁴. Therefore, one main advantage of cryoET is its ability to visualize pleomorphic samples, such as membranous viruses, including influenza^5,6,7. Another promise of cryoET is its ability to image across scales. While tomograms are not typically resolved past 5-10 nm⁸, the integration of subtomogram averaging, where copies of the same particle are identified, aligned, and averaged, can result in near atomic resolution in some biological molecules such as ribosomes^9,10. However, only limited types of molecules can reach this resolution; subtomogram averages do not typically surpass 10-15 Å resolution. In contrast, single particle cryoEM routinely achieves resolutions of 3-4 Å post the resolution revolution¹¹. Recent advancements in both higher throughput of cryoET data acquisition and analysis software have allowed for subnanometer resolution structure determination of additional biological molecules within their native context^{12,13,14,15,16,17,18}.

One common usage for cryoET is to visualize virus morphology, organization, and structure. Despite the lower resolution afforded by this technique compared to single particle cryoEM or X-ray crystallography, cryoET combined with subtomogram averaging can provide information on how viral proteins behave in their native environment and provide crucial details on their organization in the context of the virion. A common target for cryoET of viruses is the surface glycoproteins that are commonly used for host cell attachment and fusion, as they are often the main antigens and targets for therapeutics or vaccines. With recent advancements in cryoET processing packages, it has become increasingly feasible to achieve subnanometer resolution averages of these glycoproteins^19,20,21,22. One such example is hemagglutinin (HA), the major protein on the surface of influenza virions. Not only does this protein conduct both receptor binding and membrane fusion, but it also covers the virion in an incredibly dense manner, with hundreds to thousands of HAs on a singular virion⁵. The protocol presented here (Figure 1) integrates several commonly used packages with in-house scripts to delineate stages from pre-processing to model refinement for a subtomogram average of influenza HA.

NOTE: Sample datasets used for this protocol can be accessed at EMPIAR-12864, which includes the two sets of tilt series used for this protocol. The tilt series are collected from manually plunged grids of purified influenza A virus at a physical pixel size of 2.09 Å/pixel to ensure a large enough field of view so that each tilt series contains several virions, and also to render the highest resolution reconstructions possible. For users' own datasets, it is recommended to start the workflow with raw tilt movies. These datasets were processed and visualized using high-performance workstations. The Table of Materials lists the hardware and software used for this protocol. All software packages used in this protocol are open source and available for download; installation links and instructions are listed in the Table of Materials. The recommended workstation for processing cryoET datasets should be equipped with at least an 8-core processor, a dedicated GPU card with 6 GB of VRAM, 64 GB of RAM, and 2 TB of local storage.

1. Data preprocessing of tilt movies and reconstruction of cryo-electron tomograms in Warp ²³ and IMOD ²⁴

1. In a terminal, activate a conda environment that has Warp installed by using the following command:
   conda activate warp_environment
2. Create a frame series settings file for processing frame series. This is a config file that contains metadata about the microscope, the location of data to process, and where to store output files.
   WarpTools create_settings --folder_data path/to/.tif --folder_processing warp_frameseries --output warp_frameseries.settings --extension “*.tif” --angpix 1.04 --gain_path gain_file.mrc --exposure 3.07
   NOTE: These data were collected at super-resolution mode.
3. Conduct contrast transfer function (CTF) estimation and motion correction.
   WarpTools fs_motion_and_ctf --settings warp_frameseries.settings --m_grid 1x1x5 --c_grid 2x2x1 --c_range_max 7 --c_defocus_max 10 --c_defocus_min 4 --c_use_sum --out_averages
   NOTE: m_grid parameter should correspond to the number of subframes within a tilt movie - our datasets consisted of 5 subframes per tilt movie.
4. Import tilt series metadata so that WarpTools can identify which movies belong to which tilt series.
   WarpTools ts_import --mdocs path/to/.mdoc --frameseries /path/to/frameseries --tilt_exposure 3.07 --min_intensity 0.3 --output tomostar
5. After the tomostar folder has been populated, create a tilt series settings file for tilt series processing
   WarpTools create_settings --folder_data tomostar --folder_processing warp_tiltseries --output warp_tiltseries.settings --extension “*.tomostar” --angpix 1.04 --gain_path gain_file.mrc --exposure 3.07 --tomo_dimensions NxNxN
6. Write out tilt stacks and conduct automated fiducial-based tilt-series alignment using IMOD's batchruntomo program using the Etomo GUI²⁴.
   1. Run the following command in the terminal:
      Warptools ts_stack --settings warp_tiltseries.settings --angpix 8.35
      NOTE: The tiltseries were exported at 4x the physical pixel size to decrease alignment time and computational resources.
   2. Choose a sample dataset to first set alignment parameters before applying them as a template for batchruntomo.
   3. Import alignment parameters from IMOD if using batchruntomo.
      WarpTools ts_import_alignments --settings warp_tiltseries.settings --alignments warp_tiltseries/tiltstack/ --alignment_angpix 8.35
      NOTE: For this dataset, using the Etomo GUI produced better results than the wrappers within WarpTools.
   4. Alternatively, conduct automated fiducial-less tiltseries alignment using the AreTomo²⁵ wrapper within WarpTools.
      WarpTools ts_aretomo --settings warp_tiltseries.settings --angpix 8.35 --alignz 1000 --axis_iter 3 --exe AreTomo_executive
      NOTE: The AreTomo wrapper was tested with AreTomo1.3.4.
7. Run the following command in the terminal to estimate CTF parameters for tilt series:
   WarpTools ts_ctf --settings warp_tiltseries.settings --range_high 7 --defocus_min 2 --defocus_max 10 --auto_hand 4
8. Reconstruct tomograms using WarpTools.
   1. Set environment variable:
      export WARP_FORCE_MRC_FLOAT32=1
      NOTE: This step ensures tomograms will be compatible with subsequent stages of image processing and visualization in other softwares used in this protocol.
   2. To reconstruct tomograms, execute in a terminal: WarpTools ts_reconstruct --settings  warp_tiltseries.settings --input_data input file names --angpix 8.35 --dont_invert
9. Repeat steps 1.2 - 1.8.2 for all datasets.

2. Tomogram preprocessing and particle picking

1. In a terminal, activate a conda environment with IsoNet²⁶ installed.
   conda activate isonet_env
2. Prepare for tomograms processing by first creating a subfolder and moving all tomograms into that folder.
   mkdir tomo_folder
mv tomograms*.mrc tomo_folder/
3. Generate a star file in the project folder.
   isonet.py prepare_star tomo_folder --output_star tomograms.star --pixel_size 8.35
4. Using a text editor, open the generated star file and enter the approximate defocus value for the 0-degree tilt images in the fourth column (_rlnDefocus) for each tomogram, which can be located in the processed_items.json file in the warp_tiltseries folder.
   NOTE: The defocus value should be in angstroms for IsoNet. Warp writes out defocus values in µm, which requires a multiplication factor of 10,000.
5. Run the CTF deconvolve command in the terminal.
   isonet.py deconv tomograms.star --snrfalloff 0.7 --deconv_folder deconvolve
6. After deconvolving, initiate EMAN2 GUI²⁷ to begin preprocessing tomograms for particle picking.
   conda activate eman_env
e2projectmanager.py
7. Under Tomography, left-click on the arrow next to Raw Data and select Import Tomograms from the drop-down menu.
8. Left-click on the arrow adjacent to Segmentation and select Preprocess tomograms.
   NOTE: The default parameters are likely suitable for many datasets. For these tomograms, a lowpass filter at 4 Å was applied and tilt images were normalized. After preprocessing, an info directory will be automatically generated by EMAN2 with blank .json files (with similar base names to pre-processed files). The angpix needs to be added to the json files before particle picking.
9. Train convolutional neural network (CNN) to recognize HA glycoprotein.
   1. Change working directory to the location of the tomograms to be used for CNN training and particle picking:
      cd path/to/tomograms
   2. Use the terminal to open the CNN training window.
      e2spt_boxer_convnet.py --label label_name
   3. Confirm that four windows are open: one containing the information on CNN parameters and tomograms in the directory, and the other three windows containing images of good references (Positive), bad references (Negative), and particles selected by the CNN (Particles).
   4. Left-click the New button on the CNN GUI to initialize a new CNN. Set default learning rate to 0.0001 and box size to 8.
   5. Under the FileName column, left-click on a representative tomogram to open it in a new window.
      NOTE: Only one tomogram can be open at a time.
      1. To move through the z-axis of a tomogram, place the cursor over the open tomogram and press the scroll wheel of the computer mouse to open a new window. The slider labelled N# will change the z-axis.
   6. On the open tomogram, select 10-20 features corresponding to HA using the left mouse button in the Positive panel.
      NOTE: Select both top and side views of HA, which either appear as cylinders or triangles above the membrane, as good references for a more robust network.
      1. Images of positive references will appear in the Positive window. To remove a reference, hold the Shift key and left-click on a reference image.
   7. Switch to Negative panel and select 10-20 references corresponding to features such as empty membrane patches, vRNP density, fiducial markers, and debris.
      NOTE: References will appear in the Negative window.
   8. Initiate CNN training by left-clicking on the Train button in the main window.
      NOTE: The number of iterations (Niter) in the main window changes how many iterations of training the CNN undergoes. Recommend setting the parameter to 50.
   9. Once the CNN has been trained on the selected references, left-click Apply to use the CNN to pick particles in the open tomogram.
      NOTE: Images of picked particles can be seen in the Particles window. Selected particles will also be displayed on the tomogram in blue circles.
   10. Continue to select positive and negative references based on the applied network, periodically saving progress through the Save button.
   11. Once satisfactory, select Apply All to have the CNN select particles in all tomograms in the directory and evaluate the results on several tomograms; perform additional training as needed.
10. Save coordinates to a text file for each tomogram.
    1. Open the EMAN2 GUI using the e2projectmanager.py command.
    2. Click on the arrow next to Subtomogram Average and click Manual Boxing.
    3. Input the name of a tomogram and click Launch.
    4. Save the coordinates to the text file by selecting File > Save Box Coord > tomogram_ha.txt.
    5. Navigate to the neuralnets subfolder and info subfolder to back up nnet_save.hdf, trainouts.hdf, segouts.hdf, and boxes3dref.hdf.
11. To ensure analysis is only performed on fully assembled virions where the presence of the matrix (M1) protein layer and viral ribonucleoprotein (vRNP) complex is apparent, train a second convolutional neural network to recognize M1 protein.
    1. Follow the same protocol as outlined in step 2.9, changing the training box size from 8 to 14.

3. Particle curation

1. Download notebooks from https://github.com/jqyhuang/influenza-analysis.
2. Open CNN_Particle_Cleaning.ipynb notebook and load the required modules.
   NOTE: The script uses Open3D²⁸ as a package to visualize particle coordinates as 3D point clouds.
3. Load in text files corresponding to HA and M1 coordinates and view as 3D point clouds using the Open3D package.
4. Filter out HA coordinate outliers using statistical outlier removal as implemented in Open3D.
   NOTE: Suggested initial values for HA are nb_neighbors = 50 (number of neighbouring coordinates around a coordinate) and std_ratio = 0.5.
5. Calculate the distance between the HA and M1 point clouds. All HA coordinates more than 20 pixels away from the M1 point cloud will then be identified as an outlier. All particle coordinates not identified as an outlier will be saved in an output .txt file.
   NOTE: 20 pixels correspond to an approximate distance of 16 nm with a pixel size of 8.35 Å/pixel. The center of an HA trimer to the M1 layer is approximately 15 nm in our tomograms.
6. Concatenate all particles and save as starfiles using pts2starfile.ipynb.

4. Iterative subtomogram averaging and classification

1. Using WarpTools, extract particles at a binning factor of 4 and conduct initial rounds of subtomogram averaging in RELION4²⁹.
   WarpTools ts_export_particles --settings warp_tiltseries.setting --input_star pts2star.star --coords_angpix 8.35 --output_star bin4_export.star --output_angpix 8.35 --box 48 --diameter 140 --3d
   NOTE: The suggested box size of 48 x 48 x 48 pixels, corresponding to a ~360 ųbox that would be large enough to contain an array of 7-8 HA.
   1. Convert warp starfile to RELION 4 compatible file
      relion_convert_star --i bin4_export.star --o bin4_conv.star
   2. Generate an initial reference using relion_refine_mpi on a subset of particles.
      head -n 30 bin4_conv.star >> subset.star & tail -n +31 bin4_conv.star | shuf -n 2000 >> subset.star
mpiexec -n 3 relion_refine_mpi --o init_ref/job001/run --auto_refine --split_random_halves --i subset.star --firstiter_cc --ini_high 20 --dont_combine_weights_via_disc --pool 3 --pad 2 --ctf --particle_diameter 300 --flatten_solvent --zero_mask --oversampling 1 --healpix_order 2 --auto_local_healpix_order 4 --offset_range 14 --offset_step 4 --sym C1 --low_resol_join_halves 40 --norm --scale --j 12 --gpu 0:1 --pipeline_control init_ref/job001
   3. Use relion_refine_mpi to conduct 3D autorefine on the subtomograms.
      1. Sample command: mpiexec -n 3 relion_refine_mpi --o Refine3D/job001/run --auto_refine --split_random_halves --i bin4_conv.star --ref init_ref/job001/run_class001.mrc --firstiter_cc --ini_high 20 --dont_combine_weights_via_disc --pool 3 --pad 2 --ctf --particle_diameter 400 --flatten_solvent --zero_mask --oversampling 1 --healpix_order 2 --auto_local_healpix_order 4 --offset_range 16 --offset_step 4 --sym C1 --low_resol_join_halves 40 --norm --scale --j 12 --gpu 0:1 --pipeline_control Refine3D/job001
         NOTE: The initial global and local sampling was set to 7.5 and 1.8 deg, corresponding to a healpix order of 4 and local healpix of 2. The initial refinements mostly leverage the strong density of the membrane, M1 protein, and the HA array. The translation offset range is allowed to be larger to encompass any shifts that may occur when aligning the array.
   4. Repeat refinement after the first refinement run converges, changing the translation to 8 and step to 2.
2. Leverage 2D classification to discard junk particles using relion_refine.
   1. Sample command: relion_refine --o Class2D/job003/run --grad --class_inactivity_threshold 0.1 --grad_write_iter 200 --iter 200 --i Refine3D/job002/run_data.star --dont_combine_weights_via_disc --pool 3 --pad 2 --ctf --tau2_fudge 2 --particle_diameter 300 --K 20 --flatten_solvent --zero_mask -- strict_highres_exp 14 --center_classes --oversampling 1 --norm --scale --j 24 --skip_align --pipeline_control Class2D/job002.
      NOTE:2D classification was used instead of 3D classification as it is more computationally efficient. Subtomograms can be directly used as input to the 2D classification job without additional image processing.
   2. Combine classes containing an identifiable HA array containing cylindrical HA, membrane, and M1 density using relion_star_handler.
   3. First, create star files containing good classes.
      relion_star_handler --i input_file2.star --o output_good_class.star --select rlnClassNumber --minval goodclassnumber -maxval goodclassnumber
   4. Then, use the following relion_star_handler command to combine the individual star files.
      relion_star_handler --i "output_good_class1.star output_good_class2.star … output_good_classn.star" --o bin4_keep.star --combine
3. Re-extract at bin2 and unbinned particles for iterative local refinement.
   1. For bin2 refinement, extract particles using a smaller box size of 80 and conduct local search with both healpix and local healpix set to 4 (1.8 deg local angular sampling).
   2. At this stage, combine particle sets from different datasets using relion_star_handler.
      relion_star_handler --i input_file.star --o output_file.star --combine
   3. After the first round of refinement, create a cylindrical mask using e2filtertool.py within the EMAN2 environment that encompasses the central HA plus membrane and M1 density.
      NOTE: Parameters used for the cylindrical mask: cx=24, cy=24, outer_radius=8, zmax=40, zmin=12; all other parameters are unchecked.
   4. Conduct one additional round of refinement with the mask applied.
   5. Conduct another round of 2D classification to eliminate outliers that do not align with the central HA and additional debris.
   6. Repeat process with unbinned particles with a box size of 120 and create a soft mask that covers the central HA, align reference to C3 symmetry, and apply symmetry at this refinement stage.
      relion_image_handler --i bin1_ref.mrc --o bin1_c3.mrc --sym c3
      NOTE: The final resolution achieved using RELION was 6.1 Å.
4. Final refinement in MTools⁹
   1. Create refinement population (after activation of Warp conda environment).
      MTools create_population --directory refine_m --name ha_final
   2. Add data sources for each dataset.
      MTools create_source --name source_1 --population refine_m/ha_final.population --processing_settings warp_tiltseries.settings
      1. Repeat the previous step with additional datasets.
   3. Create the refinement species.
      MTools create_species --population refine_m/ha_final.population --name ha_todaysdate --diameter 160 --sym c3 --temporal_samples 1 --half1 last_relion_refine/run_half1_class001_unfil.mrc --half2 last_relion_refine/run_half2_class001_unfil.mrc --particles_relion last_relion_refine/run_data.star --mask mask.mrc
   4. Multi-particle refinement.
      1. First refine particle poses with the command: MCore --population refine_m/ha_final.population --refine_particles
      2. Subsequently, refine both particle poses and spherical aberration with the command: MCore --population refine_m/ha_final.population --refine_particles --ctf_cs
         NOTE: The refinement was stopped here as further iterations did not improve resolution. Different combinations of parameters that were not exhaustively tested, however, may continue to boost results.

5. Model refinement

1. Using ChimeraX³⁰, load in the final map and an atomic model of HA.
2. Segment map and combine all segments that correspond to the HA ectodomain, and use the Fit to Segments tool to dock in HA model.
   1. Save transformed coordinates of model.
3. Open the Phenix GUI³¹, and use the Real Space Refinement tool.
   1. When prompted to add files, add the transformed HA model and the map, and fill in the resolution of the final reconstruction.
   2. Conduct five rounds of global minimization, local rotamer fitting, occupancy refinement, and group ADP refinement.
4. Visualize final reconstruction and model using ChimeraX.

To demonstrate the utilization of this processing protocol (Figure 1), the previously outlined workflow was applied to two datasets of 25 tomograms combined, obtained from an H1N1 influenza A virus strain (A/Puerto Rico/8/1934). Data collection parameters are outlined in Table 1. Figure 2 illustrates a representative tomogram and zoomed-in views of pleomorphic influenza virions. Varied morphologies are captured in this tomogram, as virions range from spherical to oval/elongated in shape. While most influenza particles contain well-organized M1 and vRNP assemblies, some virions appear to be more disorganized and lack crucial structural components.

From this dataset, an initial set of 40,995 subtomograms was used for bin4 (8.35 Å/pix) reconstruction after particle picking and curation. The two datasets were initially processed independently in RELION4 with C1 symmetry and a wide spherical mask that encompassed the overall HA array. Three refinement cycles were carried out for these subtomograms, followed by 2D classification. Post-classification, poorly resolved subtomograms and junk particles were discarded; remaining subtomograms were extracted at bin2 (4.17 Å/pix) and the two datasets were combined. Bin2 subtomograms were first aligned together in a round of subtomogram averaging, then a cylindrical mask around the central HA was applied to focus alignment. At this stage, trimeric symmetry can be clearly visualized for the HA reconstruction. Additional rounds of refinement were carried out at bin2 and at 2.8 Å/pix. One last round of 2D classification was conducted with a small mask covering just the central HA trimer with aligned subtomograms. The major class, consisting of ~94% of remaining subtomograms, was extracted to unbinned particles and subjected to 3D refinement with C3 symmetry applied. Lastly, these subtomograms were exported to M, where particle poses and spherical aberration refinement cycles were carried out (Supplemental Figure 1).

The final subtomogram average (Figure 3A), consisting of 15,970 HA particles, reached a global resolution of 6.0 Å and a local resolution range of 5-7 Å (Figure 4). A model of PR8 HA was flexibly refined into the density; at FSC=0.5 and FSC=0.143, the map-to-model resolution was 8.1 Å and 6.6 Å, respectively. The architecture of the HA reconstruction greatly resembled previous cryoEM and cryoET maps. At this resolution, alpha-helices and beta sheets can be distinguished (Figure 3B); moreover, glycans can begin to be identified at four glycosylation sites on the HA head and stem (Figure 3C).

Our results show the suitability of cryoET in the reconstruction of HA from native influenza virions. Through the protocol, clear density for the cylindrical glycoprotein was observed perpendicular to the viral membrane, and the resolution improved through each step. For one's own data, it is suggested that the subtomogram averaging should start from a binned particle stack, and the results should be closely monitored at each refinement cycle. If glycoprotein density is not apparent from the beginning stages, it is recommended to map the subtomogram locations back onto the tomogram to verify accurate positioning. Otherwise, one can tweak alignment parameters or implement additional classification stages to achieve optimal results.

Figure 1: Overall pipeline for subtomogram averaging of HA from cryoET of influenza virus. The top panel represents a summary workflow for the two datasets used to demonstrate the protocol. The second panel corresponds to Section 1 of the protocol, the third panel corresponds to Sections 2-3, and the fourth panel corresponds to Sections 4-5. Please click here to view a larger version of this figure.

Figure 2: Representative tomogram of PR8 influenza virus. (A) Slice through reconstructed tomogram. Scale bar is 100 nm. (B-D) Zoomed in view of (B) spherical, (C) cylindrical, (D) M1-less PR8 virion. All scale bars in B-D correspond to 50 nm. Please click here to view a larger version of this figure.

Figure 3: Subtomogram average of PR8 HA. (A) Top and side view of the HA reconstruction at two contour levels flexibly fitted with a PR8 HA single particle cryoEM structure. (B) Clipped views through the HA reconstruction. Colored arrows correspond to colored boxes. (C) HA reconstruction shown at lower contour to unveil glycan density. Close-up views of glycans in stick representation in the cryoET map are also shown. Please click here to view a larger version of this figure.

Figure 4: Resolution estimate of HA subtomogram average. (A) Local resolution estimate of the HA subtomogram average mapped onto the reconstruction. Color bar depicts 5-7 Å on the blue-white-red palette. (B) FSC curves of half maps and of map-to-model resolution. The blue curve is of the unmasked reconstruction, and the red is of the masked. Please click here to view a larger version of this figure.

Table 1: Data collection parameters for cryoET datasets of PR8 influenza virus.

Supplemental Figure 1: Subtomogram averaging workflow for HA. Iterative alignment, averaging, and classification are conducted for HA subtomograms through gradual unbinning. Please click here to download this File.

Better structural understanding of critical viral proteins can accelerate the discovery of novel treatments against these viruses. In the past decade, the resolution revolution has accelerated the determination of high-resolution viral structures using single particle cryoEM, but this method is limited to either purified proteins or non-enveloped viruses with icosahedral symmetry. In contrast, cryoET is capable of visualizing morphologically diverse membranous virions in 3D, but is of limited resolution. The pipeline presented here builds on the widely used Warp-RELION-M pipeline^9,32 and integrates a set of custom scripts with the semi-automated convolutional neural network-based particle picker from EMAN2, in order to establish an end-to-end protocol and generate a subnanometer resolution reconstruction of influenza HA.

One of the more difficult problems in cryoET is choosing the appropriate suite of analysis tools after data collection, as numerous software programs exist for each step of tomographic analysis, from initial preprocessing to model building. The pipeline was adapted from the Warp-RELION-M^9,32 pipeline as it is both well-documented, with several established protocols contributed by the developers and the community, and has also produced numerous high-resolution subtomogram averages in the past years, including our own work^19,32,33,34. Additionally, this pipeline limits the transfer of metadata between software packages during the various processing stages, as different packages use varied conventions that may result in additional errors. Warp and M can connect the initial preprocessing stages with the final multi-particle refinement, which ensures the initial calculations of motion and CTF parameters are integrated with final refinement and map reconstruction. To this end, RELION4 was used for the majority of the subtomogram averaging stages, as Warp can generate particle stacks and metadata for refinement and classification in RELION. Post subtomogram refinement, results from RELION can be directly imported into M for multi-particle refinement. Lastly, while our particle alignment workflow was relatively straightforward, RELION is capable of more extensive classification and flexible refinement protocols if needed.

A particularly challenging aspect of influenza structural biology is the dense layer of glycoproteins that covers the surface of extensively pleomorphic virions. Compared to common cryoET targets such as the ribosome, they are too small to be accurately identified using 3D template matching. Moreover, hundreds to thousands of HA glycoproteins cover the viral surface; therefore, manual identification of particles of interest is time consuming and tedious. This aspect of influenza is unlike other viruses of a similar size, like SARS-CoV-2⁶ or HIV-1^7,35 that only have a few dozen glycoproteins on each particle. Lastly, as these virions are not uniformly spherical or oval but of mixed morphology, the oversampling approach (most commonly implemented using the Dynamo package^36,37,38) that is commonly employed for evenly shaped virions or viral-like particles may require additional effort to create a surface model for each virion. This aspect is especially relevant when working with large datasets or when combining datasets, as the amount of time spent on particle picking would scale with the number of tomograms analyzed. Of course, if working with a virion that is more uniform in morphology with a known radius, the oversampling approach can be more efficiently applied compared to this influenza dataset. To simplify manual input and decrease the time spent on particle picking, the CNN-mediated particle picking approach employed in this protocol aims to recognize the shape of HA in its local environment and can be applied to entire datasets of tomograms with limited training. Only several dozen particle images for both 'good' and 'bad' references are required to train a neural network; this process typically takes only a few hours for a user familiar with HA. Using this approach, an initial dataset of more than 40,000 subtomograms was curated. While not an absolute metric for final resolution, obtaining a large particle set is required if the goal is to obtain a subnanometer resolution reconstruction.

Another metric to enhance resolution is particle set curation, which was implemented in this protocol through initial coordinate-based filtering as well as subsequent rounds of classification. For influenza virions, not only was a neural net trained to recognize HA, the protein of interest, but a separate neural net was also trained to recognize the M1 protein layer beneath the membrane. This step is to ensure that the HA subtomograms that end up being integrated into our dataset are from intact virions. To further enhance the resolution of our final reconstruction, several 2D classification steps were included in this protocol, which is a computationally efficient way to parse and classify several tens of thousands of subtomograms. If one singular dataset does not generate adequate numbers of subtomograms, particles derived from several collections can be combined. To best ensure the compatibility of different datasets, particles were combined after preliminary rounds of subtomogram alignment and averaging. Each dataset was first pre-processed independently through Warp, as collection parameters across sessions can be different. However, these metrics are considered during final particle alignment in M, when final iterations of particle alignment are conducted to optimize for image deformation and electron-optical aberration parameters.

While our proposed protocol greatly improved subtomogram averaging resolution for HA, the Nyquist limit of ~4.2 Å was not reached. This limitation might be due to the final number of subtomograms incorporated into our reconstruction. While a dataset of ~15,000 subtomograms was curated, it remains possible that increasing the number of total tilt series would further improve the quality and resolution of the HA subtomogram average. Moreover, these data were acquired at a high defocus of 4-8 µm to ensure tomographic contrast. Tilt series acquisition at a lower defocus range would better retain high-frequency information that could improve the resolution of subtomogram averages.

This protocol shows an example of using cryoET and integrated analysis to determine a subtomogram averaging of influenza.