Systematic multi-trait AAV capsid engineering for efficient gene delivery

Ethical StatementAll mouse procedures were performed as approved by the Broad Institute Institutional Animal Care and Use Committee (IACUC), approval number 0213-06-18-1. For the cynomolgus macaque experiment (n = 1), the study plan involving the care and use of animals was reviewed and approved by the Charles River CR-LAV Institutional Animal Care and Use Committee (IACUC). The rhesus macaque study (n = 2) was conducted at the NIH Nonhuman Primate Testing Center for Evaluation of Somatic Cell Genome Editing Tools at the University of California, Davis (UC Davis). All procedures conformed to the requirements of the Animal Welfare Act, and protocols were approved prior to implementation by the UC Davis IACUC.AnimalsAll mouse procedures were performed as approved by the Broad Institute Institutional Animal Care and Use Committee (IACUC), approval number 0213-06-18-1. Unless otherwise stated, all mice were female C57BL/6J (000664) mice from the Jackson Laboratory (JAX). During individual variant characterization, recombinant AAV vectors were administered intravenously via the retro-orbital sinus in young adult female (7–8-week-old) C57BL/6J animals (n = 5 mice per vector group). Mice were randomly assigned to groups based on predetermined sample sizes. No mice were excluded from the analyses. For all assays, mice were euthanized with EUTHASOL™ (Virbac) and transcardially perfused with phosphate buffer saline, pH 7.4, at room temperature (RT). Experimenters were not blinded to the sample groups.For the cynomolgus macaque experiment (n = 1), the study plan involving the care and use of animals was reviewed and approved by the Charles River CR-LAV Institutional Animal Care and Use Committee (IACUC). During the study, the care and use of animals were conducted by CR-LAV with guidance from the USA National Research Council and the Canadian Council on Animal Care (CCAC). The Test Facility is accredited by the CCAC and AAALAC. Per the CCAC guidelines, this study was considered as a category of invasiveness C.The rhesus macaque study (n = 2) was conducted at the NIH Nonhuman Primate Testing Center for Evaluation of Somatic Cell Genome Editing Tools at the University of California, Davis (UC Davis). All procedures conformed to the requirements of the Animal Welfare Act, and protocols were approved prior to implementation by the UC Davis IACUC.Modeling and assessment library designThe modeling and assessment libraries were designed to contain 150K nucleotide sequences each. The libraries were composed of 64.5K unique and 10K shared amino acid sequences generated by uniformly sampling all 20 amino acids at each position. The combined 74.5K variants were duplicated via codon replication. 1K nucleotide sequences containing stop codons were included to detect any potential problems with cross-packaging.Capsid library synthesisTo produce synthetic library inserts, lyophilized DNA oligonucleotide libraries (Agilent G7223A) or NNK hand-mixed primers (IDT) were spun down at 8000 RCF for 1 min, resuspended in 10 μL UltraPure DNase/RNase-Free Distilled Water (Thermo Fisher Scientific, 10977015), and incubated at 37 °C for 20 min. For pooled synthetic oligonucleotide libraries, the following primer format was used: 5′-GTATTCCTTGGTTTTGAACCCAACCGGTCTGCGCCTGTGC-(NNN)7-TTGGGCACTCTGGTGGTTTGTGGCCAC-3′. To produce NNK inserts, the AAV9_K449R_Forward (5′-CGGACTCAGACTATCAGCTCCC-3′) and AAV9_K449R_NNK_Reverse (5′-GTATTCCTTGGTTTTGAACCCAACCGGTCTGCGCCTGTGC-(MNN)7-TTGGGCACTCTGGTGGTTTGTG-3′) primers were used.To amplify the oligonucleotide libraries and incorporate them into an AAV9 (K449R) template, 2 μL of the resuspended pooled oligonucleotide library or NNK-based library was used as an initial reverse primer along with 0.5 µM AAV9_K449R_Forward primer in a 25 μL PCR amplification reaction using Q5 Hot Start High-Fidelity 2× Master Mix (NEB, M0494S). 50 ng of a plasmid containing only AAV9 (K449R) VP1 amino acids 347–586 was used as a PCR template. PCR was performed following the manufacturer’s protocol with an annealing temperature of 65 °C for 20 s and an extension time of 90 s. After six PCR cycles, 0.5 µM AAV9_K449R_Reverse (5′-GTATTCCTTGGTTTTGAACCCAACCG-3′) was spiked into the reaction as a reverse primer to further amplify sequences containing the oligonucleotide library for an additional 25 cycles. To remove the PCR template, 1 μL of DpnI (NEB, R0176S) was added to the PCR reaction and incubated at 37 °C for 1 h. Afterward, the PCR products were cleaned using AMPure XP beads (Beckman, A63881) following the manufacturer’s protocol.The PCR insert was assembled into 1600 ng of a linearized mRNA selection vector (AAV9-CMV-Express) with NEBuilder HiFi DNA Assembly Master Mix (NEB, E2621L) at a 3:1 insert:vector Molar ratio in an 80 μL reaction volume, incubated at 50 °C for 1 h, and then at 72 °C for 5 min. Afterward, 4 μL of Quick CIP (NEB, M0508S) was spiked into the reaction and incubated at 37 °C for 30 min to dephosphorylate unincorporated dNTPs that may inhibit downstream processes. Finally, 4 μL of T5 Exonuclease (NEB M0663S) was added to the reaction and incubated at 37 °C for 30 min to remove unassembled products. The final assembled products were cleaned using AMPure XP beads (Beckman, A63881) following the manufacturer’s protocol and their concentrations were quantified with a Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, Q32851) and a Qubit fluorometer.mRNA selection vectorThe mRNA selection vector (AAV9-CMV-Express) was designed to enrich for functional AAV capsid sequences by recovering capsid mRNA from transduced cells. AAV9-CMV-Express uses a ubiquitous CMV enhancer and AAV5 p41 gene regulatory elements to drive AAV Cap expression. The AAV-Express plasmid was constructed by cloning the following elements into an AAV genome plasmid in the following order: a cytomegalovirus (CMV) enhancer-promoter, a synthetic intron, and the AAV5 P41 promoter along with the 3′ end of the AAV2 Rep gene, which includes the splice donor sequences for the capsid RNA. The capsid gene splice donor sequence in AAV2 Rep was modified from a non-consensus donor sequence CAGGTACCA to a consensus donor sequence CAGGTAAGT. The AAV9 capsid gene sequence was synthesized with nucleotide changes at S448 (TCA to TCT, silent mutation), K449R (AAG to AGA), and G594 (GGC to GGT, silent mutation) to introduce restriction enzyme recognition sites for oligonucleotide library fragment cloning. The AAV2 polyadenylation sequence was replaced with a simian virus 40 (SV40) late polyadenylation signal to terminate the capsid RNA transcript.AAV productionFor library production, HEK293T/17 cells (ATCC, CRL-11268) were seeded at 22 million cells per 15 cm plate the day before transfection and grown in DMEM with GlutaMAX (Gibco, 10569010) supplemented with 5% FBS and 1× non-essential amino acid solution (NEAA) (Gibco, 11140050). The next day, each plate was triple transfected with 39.93 μg of total plasmid DNA encompassing pHelper, RepStop encoding the AAV2 Rep genes, and pUC19 at a ratio of 2:1:1, respectively, alongside 10 ng of assembled library DNA. The media was exchanged for fresh DMEM with 5% FBS and 1× NEAA at 20 h post-transfection. At 60 h post-transfection, the media and cell lysates were harvested and purified following the published protocol21.Individual recombinant AAVs were produced in suspension HEK293T cells, using F17 media (Thermo Fisher Scientific). Cell suspensions were incubated at 37 °C, 8% CO2, 125 RPM. At 24 h before transfection, cells were seeded into 200 mL of media at ~1 million cells/mL. The day after, the cells (~2 million cells/mL) were transfected with pHelper, RepCap, and transgene plasmids (2:1:1 ratio, 2 µg DNA per million cells) using Transport 5 transfection reagent (Polysciences) with a 2:1 PEI:DNA ratio. At three days post-transfection, cells were pelleted at 2000 RPM for 10 min into Nalgene conical bottles. The supernatant was discarded, and cell pellets were stored at −20 °C until purification. Each pellet, corresponding to 200 mL of cell culture, was resuspended in 7 mL of 500 mM NaCl, 40 mM Tris-base, and 10 mM MgCl2, with Salt Active Nuclease (ArcticZymes, #70920-202) at 100 U/mL. Afterward, the lysate was clarified at 2000 RCF for 10 min and loaded onto a density step gradient containing OptiPrep (Cosmo Bio, AXS-1114542) at 60%, 40%, 25%, and 15% at a volume of 5, 5, 6, and 6 mL respectively in OptiSeal tubes (Beckman, 361625). The step gradients were spun in a Beckman Type 70ti rotor (Beckman, 337922) in a Sorvall WX+ ultracentrifuge (Thermo Fisher Scientific, 75000090) at 340,252 × g for 1 h at 18 °C. Afterward, ~4.5 mL of the 40–60% interface was extracted using a 16-gauge needle, filtered through a 0.22 μm PES filter, buffer exchanged with 100K MWCO protein concentrators (Thermo Fisher Scientific, 88532) into PBS containing 0.001% Pluronic F-68, and concentrated down to a volume of 500 μL. The concentrated AAV was filtered through a 0.22 μm PES filter and stored at 4 °C or −80 °C.AAV titeringTo determine AAV titers, 5 μL of each purified AAV library were incubated with 100 μL of an endonuclease cocktail consisting of 11,000 U/mL Turbonuclease (Sigma T4330-50KU) with 1× DNase I reaction buffer (NEB B0303S) in UltraPure DNase/RNase-Free distilled water at 37 °C for 1 h. Next, the endonuclease solution was inactivated by adding 5 μL of 0.5 M EDTA, pH 8.0 (Thermo Fisher Scientific, 15575020) and incubated at room temperature for 5 min and then at 70 °C for 10 min. To release the encapsidated AAV genomes, 120 μL of a Proteinase K cocktail consisting of 1 M NaCl, 1% N-lauroylsarcosine, 100 μg/mL Proteinase K (Qiagen, 19131) in UltraPure DNase/RNase-Free distilled water was added to the mixture and incubated at 56 °C for 2–16 h. The Proteinase K-treated samples were then heat-inactivated at 95 °C for 10 min. The released AAV genomes were serial diluted between 460–460,000× in dilution buffer consisting of 1× PCR Buffer (Thermo Fisher Scientific, N8080129), 2 μg/mL sheared salmon sperm DNA (Thermo Fisher Scientific, AM9680), and 0.05% Pluronic F-68 (Thermo Fisher Scientific, 24040032) in UltraPure Water (Thermo Fisher Scientific). Then, 2 μL of the diluted samples were used as input in a ddPCR supermix (Bio-Rad, 1863023). Primers and probes, targeting the ITR and CAG promoter region, were used for titration, at a final concentration of 900 nM and 250 nM, respectively (ITR2_Forward: 5′-GGAACCCCTAGTGATGGAGTT-3′; ITR2_Reverse: 5′-CGGCCTCAGTGAGCGA-3′; ITR2_Probe: 5′-CACTCCCTCTCTGCGCGCTCG-3′ [FAM/Iowa Black FQ Zen]; CAG_Forward: 5′-TGTTCCCATAGTAACGCCAATAG-3′; CAG_Reverse: GTACTTGGCATATGATACACTTGATG-3′; CAG_Probe: 5′-TTACGGTAAACTGCCCACTTGGCA-3′ [FAM/Iowa Black FQ Zen]). Droplets were generated using a QX100 Droplet Generator following the manufacturer’s protocol. The droplets were transferred to a thermocycler and cycled according to the manufacturer’s protocol with an annealing/extension of 58 °C for 1 min. Finally, droplets were read on a QX100 Droplet Digital System to determine titers.Assessing production fitnessTo recover only encapsidated AAV genomes for downstream analysis, 1011 viral genomes were extracted using the endonuclease and Proteinase K steps outlined above (AAV Titering). After Proteinase K treatment, samples were column purified using a DNA Clean and Concentrator Kit (Zymo Research, D4033) and eluted in 25 μL elution buffer for NGS preparation.NGS sample preparationTo prepare AAV libraries for sequencing, qPCR was performed on the extracted AAV genomes or cDNA to determine the cycle thresholds for each sample type to prevent overamplification. PCR amplification using equal primer pairs (1–8) (described in ref. 22) was used to attach partial Illumina Read 1 and Read 2 sequences using Q5 Hot Start High-Fidelity 2× Master Mix with an annealing temperature of 65 °C for 20 s and an extension time of 60 s. Round one PCR products were purified using AMPure XP beads following the manufacturer’s protocol and eluted in 25 μL UltraPure Water (Thermo Fisher Scientific). Then, 2 μL was used as input in a second round of PCR to attach on Illumina adapters and dual index primers (NEB, E7600S) for five PCR cycles using Q5 Hot Start-High-Fidelity 2× Master Mix with an annealing temperature of 65 °C for 20 s and an extension time of 60 s. The round two PCR products were purified using AMPure XP beads following the manufacturer’s protocol and eluted in 25 μL UltraPure DNase/RNase-Free distilled water (Thermo Fisher Scientific).To quantify the amount of PCR products for NGS, an Agilent High Sensitivity DNA Kit (Agilent, 5067-4626) was used with an Agilent 2100 Bioanalyzer. PCR products were pooled and diluted to 2–4 nM in 10 mM Tris-HCl, pH 8.5 and sequenced on an Illumina NextSeq 550 following the manufacturer’s instructions using a NextSeq 500/550 Mid or High Output Kit (Illumina, 20024904 or 20024907), or on an Illumina NextSeq 1000 following the manufacturer’s instructions using NextSeq P2 v3 kits (Illumina, 20046812). Reads were allocated as follows: I1: 8, I2: 8, R1: 150, R2: 0.NGS data processingSequencing data was de-multiplexed with bcl2fastq (version v2.20.0.422) using the default parameters. The Read 1 sequence (excluding Illumina barcodes) was aligned to a short reference sequence of AAV9:CCAACGAAGAAGAAATTAAAACTACTAACCCGGTAGCAACGGAGTCCTATGGACAAGTGGCCACAAACCACCAGAGTGCCCAANNNNNNNNNNNNNNNNNNNNNGCACAGGCGCAGACCGGTTGGGTTCAAAACCAAGGAATACTTCCGAlignment was performed with bowtie2 (version 2.4.1)23 with the following parameters:–end-to-end –very-sensitive –np 0 –n-ceil L,21,0.5 –xeq -N 1 –reorder –score-min L,-0.6,-0.6 -5 8 -3 8Resulting sam files from bowtie2 were sorted by read and compressed to bam files with samtools (version 1.11-2-g26d7c73, htslib version 1.11-9-g2264113)24,25.Python (version 3.8.3) scripts and pysam (version 0.15.4) were used to extract the 21-nucleotide insertion from each amplicon read. Each read was assigned to one of the following bins: Failed, Invalid, or Valid. Failed reads were defined as reads that did not align to the reference sequence, or that had an indel in the insertion region (i.e., 20 bases instead of 21 bases). Invalid reads were defined as reads whose 21 bases were successfully extracted but matched any of the following conditions: (1) Any one base of the 21 bases had a quality score (AKA Phred score, QScore) below 20, i.e., error probability >1/100, (2) any one base was undetermined, i.e., “N”, (3) the 21 base sequence was not from the synthetic library (this condition does not apply to NNK library variants). Valid reads were defined as reads that did not fit into either the Failed or Invalid bins. The Failed and Invalid reads were collected and analyzed for quality control purposes, and all subsequent analyses were performed on the Valid reads.Count data for valid reads was aggregated per sequence, per sample, and was stored in a pivot table format, with nucleotide sequences on the rows, and samples (Illumina barcodes) on the columns. Sequences not detected in samples were assigned a count of 0.Data normalizationCount data was reads-per-million (RPM) normalized to the sequencing depth of each sample (Illumina barcode) with:$${r}_{i,j}=\frac{{k}_{i,j}}{{\sum }_{l=1}^{n}{k}_{l,j}}\times 1,000,000$$
(1)
Where r is the RPM-normalized count, k is the raw count, i = 1 … n sequences, and j = 1 … m samples.As each biological sample was run in triplicate, we aggregated data for each sample by taking the mean of the RPMs across p replicates of sample s:$${\mu }_{i,s}=\frac{{\sum }_{l=1}^{p}{r}_{i,l}}{p}$$
(2)
We estimated normalized variance across replicates by taking the coefficient of variation (CV):$${{CV}}_{i,s}=\frac{{\mu }_{i,s}}{{\sigma }_{i,s}}$$
(3)
where \({\sigma }_{i,s}\) is the standard deviation for variant i in sample s over p replicates.Log2 enrichment for each sequence was defined as:$${e}_{i,s}={\log }_{2} \left(\frac{{\mu }_{i,s}}{{\mu }_{{corrected},i,t}}\right)$$
(4)
where e is the \({\log }_{2}\) enrichment, \(\mu\) is the mean of the replicate RPMs, and t is the normalization sample. For production fitness, the sample s is the variant abundance after AAV production, and the normalization sample t is the variant abundance in the plasmid pool. For functional screens, the sample s is the variant abundance of the screen, and the normalization factor t is the variant abundance after AAV production.To avoid dividing by 0 in e (for NNK library processing), \({\mu }_{{corrected}}\) is defined as:$${\mu }_{{corrected},i,t}=\left\{\begin{array}{cc}{\mu }_{i,t},& {\mu }_{i,t} \, > \, 0 \\ \frac{1}{{\sum }_{l=1}^{n}{k}_{l,t}},& {\mu }_{i,t} \,=\, 0\end{array}\right.$$
(5)
Counts of 0 across all three replicates for the normalization sample were adjusted to a pseudocount of 1 across all three replicates.Production fitness training and assessmentWe designed a robust ML framework for the production fitness and Fit4Function functional mappings. A long short-term memory (LSTM) regression model with two hidden layers of 140 and 20 nodes was implemented in Keras26. RNNs, and LSTMs in particular, have been successfully applied for learning functions from biological sequence data as they are designed to capture local and distant relationships across different parts of the input sequences13,27. Model parameters and hyperparameters were subject to fine-tuning processes but no significant performance was gained across all different functional fitness models implemented in this study. Therefore, we kept a simple model architecture across all modeling throughout this study. The input layer was 7-mer amino acid sequences one-hot encoded into a 20 × 7 matrix. The target/output is the relative production or functional fitness score. The loss was optimized by mean-squared error with the Adam optimizer running on a learning rate of 0.00128. The batch size was set to 500 observations. To avoid overfitting, model training was controlled by a custom early stopping procedure where the training process was terminated if the ratio of training error to validation error dropped below 0.90.For production fitness learning, the training size was optimized by training the framework on increments of 1K variants. Variants that were not detected (n = 5693) after AAV production were filtered out from training. Model validation performance was reported at each training size, and a size of 24K variants was arbitrarily selected for final model training given that the model performance reached a plateau after a training size of ~5K. The modeling library core variants (N ~ 60K after removing the non-detected sequences) were then randomly divided into training (24K), validation (5K), and testing subsets (30.6K), all from the modeling library. The model was trained on the training set, validated during the training process on the validation set, and tested on the testing set. The model was further tested on the unique variants from the assessment library to assess its generalization across libraries.Fit4Function library samplingThe Fit4Function libraries are intended to be sampled from the production-fit space. For the Fit4Function library utilized in this study, we first uniformly sampled a set of 7-mer amino acid sequences 100× the required library size (240K Fit4Function variants * 100 = 24M variants), by equally sampling each amino acid at each of the seven positions. Duplicates were removed and the remaining sequences were scored using the production fitness model. Then, the 240K Fit4Function library variants were probabilistically sampled from the parametrized production-fit distribution. In addition to the 240K production-fit variants, we added 1K stop-codon-containing variants and 3K variants from the 10K variants shared between the modeling and assessment libraries as a control set.Fit4Function library validationFitness enrichment scores are relative across library variants due to normalization calculations; calibration is needed to make the fitness scores of libraries of different compositions comparable for assessment or integration purposes. Therefore, we used the 3K control set in the Fit4Function library to fit an ordinary linear regression model of the measured production fitness scores between the Fit4Function library and the modeling library. These regression parameters were applied to the production fitness measured scores of the 240K Fit4Function variants to obtain calibrated production fitness scores. After synthesizing the Fit4Function library, we compared, by means of correlation, the predicted fitness scores to the calibrated measured fitness.AAV mouse in vivo biodistribution assaysPurified AAV libraries were injected intravenously (retro-orbital sinus) into four C57BL/6J adult female mice at a dose of 1 × 1012 vg/mouse (6–8 weeks of age). Two hours post-injection, serum was collected and tissues were harvested using disposable 3 mm biopsy punches (Integra, 33-32-P/25) with a new biopsy punch used per organ per replicate. Harvested tissues were immediately frozen on dry ice. AAV genomes were recovered using a DNeasy kit (Qiagen, 69504) following the manufacturer’s protocol, and samples were eluted in 200 μL elution buffer for NGS preparation.AAV cynomolgus macaque in vivo biodistribution assaysThe library administered had 100K unique amino acid variants following the Fit4Function criteria (uniformly sampled from the production-fit sequence space) in addition to a calibration set (3K), control variants, and AAV9. Each variant in the Fit4Function distribution was represented by either two or six codon replicates; AAV9 was represented by two codon replicates. The purified AAV library was injected intravenously at a dose of 4.6 × 1012 vg/kg into a 3.4 kg adult female cynomolgus macaque that was pre-screened for neutralizing antibodies (NAb) against AAV9 (CRL). Four hours after systemic delivery, the animal was sedated with an intramuscular injection of a combination of ketamine hydrochloride and acepromazine, and anesthetized by a mixture of isoflurane and oxygen. Next, the animal was perfused with 0.9% sodium chloride, and tissues were harvested and snap-frozen on dry ice. DNA was extracted using a DNeasy kit in a Qiagen QIAcube Connect. Samples were then processed as detailed in the “NGS sample preparation” section.AAV rhesus macaque in vivo transduction assaysApproximately 3-month-old rhesus macaques (~1 kg; one male, one female) were screened and assigned to the project after confirming seronegative status for AAV9 antibodies using standardized Testing Center assays. Sedation with Telazol (IM) was performed prior to intravenous administration of a purified AAV library (1 × 1013 vg/kg) with blood samples collected (~4 mL; hematology, clinical chemistry, serum, plasma; pre-administration and weekly post-administration). Animals were monitored closely during the study period and until the endpoint (four weeks post-administration) and remained robust and healthy with no evidence of adverse findings. Four weeks after systemic delivery, tissues were collected and snap-frozen over liquid nitrogen then placed on dry ice immediately prior to storage at ≤ − 80 °C. RNA and DNA were extracted using TRIzol (Invitrogen, 15596026) following the manufacturer’s instructions. Total RNA was then processed through an RNeasy kit (Qiagen, 74106) followed by on-column DNA digestion. RNA was converted to cDNA using Maxima H Minus Reverse Transcriptase (Thermo Fisher Scientific, EP0751) according to the manufacturer’s instructions. Samples were then processed as detailed in the “NGS sample preparation” section.Rhesus macaque serum screening for anti-AAV9 neutralizing antibodiesNeutralization assays were performed at 500 or 1000 vector genomes (vg)/cell in Perkin–Elmer white 96-well plates. Four-fold serial dilutions (1:4 to 1:16,384) of macaque serum samples were prepared in 96-well plates using DMEM supplemented with 5% fetal bovine serum (FBS). Then, 40 μL of each dilution was transferred to a separate 96-well plate, mixed with an equal volume of AAV9:CAG-GFP-P2A-Luciferase-WPRE-SV40 vector (4–8 × 107 vg per 40 μL, diluted in DMEM-5% FBS), and incubated for 1 h at 37 °C. Following the incubation, AAV-serum samples were transferred into a new 96-well plate (20 μL triplicates) and a total of 80 μL of DMEM-5% fetal bovine serum, containing 20,000 HEK293T cells, was added to each well (final volume of 100 μL). 96-well plates were incubated for 48 h at 37 °C, 5% CO2. Luminescence levels were read using a Perkin–Elmer Victor Luminescence Plate Reader using the britelite plus Reporter Gene Assay System (Perkin–Elmer, #6066761). Data was analyzed using the neutcurve Python package developed by the Bloom laboratory. The NAb titer was measured as the concentration that resulted in a 50% reduction in luciferase activity relative to the no-serum control. Animals used in the transduction study had NAb titers of <1:12 in this set of antibody screens.In vitro binding and transductionHEK293T/17 (ATCC® CRL-11268™), HEPG2 (ATCC® HB-8065™), THLE-2 (ATCC® CRL-2706™), hCMEC/D3 (Millipore, SCC066), and human and mouse BMVECs (Cell Biologics, H-6023 and C57-H6023) were used in this study. Among these, HEK293T/17 is the only cell line known to be potentially cross-contaminated with HeLa cells29. We used morphology checks under light microscopy with a Leica DM IL LED Inverted Laboratory Microscope to rule out cross-contamination with HeLa cells. Cells were grown in 100 mm dishes and exposed to the Fit4Function or NNK 7-mer library (1 × 104 vg/cell for HEK293T/17, 3 × 104 vg/cell for hCMEC/D3, 6 × 104 vg/cell for primary human and mouse BMVECs, and 5 × 103 vg/cell for HEPG2 and THLE-2) diluted in 10 mL of growth media at 4 °C with gentle rocking for 2 h. Cells were then washed 3× with DPBS, and total DNA was extracted with the DNeasy kit (Qiagen) according to the manufacturer’s instructions. Half of the recovered DNA was used in PCR amplification for viral genome sequence recovery.Transduction assays were performed as described above with the following exceptions: The cells were cultured in growth media containing AAV for 60 h and total RNA was then extracted with the RNeasy kit (Qiagen). From the total RNA, 5 μg was converted to cDNA using the Maxima H Minus Reverse Transcriptase according to the manufacturer’s instructions.Sequence-to-function mappingFunctional scores were quantified as the log2 of the fold-change enrichment of the variant reads-per-million (RPM) after the screen relative to its RPM in the AAV library, i.e., log2 (Assay RPM/Virus RPM). Fit4Function models utilized the same design of the ML framework utilized for production fitness mapping (two-layer LSTM, custom early stopping, batch size of 500 variants, MSE error, and Adam optimizer). Out of the 240K variants in the Fit4Function library, 90K were allocated for training and testing the ML function models (model construction), and 150K variants were held out for validation of the MultiFunction approach. The training size for each model was optimized independently. As with the production fitness model, the liver-related functional fitness models were assessed by correlation between the predicted and measured functional scores.MultiFunction library designUsing the previously generated fitness models of the production fitness and the five liver-related functional fitness models described in the main text, we conducted an in silico screen of 10M randomly sampled 7-mer sequences to identify variants that are highly fit for all six traits. The threshold of high fitness for each function was arbitrarily set to the 50th percentile of each functional fitness distribution from the Fit4Function screening data. The percentiles were calculated on the detected variants of each functional assay from the 90K model construction data set. To reduce false positive predictions, we arbitrarily increased the threshold for predicted hits by 5% of the fitness dynamic range for each function. These thresholds were used to filter the 10M variants that were run through the six functional prediction models for variants predicted to pass the six modified thresholds. Of these, we sampled 30K variants to be included in the MultiFunction library. The 30K variants were each represented by two codon replicates. The MultiFunction library also included (1) a positive control set (3K) that was drawn from the subset of the 150K Fit4Function validation set that met the six conditions on the actual measurements, (2) a set of 10K variants randomly sampled from the Fit4Function 240K core variants as background controls representing the production-fit space, (3) a set of 3K calibration variants present in both the Fit4Function library and the modeling library to be used as background controls representing the entire (unbiased) sequence space, and (4) 1K stop-codon-containing sequences.MultiFunction library validationThe MultiFunction library was synthesized, AAV was produced, and the five liver-related functions were screened for as was done with the Fit4Function library. We quantified the success rate of the MultiFunction library in terms of hit rate, i.e., out of the 30K variants predicted to meet the six criteria, what percentage satisfied the six criteria when the MultiFunction library was screened on those functions (predicted positive versus measured positive). To determine whether a variant meets specific functional criteria, we compared the distribution of that function for the MultiFunction variants against the positive control set. For a variant to be considered a hit for a given function, its measured enrichment score had to be greater than the mean − 2 SD (standard deviations) of the enrichment scores of the positive control set measured in the same experiment. This thresholding was applied to avoid the overestimation of the hit rate due to outliers in the positive control set. A variant is considered a hit in calculating the MultiFunction hit rate only if it is a hit for all six functions; a variant that meets only five or fewer conditions is not considered a hit.As the sample sizes for the Fit4Function and the Uniform control sets in the MultiFunction library are small (3K and 10K, respectively), we estimated the multi-trait hit rate from the 240K Fit4Function library to provide a more confident assessment of the percentage of variants in the production-fit space and the uniform space that meet all six criteria. We calculated the hit rate of the Fit4Function space as the percentage of the production-fit variants that are measured to pass the thresholds for all six functions simultaneously (17.14K out of the 240K Fit4Function variants). The thresholds used are the same as those used to define the multi-trait variant set from which the positive control set was drawn. As the uniform libraries in the study (the Modeling and Assessment libraries) were not used in screens for the five liver-related functions, we could not use them to assess the Uniform hit rate. Instead, we estimated the hit rate for the uniform space as the hit rate in the production-fit space scaled by the percentage of the uniform space occupied by the production-fit space. In other words, Uniform hit rate = Fit4Function hit rate × production-fit ratio = 7.0% × 37.3% = 2.6%, given that no non-fit variants meet the criteria for production fitness. For comparison, when the Uniform multi-trait hit rate was calculated from the small Uniform control sets in the Fit4Function and the MultiFunction libraries, we obtained hit rates of 0.0% and 4.6%, respectively. When the Fit4Function hit rate was calculated from the Fit4Function control set in the MultiFunction library, we obtained a hit rate of 15.4%. We believe that the hit rates calculated with the larger Fit4Function library are more accurate.Individual capsid characterizationIndividual capsids were cloned into iCAP-AAV9 (K449R) backbone (GenScript), produced as described above, and administered to C57BL/6J (The Jackson Laboratory, 000664) mice at a dose of 1 × 1010 vg/mouse (n = 5/group). Three weeks later, three separate lobes of the liver were collected for RNA extraction and a single lobe per mouse was fixed in 4% PFA.For microscopy, fixed liver tissues were sectioned at 100 µm using a Leica VT1200 vibratome. Sections were mounted with ProLong™ Gold Antifade Mountant with DAPI (Thermo Fisher Scientific, P36931). Liver images were collected using the optical sectioning module on a Keyence BZ-X800 with a Plan Apochromat 20× objective (Keyence, BZ-PA20). Three images were taken for each animal (n = 5/group) and compared to a no-injection control (n = 3 animals). In CellProfiler, nuclei were segmented and DAPI+ nuclei were identified using a threshold on DAPI intensity determined from the no-injection control. Each DAPI+ nuclei was then quantified with the median pixel intensity in the GFP channel.For assessment of liver transduction by quantitative RT-PCR, total RNA was recovered using TRIzol (Invitrogen, 15596026) following the manufacturer’s instructions. Total RNA was then processed (RNeasy kit, Qiagen, 74106) followed by on-column DNA digestion. RNA was converted to cDNA using Maxima H Minus Reverse Transcriptase (Thermo Fisher Scientific, EP0751) according to manufacturer instructions. Afterward, qPCR was used to detect AAV-encoded RNA transcripts with the following primer pair (5′-GCACAAGCTGGAGTACAACTA-3′) and (5′-TGTTGTGGCGGATCTTGAA-3′) and the following primer pair for GAPDH (5′-ACCACAGTCCATGCCATCAC-3′) and (5′-TCCACCACCCTGTTGCTGTA-3′).THLE-2 and HEPG2 cells were seeded in a 96-well plate the day before adding the AAVs at 5000 vg/cell. For binding assays, AAVs were diluted in media and incubated with cells at 4 °C with gentle shaking for 1 h. After incubation, cells were washed 3× with PBS to remove unbound AAV particles and treated with Proteinase K to release AAV genomes for qPCR quantification. For transduction assays, cells were incubated with the AAVs for 24 h at 37 °C and assayed with Britelite plus (Perkin–Elmer, 6066766) following the manufacturer’s protocol.Interpreting the production fitness and liver-related functional fitness modelsWe first assessed the effect of each amino acid residue on each predicted function. We fixed each amino acid at each position in the 7-mer one at a time while randomizing the residues at the rest of the positions as in ref. 14. However, we used regression models whereas the referenced study used classification models. Therefore, we used a different importance summarization metric. We used the production fitness and liver-related functional fitness models to score the functional fitness of the randomized variants (n = 100K) to create a distribution in each scenario. The randomized variants were filtered for only production-fit variants using our production fitness model; this was done to eliminate non-production-fit variants that would impact the accuracy of the models trained to only make predictions in the context of production-fit variants (this is a key feature of the Fit4Function approach). We used the mean of the distributions in each case to indicate the overall effect of each fixed residue on each function. The mean predictions for each distribution were z-score normalized to predictions made for randomized 7-mer sequences with no fixed residues at any positions.Next, we assessed the effect of different physicochemical properties on the predicted functions. For this, we built six surrogate concept models30 to interpret each of the six functional models. Unlike the original models that are trained on the amino acid sequence of each variant, each concept surrogate model is trained on features representing the physicochemical properties (concepts) of the training variants. To learn how much of those concepts are captured by the original models, the surrogate models are trained to predict the predictions of the original models. In other words, while the original models are trained to learn to predict functional enrichment (fitness) from the sequences of variants, the surrogate models are trained to predict what the original models would predict from the concept features of variants. If the surrogate model can closely learn (replicate) the predictions of the original model by training on the concept features (not the sequences), this means that the original model made predictions largely from those concepts. We built the surrogate concept models as linear regression lasso models, with 10-fold cross-validation for a training set of an arbitrary size of 120K variants, trained on each of twelve physicochemical properties of the 7-mer sequences (i.e., a total of 84 features): the volume of the side chain (normalized), hydropathy (normalized), polarity, hydrogen donor, hydrogen acceptor, positive charge, negative charge, aliphatic, aromatic, sulfur, hydroxyl, and amide. We used the same properties considered in a prior interpretability study14 but with the assumption that histidine is not positively charged in the pH of the environments relevant to our functions. We used R2 to assess the extent to which the variation of the original model’s predictions was explained by those of the surrogate models.Finally, we assessed the effect of first- and second-degree interactions between residues in the 7-mer on the predicted functions. We built a set of surrogate concept linear regression lasso models with the hot-encoding of first-degree interactions (20 amino acids × 7 positions = 140 features) and another set of models with the hot-encoding of both first- and second-degree interactions (400 combinations of amino acids at two positions × 21 possible dual-position combinations = 8400 features). In each case, data from the 150K Fit4Function variants in the screens of the six functions, excluding non-detected variants, were used for training. R2 was used to assess the extent to which the variation of the original model’s predictions was explained by those of the surrogate models.Statistics and reproducibilityMicrosoft Excel and GraphPad Prism 9 were used for the analysis of experimental data. Python was used for large computational analyses and modeling. The reproducibility of the production fitness scores was tested via 10K variants common to the assessment library and modeling library as shown in Fig. 2d. Biological triplicates were performed for the binding or transduction of different cell lines using the Fit4Function library versus an NNK library and their reproducibility was quantified via pairwise Pearson correlation as shown in Fig. 3c and Supplementary Fig. 5. The Fit4Function library was screened in four mice, and the reproducibility of the biodistribution in eight tissues was assessed in the form of replication quality between pairs of animals as shown in Fig. 3e and Supplementary Fig. 4. The MultiFunction library was screened across in vitro (three replicates for production fitness, four replicates for binding and transduction) and in vivo (n = 3 mice) assays, and the reproducibility within each assay was assessed in the form of replication quality between pairs of replicates as shown in Supplementary Fig. 6. For the comparison between AAV9 and the seven MultiFunction variants, in vitro cell transduction was assessed with four replicates per group (Fig. 4e and Supplementary Fig. 7d; error bars show the standard deviation from the mean) and unpaired, one-sided t-tests were conducted on log-transformed values, with Bonferroni correction for multiple hypotheses. The in vivo characterization was performed using five female mice per group except for the no-injection control group (Supplementary Fig. 7b and c, n = 3 mice; error bars in Supplementary Fig. 7c show the standard deviation from the mean), experimenters were not blinded to the sample groups, and unpaired, one-sided t-tests were conducted on log-transformed values, with Bonferroni correction for multiple hypotheses. Other statistical tests are described in the text. No animals or samples were excluded from the analysis. The variant BI152 was unintentionally excluded from the rhesus macaque study due to a library assembly error. The number of replicates was chosen based on prior data that indicated large effect sizes. ML models were trained in repetition with randomized subsampling and parameterization to ensure reproducibility. Training and testing sample sizes are described in the “Methods” section for each model. Models were tested using independent (blind) datasets where possible, i.e., testing was conducted using an independent assessment library for production fitness (Fig. 2f) and independent measurements in a separate animal (Fig. 3f).Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.