extfieldsofstudy: Biology Medicine

text:
Multifunctional CRISPR-Cas9 with engineered immunosilenced human T cell epitopes

The CRISPR-Cas9 system has raised hopes for developing personalized gene therapies for complex diseases. Its application for genetic and epigenetic therapies in humans raises concerns over immunogenicity of the bacterially derived Cas9 protein. Here we detect antibodies to Streptococcus pyogenes Cas9 (SpCas9) in at least 5% of 143 healthy individuals. We also report pre-existing human CD8+T cell immunity in the majority of healthy individuals screened. We identify two immunodominant SpCas9 T cell epitopes for HLA-A*02:01 using an enhanced prediction algorithm that incorporates T cell receptor contact residue hydrophobicity and HLA binding and evaluated them by T cell assays using healthy donor PBMCs. In a proof-of-principle study, we demonstrate that Cas9 protein can be modified to eliminate immunodominant epitopes through targeted mutation while preserving its function and specificity. Our study highlights the problem of pre-existing immunity against CRISPR-associated nucleases and offers a potential solution to mitigate the T cell immune response.

T he Clustered Regularly Interspaced Short Palindromic Repeat (CRISPR)-Cas9 technology has raised hopes for developing personalized gene therapies for complex diseases such as cancer as well as genetic disorders, and is currently entering clinical trials 1,2 . The history of gene therapy has included both impressive success stories and serious immunologic adverse events [3][4][5][6][7][8] . The expression of Streptococcus pyogenes Cas9 protein (SpCas9) in mice has evoked both cellular and humoral immune responses 9,10 , which raises concerns regarding its safety and efficacy as a gene or epi-gene therapy in humans. These preclinical models and host immune reactions to other exogenous gene delivery systems [11][12][13] suggest that the pathogenic, non-self origin of Cas9 may be immunogenic in humans.
Both B cell and T cell host responses specific to either the transgene or the viral components of adenoviral 14,15 and adenoassociated viral (AAV) 11,12 vectors have been detected, despite relatively low immunogenicity of AAV vectors. In the case of AAV, specific neutralizing antibodies (Abs) and T cells are frequently detected in healthy donors [16][17][18][19] and specific CD8+T cells have been shown to expand following gene delivery 18 . There has been recent progress in developing strategies to overcome this problem, such as capsid engineering and transient immunosuppression [20][21][22] . The potential consequences of immune responses to expressed proteins from viral vectors or transgenes include neutralization of the gene product; destruction of the cells expressing it, leading to loss of therapeutic activity or tissue destruction; induction of immune memory that prevents readministration; and fulminant innate inflammatory responses 23,24 . More potent immune responses to gene therapies have been observed in humans and non-human primate models compared to mice 8,25 .
Of the Cas9 orthologs derived from bacterial species, the SpCas9 is the best characterized. S. pyogenes is a ubiquitous pathogen, with an annual incidence of 700 million worldwide 26 , but the field is only now beginning to explore potential immunity to SpCas9 in humans 27,28 . CRISPR application for human therapies will span its use both for gene editing (through DNA double-strand breaks) or epigenetic therapies (without DNA double-strand breaks). In fact, recent reports shed light on CRISPR's ability to activate or repress gene expression in mice [29][30][31] , which opens the door to a variety of new therapeutic applications such as activating silent genes, compensating for disrupted genes, cell fate reprogramming, or silencing disrupted genes, without the concern over permanent change in DNA sequence. However, unlike the use of Cas9 for gene editing, which may only require Cas9 presence in cells for a few hours, current techniques for CRISPR-based epigenetic therapies require longer term expression of Cas9 in vivo, possibly for weeks and months 30,31 , which poses the challenge of combating pre-existing immune response towards Cas9. This challenge will need to be addressed before CRISPR application for human therapies, especially for epigenetic therapies, can be fully implemented. Delivery of CRISPR in vivo by incorporating its expression cassette in adeno-associated virus (AAV), will most likely shape many of the initial clinical trials as AAV-based gene delivery is one of the safest and most prevalent forms of gene therapies in human. AAV will enable longer term expression of Cas9, desirable for epigenetic therapies. Therefore, it is highly likely that CRISPR delivery through AAV and its expression within target cells will engage CD8+T cell immunity.
Here, we seek to characterize the pre-existing immune response to SpCas9 in healthy individuals and to identify the immunodominant T cell epitopes with the aim of developing SpCas9 proteins that have diminished capacity to invoke human adaptive response. We identify two immunodominant SpCas9 T cell epitopes for HLA-A*02:01 by an improved prediction algorithm and T cell assays using healthy donor PBMCs. We demonstrate that Cas9 protein can be modified to eliminate immunodominant epitopes through targeted mutation while preserving its function and specificity.

Results
Cas9-specific serum antibodies in healthy controls. We first determined whether healthy individuals have detectable IgG Abs to SpCas9. Of 143 healthy control sera screened, 82 (57.3%) had detectable Abs against S. pyogenes lysate using ELISA (Fig. 1a). Sera with the highest reactivity to S. pyogenes lysate (n = 80) were screened for Abs against recombinant SpCas9, of which 7 (8.8%) were positive (p < 0.0001, two-tailed t-test). At least 5.0% of healthy individuals screened in this study had Cas9-specific Abs (Fig. 1a).
Cas9 candidate T cell epitope prediction. Whether Cas9-specific antibodies impact the efficacy or safety of CRISPR application in human remains to be seen. However, cellular immunity is expected to have a more significant impact in case of CRISPR delivery through viruses. The Cas9 expression cassette, delivered by a viral vector, leads to intracellular expression of this protein in target cells, which could evoke a cellular immune response. We thus focused on investigating the T cell immune response against SpCas9. We predicted HLA-A*02:01-restricted T cell epitopes derived from SpCas9 using a model that uses both MHC binding affinity and biochemical properties of immunogenicity 32 (Supplementary Table 1; the top 5 are shown in Fig. 1b). This model incorporates T cell receptor contact residue hydrophobicity and HLA binding prediction, which enhances the efficiency of epitope identification, as we previously reported 32 . We plotted the calculated normalized binding (S b ) and immunogenicity (S i ) scores for each peptide (Fig. 1c) to predict the more immunogenic epitopes, which are expected to have both high HLA binding (low S b ) and more hydrophobicity (high S i ). We chose HLA-A*02:01 because it is the most common HLA type in European/North American Caucasians.
T cell epitope mapping of Cas9. We then investigated whether peripheral blood mononuclear cells (PBMCs) derived from healthy individuals had measurable T cell reactivity against the predicted SpCas9 MHC class I epitopes. We synthesized 38 peptides (Supplementary Table 1) and grouped them into 10 pools of 3-4 peptides each. We measured peptide-specific T cell immunity using IFN-γ secretion ELISpot assays with PBMCs derived from 12 healthy individuals (HLA-A*02:01, n = 10; non-HLA-A*02:01, n = 2) and identified immunoreactive epitopes within pools 3 or 5 in 83.0% of the donors tested (90% of the HLA-A*02:01 donors; Fig. 1d). The seven individual peptides from pools 3 and 5 were evaluated by IFN-γ ELISpot and the dominant immunogenic epitopes were SpCas9_240-248 and SpCas9_615-623, designated peptides α and β, from pools 5 and 3, respectively. The subdominant epitopes were found to be γ and δ from pools 3 and 5, respectively. Both peptides α and β are located in the REC lobe of the Cas9 protein (Fig. 1e) that binds the sgRNA and the target DNA heteroduplex 33 . The individual peptides within pools that were positive for any donor were evaluated for this donor by IFN-γ ELISpot. The immunoreactivity and position of the 38 predicted peptides (a few of which are overlapping) within the Cas9 protein are shown in Fig. 1e.
Peptides α and β are shown as red dots on the epitope prediction plot (Fig. 1c) and their sequences and predicted ranking are shown in Fig. 1b and Supplementary Table 1. As predicted, these peptides had low S b and high S i values. Both the immunodominant (α and β) and subdominant (γ and δ) T cell epitopes identified by IFN-γ ELISpot were within the top 5 most immunogenic epitopes predicted by our immunogenicity model 32 . Their ranking as predicted by the consensus method hosted on the IEDB server using default settings was 14, 5, 18, and 4, respectively. Sequence similarity of peptides α and β to amino acid sequences in known proteins was investigated using Protein BLAST and the IEDB epitope database 34 . This was done to investigate whether there is any chance that the T cell immune response that we are detecting in healthy individuals could be due to previous exposure to another protein of similar sequence. A peptide was considered 'similar' to α or β if no more than 2 of 9 amino acid residues (that are not the second or ninth) were not matching (78% similarity). None of these two peptides resembled known epitopes in the IEDB database, but similarity to other Cas9 orthologs and other bacterial proteins was detected (Supplementary Tables 2 and 3). Epitope β has sequence   Mutated Cas9 proteins have lower immune recognition. We next hypothesized that mutation of the MHC-binding anchor residues of the identified immunogenic epitopes would abolish specific T cell recognition (Fig. 2a). The epitope anchor residues (2nd and 9th) are not only necessary for peptide binding to the MHC groove, but are also crucial for recognition by the T cell receptor 32 . The percentage of CD8+β pentamer+T cells decreased to 0.3% when APCs were pulsed with the mutated peptide (β2; Fig. 2b) compared with 3.09% with the wild type peptide (β; Fig. 2a). We then examined the reactivity of healthy donor T cells to modified peptides α or β with mutations in residues 2, 9, or both (sequences are shown in Fig. 2c) using IFN-γ ELISpot assay. The epitope-specific T cell reactivity was W T -C a s 9 C a s 9 -α 2 C a s 9 -β 2 C a s 9 -α 2 β 2 DN DN Fig. 2 SpCas9 immunodominant epitope-specific CD8+T cell recognition is abolished after anchor residue mutation. a Epitope β-specific CD8+T cell response detected using β-specific pentamer in PBMCs stimulated with peptide β-pulsed antigen presenting cells. b The percentage of CD8+ pentamerβ+T cells was reduced to 0.3% when healthy donor B cell APCs were pulsed with the mutated peptide β2. c Positions, sequences, and IEDB HLA binding percentile rank of epitopes α and β before and after mutation of the anchor (2nd and/or 9th) residues. S b , normalized binding score; S i , normalized immunogenicity score. d Representative IFN-γ ELISpot assay in triplicate wells comparing T cell reactivity to wild type or mutated epitopes α and β. These results are representative of 12 donors and two independent replicates (data from all 12 donors are shown in Supplementary Fig. 1). e, f IFN-γ ELISpot comparing T cell reactivity to APCs expressing WT or modified Cas9 proteins. APCs expressing FluM1 were used as a positive control. APCs expressing GAPDH or spiked with peptide α2 were used as negative controls. Data represent mean ± SEM of 5 replicates (right). Statistical analysis was performed post hoc and results are exploratory. Source data are available in the Source Data file markedly reduced with the mutant peptides (Supplementary Fig. 1; representative ELISPOT well images are shown in Fig. 2d). The average reduction for the responsive HLA-A*02:01 donors was 25-fold from α to α29 (n = 7; p < 0.03) and 30-fold from β to β29 (n = 8; p < 0.03; Benjamini-Hochberg; Supplementary Fig. 1; Supplementary Table 4). The predicted binding affinity to MHC class II was also decreased for α2 and β2 epitopes, although the experimental significance of this alteration is unknown. We then generated modified Cas9 constructs by mutating the second residue of peptide α (L241G; Cas9-α2), peptide β (L616G; Cas9-β2), or both (Cas9-α2β2). To measure the effect of mutating the anchor residue of the immunogenic epitopes on T cell recognition of the Cas9 protein, we transiently transfected healthy donor B cell APCs with mRNA encoding wild type Cas9 (WT-Cas9), Cas9-α2, Cas9-β2, or Cas9-α2β2. The T cell response measured by IFN-γ ELISpot after coculturing of Cas9 transfected APCs with autologous PBMCs was significantly decreased for the modified Cas9 proteins (Fig. 2e, f). Introduction of the β2 mutation was the most effective in reducing T cell immunogenicity (5.5-fold, p < 0.0001, two-tailed t-test). This mutation in the REC1 domain (Figs. 1e and 3a) is not located in any of the two regions that are absolutely essential for DNA cleavage, the repeatinteracting (97-150) and the anti-repeat-interacting (312-409) regions 33 . These results demonstrate that mutating the anchor amino acid of a highly immunogenic epitope can influence the overall immunogenicity of Cas9. Thus, engineering Cas9 variants with reduced immunogenicity potential can be used in conjunction with other strategies for safer CRISPR therapies and even possibly reduce the dosage of systemic immunosuppression needed for patients. We therefore tested the function of Cas9-β2 in comparison with WT-Cas9 in the context of DNA cleavage and transcriptional modulation. To examine the nuclease activity of Cas9-β2 and compare with WT-Cas9, we targeted Cas9-β2 or WT-Cas9 to an endogenous locus (EMX-1) and measured percent indel formation (Fig. 3b, c). Our data demonstrate that Cas9-β2 retains nuclease capacity in the locus we studied as well as on a synthetic promoter (Fig. 3c, Supplementary Fig. 2A, B).
Next, we determined whether Cas9-β2 can successfully recognize and bind its target DNA leading to transcriptional modulation. We first tested this in the context of enhanced transgene expression from a synthetic CRISPR responsive promoter in HEK293 cells using 14nt gRNAs and aptamermediated recruitment of transcriptional modulators similar to what we had shown before ( Supplementary Fig. 2C, D). Having shown successful transgene activation, we then investigated whether this variant retains such capacity within the chromosomal contexts of endogenous genes. We transfected the cells with plasmids encoding Cas9-β2 or WT-Cas9 and 14nt gRNAs against two different endogenous genes (TTN and MIAT). qRT-PCR analysis showed that this variant successfully led to target gene expression (Fig. 3d-f). To further characterize Cas9-β2 specificity, we performed genome-wide RNA sequencing after targeting Cas9-β2 or WT-Cas9 to the MIAT locus for transcriptional activation. The results demonstrated no significant increase in undesired off-target activity by Cas9-β2 as compared to WT-Cas9 (Fig. 3g).
To show the extensibility of our approach, we tested the function of Cas9-α2, that has a mutation located in the REC2 domain (Figs. 1e and 3a). Cas9-α2 also demonstrated DNA cleavage and transcriptional modulation functionality comparable with WT-Cas9 ( Supplementary Fig. 3A-E). This is consistent with a previous study which showed that Cas9 with a deleted REC2 domain retains its nuclease activity 33 . When T cells were stimulated with APCs spiked with peptide α2, the percentage of CD8+CD137+T cells (a marker of T cell activation 35 ) was decreased by 2.3-fold as compared to WT peptide α stimulation ( Supplementary Fig. 3F).
Immune responses to MHC Class II Cas9 epitopes. We also predicted MHC class II binding epitopes for the SpCas9 protein to HLA-DRB1 (10 alleles), HLA-DQ (5 alleles), and HLA-DP (8 alleles) using the IEDB analysis tool (Supplementary Data 2). For MHC class II, epitope α is predicted to be a top binder to HLA-DRB1*01:02 and epitope β a top binder to HLA-DPA1*01:03 and DPB1*02:01 (Supplementary Data 2). We selected and synthesized SpCas9 long peptides that include epitopes in the top 2% of predicted MHC class II binders (sequences and alleles shown in Supplementary Table 6 and highlighted in Supplementary Data 2) and measured peptide-specific CD4 + T cell immunity using IFN-γ secretion ELISpot assays with CD8depleted PBMCs derived from three healthy individuals. We detected limited CD4 + immune responses against the peptides tested (Fig. 4b). This prompted us to evaluate the CD4 + T cell immune response against the whole Cas9 protein compared with a positive control protein. We stimulated healthy donor CD8depleted PBMCs with recombinant SpCas9 or EBNA proteins for 10 days and detected a modest response (less than twofold of that of unstimulated cells; Fig. 4c). We then sought to investigate whether the immune response that we detected against peptide β was primarily derived by CD4 + or CD8 + T lymphocytes. We stimulated MACS sorted CD4 + or CD8 + T cells isolated from PBMCs from 3 donors with peptide β and detected a primarily CD8 + response in all 3 donors (Fig. 4d).

Discussion
The detection of pre-existing B cell and T cell immunity to the most widely used nuclease ortholog of the CRISPR-Cas9 tool in a significant proportion of healthy humans confirms previous studies in mice 9,10 and a recent study in humans 27 and sheds light on the need for more studies of the immunological risks of this system. A recent study reported SpCas9-specific T-cell-mediated killing of autologous CRISPR-treated APCs in vitro 28 . However, the consequences of Cas9-specific immunity in vivo remain to be seen. The CD8+T cell immunity we observed is likely memory responses, as they are observed without ex vivo stimulation. However, following 18 days of T cell stimulation by peptides α or β, expansion of naïve T cells is not precluded. This suggests that even in the absence of a pre-existing immune response, the expression of Cas9 in naïve individuals may trigger a T cell response that could prevent subsequent administration. This could be avoided by switching to Cas9 orthologs from other bacterial species, but can be difficult given the epitope conservation across Cas9 proteins from multiple Streptococcus species and resemblance to sequences from other bacterial proteins such as the common pathogen N. meningitidis, that asymptomatically colonizes the nasopharynx in 10% of the population 36 . Therefore, selective deimmunization (immunosilencing) of Cas9 can represent an attractive alternative, particularly in patients where high-dose systemic immunosuppression is contraindicated, such as in patients with chronic infectious diseases. This strategy can be important in particular when longer term expression of Cas9 will be desired.
Using IFN-γ ELISpot, we detected a modest CD4 + immune response against recombinant SpCas9 protein and no response to any of the class II SpCas9 peptides that we used. This is consistent with a recent study that found the majority of the CD4 + response against SpCas9 to be restricted to the Treg compartment with minimum IFN-γ and TNF-α secretion 28 . While a more comprehensive epitope mapping study is needed, we did not experimentally detect immunodominant epitopes among the top binders of the non-HLA-A*02:01 alleles. We show here that silencing one epitope for an HLA-A*02:01 donor (who was also positive for A*11:01, B*39:01, and B*46:01) was sufficient to significantly reduce Cas9 immunogenicity. However, whether other mutations need to be introduced to the protein for complete silencing in HLA-A*02:01 negative individuals needs further investigation. Conventional methods of deimmunizing nonhuman therapeutic proteins rely on trial-and-error mutagenesis or machine learning and often includes deletion of whole regions of the protein [37][38][39][40][41] . Here, as a general principle, we show that alteration of one of the anchor residues of an immunodominant epitope abolished specific T cell recognition. However, HLA allotype diversity and the existence of numerous epitopes in the large Cas9 protein may complicate the process of complete deimmunization. The overall impact of removal of select immunodominant epitopes remains to be seen; both reduction 42 and enhancement 43 of the immunogenicity of subdominant epitopes have been reported with similar approaches for other proteins.
The top binding T cell epitopes within Cas9 that are most promiscuous for common HLA class I and class II alleles have been recently predicted in silico using IEDB 44 . However, this is the first study that experimentally validates predicted immunodominant epitopes. None of the HLA-A*02:01 epitopes that we report overlap with the peptides previously predicted 44 . This is not surprising since our prediction model is optimized for HLA-A*02:01. Thus, improved algorithms are needed to predict epitopes that hold up in experimental validation, as we show here. The use of CRISPR-Cas9 in humans may eventually necessitate creating HLA type-specific Cas9 variants, particularly for applications that require long-term Cas9 expression.
Non-specific localized immune suppressive approaches, such as those used by tumor cells and some viruses may complement these strategies for complete deimmunization. Antigen presentation can be blocked by viral proteins interfering with antigen presentation (VIPRs), such as the adenoviral E319K or US2 and US11 from the human cytomegalovirus 45 or molecules that inhibit proteasomal antigen processing such as the Epstein-Barr virus Gly-Ala repeat 46 . Deimmunized Cas9 may be useful in reduction of the dosage of other immunomodulatory measures needed to be co-administered in patients, thus facilitating therapeutic CRISPR applications as we develop better understanding of the immunological consequences of this system.

Methods
Detection of Cas9-specific serum antibodies. Healthy control sera (n = 143) used in this study, and previously 47 , are a subset of a molecular epidemiology study of head and neck cancer at the MD Anderson Cancer Center, collected between January 2006 and September 2008. Samples were collected using a standardized sample collection protocol and stored at -80°C until use. All relevant ethical regulations for work with human participants were followed and written informed consent was obtained from all participants under the Arizona State University institutional review board approval. Participants were ethnically diverse. Sample   48 , incubated with shaking for 2 h at room temperature, and added to the specified wells in duplicate. Horseradish peroxidase (HRP) anti-human IgG Abs (Jackson ImmunoResearch Laboratories, West Grove, PA) were added at 1:10,000, and detected using Supersignal ELISA Femto Chemiluminescent substrate (Thermo Fisher Scientific, Waltham, MA). Luminescence was detected as relative light units (RLU) on a Glomax 96 Microplate Luminometer (Promega, Madison, WI) at 425 nm. The cutoff value was defined as any reactivity higher than the top 99% of RLU values for human hemoglobin (Fig. 1a, dotted line).
Cas9 candidate T cell epitope prediction. We predicted MHC class I restricted 9mer and 10-mer candidate epitopes derived from the Cas9 protein (Uniprot -Q99ZW2) for HLA A*02:01. The protein reference sequence was entered into 5 different prediction algorithms; 3 MHC-binding: IEDB-consensus binding 49 , NetMHCpan binding 50 , Syfpeithi 51 , and 2 antigen-processing algorithms: IEDBconsensus processing, ANN processing 52 . The individual scores from each of the prediction algorithms were then normalized within the pool of predicted peptides after exclusion of poor binders, and the average normalized binding scores were used to re-rank the candidate peptides. The top 38 candidate peptides (Supplementary Table 1) were selected for experimental testing. The IEDB consensus MHC-binding prediction algorithm (http://www.iedb.org/) was applied to obtain a list of high binding Cas9 peptides, each of which was assigned a normalized binding score (S b ). The immunogenicity score (S i ) was calculated for each peptide based on its amino acid hydrophobicity (ANN-Hydro) 32 . MHC class II epitopes were predicted using the recommended setting on IEDB. We used the Visual Molecular Dynamics software 53 (http://www.ks.uiuc.edu/Research/vmd/) to generate the Cas9 protein structure (PDB ID: 4CMP, Fig. 3a).
Ex vivo stimulation and epitope mapping of Cas9 by ELISpot. All peripheral blood mononuclear cells (PBMCs) were obtained from healthy individuals with written informed consent under the Arizona State University institutional review board and all relevant ethical regulations for work with human participants were followed. Participants were ethnically diverse. Sample size was determined based on prior experience and similar experiments in the literature. There was no group allocation and the investigators were blinded to participants' information. If PBMCs from a given donor were not reactive to the positive control, the donor was excluded from the study. PBMCs were isolated from fresh heparinized blood by Ficoll-Hypaque (GE Healthcare, UK) density gradient centrifugation and stimulated. Briefly, predicted Cas9 peptides with S b < 0.148 (n = 38) were synthesized (>80% purity) by Proimmune, UK. Each peptide was reconstituted at 1 mg/mL in sterile PBS and pools were created by mixing 3-4 candidate peptides. Sterile multiscreen ELISpot plates (Merck Millipore, Billerica, MA, USA) were coated overnight with 5 μg/well of anti-IFN-γ capture antibody (clone D1K, # 3420-3-250, Mabtech, USA) diluted in sterile PBS. Frozen PBMCs were thawed rapidly and recombinant human IL-2 (20 U/mL, R&D Systems) was added. They were then stimulated in triplicates with 10 μg/mL Cas9 peptide pools (or individual peptides), pre-mixed CEF pool as a positive control (ProImmune, UK), or DMSO as a negative control in the anti-IFN-γ-coated ELISpot plates, (Merck Millipore, Billerica, MA, USA) and incubated in a 37°C, 5% CO 2 incubator for 48 h. Plates were washed three times for 5 min each with ELISpot buffer (PBS+0.5% FBS) and incubated with 1 μg/mL anti-IFN-γ secondary detection antibody (clone 7-B6-1, # 3420-6-250, Mabtech, USA) for 2 h at room temperature, washed and incubated with 1 μg/mL Streptavidin ALP conjugate for 1 h at room temperature. The wells were washed again with ELISpot buffer and spots were developed by incubating for 8-10 min with detection buffer (33 μL NBT, 16.5 μL BCIP, in 100 mM Tris-HCl pH 9, 1 mM MgCl 2 , 150 mM NaCl). Plates were left to dry for 2 days and spots were read using the AID ELISpot reader (Autoimmun Diagnostika GmbH, Germany). The average number of spot forming units for each triplicate was calculated for each test peptide or peptide pool and subtracted from the background signal. . APCs were re-stimulated with fresh irradiated KCD40Ls on days 5 and 10, after washing with PBS and expanding into a whole 24-well plate. After two weeks, APC purity was assessed by CD19+ CD86+ expressing cells using flow cytometry, and were used for T cell stimulation after >90% purity. APCs were either restimulated up to 4 weeks or cryopreserved for re-expansion as necessary.
T cell stimulation by autologous APCs. Antigen-specific T cells were detected by stimulating healthy donor B cell APCs by either peptide pulsing of specific Cas9 epitopes, or by transfecting with mRNA encoding the whole WT or modified Cas9 proteins. Peptide pulsing of APCs was done under BCM 5% human serum, with recombinant IL-4. Transfection of APCs was done with primary P3 buffer in a Lonza 4D Nucleofector and program EO117 (Lonza, MD, USA) and incubated in BCM-10% human serum and IL-4. Twenty-four hrs later, on day 1, APCs were washed and incubated with thawed whole PBMCs at a ratio of 1:2 (200,000 APCs: 400,000 PBMCs) in a 24-well plate in BCM supplemented with 20 U/mL recombinant human IL-2 (R&D Systems, MN, USA) and 5 ng/mL IL-7 (R&D Systems, MN, USA). On day 5, partial media exchange was performed by replacing half the well with fresh BCM and IL-2. On day 10, fresh APCs were peptide pulsed in a new 24-well plate. On day 11, expanded T cells were restimulated with peptide-pulsed APCs similar to day 1. T cells were used for T cell assays or immunophenotyped after day 18. Gates for expression of different markers and pentamers were determined based on flow minus one (FMO) samples for each color after doublet discrimination ( Supplementary Fig. 4). Percentages from each of the gated populations were used for the analysis. Modified Cas9 plasmids: human codon-optimized Streptococcus pyogenes Cas9 sequence was amplified from pSpCas9 (pX330; Addgene plasmid ID: 42230), using forward and reverse primers and inserted within gateway entry vectors using golden gate reaction. Desired mutations were designed within gBlocks (Integrated DNA Technologies). The gBlocks and amplicons were then cloned into gateway entry vectors using golden gate reaction. All the primer sequences are listed in Supplementary Table 7 and all the gBlocks sequences are listed in Supplementary Note 1. Next, the Cas9 vectors and CAG promoter cassettes were cloned into an appropriate gateway destination vector via LR reaction (Invitrogen).
Cell culture for endogenous target mutation and activation. HEK293FT cell line was purchased from ATCC and maintained in Dulbecco's modified Eagle's medium (DMEM -Life Technologies) containing 10% fetal bovine serum (FBS -Life Technologies), 2 mM glutamine, 1 mM sodium pyruvate (Life Technologies), and 1% penicillin-streptomycin (Life Technologies) in incubators at 37°C and 5% CO 2 . Polyethylenimine (PEI) was used to transfect HEK293FT cells seeded into 24well plates. Transfection complexes were prepared according to manufacturer's instructions (Polysciences).
Fluorescent reporter assay for quantifying Cas9 function. For the experiment assessing Cas9 cleavage capacity at a synthetic promoter, HEK293FT cells were cotransfected with 200 ng gRNA, 200 ng Cas9 constructs, 50 ng reporter plasmid, and 25 ng enhanced blue florescent protein (EBFP) expressing plasmid as the transfection control. For the experiment assessing Cas9 transcriptional activation capacity at a synthetic promoter, HEK293FT cells were co-transfected with 50 ng gRNA, 70 ng Cas9 constructs, 100 ng MS2-P65-HSF1-GFP (Addgene plasmid ID: 61423), 200 ng reporter plasmid. Fluorescent reporter experiments were performed 48 h after transfection. Nonviable cells were excluded from the analysis using 7-AAD (7-amino-actinomycin D) conjugated with PerCP. Next, we selected cells expressing EBFP >2 × 10 2 A.U. or GFP >2 × 10 2 A.U. (transfection markers) in the cleavage and activation experiments, respectively to exclude untransfected cells ( Supplementary Fig. 4). Flow cytometry was performed using a FACSCelesta flow cytometer (Becton Dickson) with HTS. Flow cytometry data were analyzed using FlowJo software. Untransfected controls were included in each experiment. Experiments underwent initial validation in duplicate and then were repeated in triplicate for the final manuscript.
Quantitative RT-PCR analysis. HEK293FT cells were co-transfected with 10 ng gRNA, 200 ng Cas9 constructs, 100 ng MS2-P65-HSF1 (Addgene plasmid ID: 61423) and 25 ng EBFP plasmid as the transfection control. Cells were lysed, and RNA was extracted using RNeasy Plus mini kit (Qiagen) 72 h post transfection, followed by cDNA synthesis using the High-Capacity RNA-to-cDNA Kit (Thermo Fisher). qRT-PCR was performed using SYBR Green PCR Master Mix (Thermo Fisher). All analyses were normalized to 18S rRNA (ΔCt) and fold changes were calculated against untransfected controls (2 −ΔΔCt ). Primer sequences for qPCR are listed in Supplementary Table 7.
Endogenous indel analysis. HEK293FT cells were co-transfected with 200 ng of Cas9 plasmids, 10 ng of gRNA coding cassette, and 25 ng EBFP plasmid as the transfection control. 72 h later, transfected cells were dissociated and spun down at 200 g for 5 min at room temperature. Genomic DNA was extracted using 50 µl of QuickExtract DNA extraction solution (Epicentre) according to the manufacturer's instructions. Genomic DNA was amplified by PCR using primers flanking the targeted region. Illumina Tru-Seq library was created by ligating partial adapters and a unique barcode to the DNA samples. Next, a small number of PCR cycles was performed to complete the partial adapters. Equal amounts of each sample were then pooled and sequenced on Illumina Tru-Seq platform with 2 × 150 run parameters, which yielded approximately 80,000 reads per sample. Sequencing was performed using a 2 × 150 paired-end (PE) configuration by CCIB DNA Core Facility at Massachusetts General Hospital (Cambridge, MA, USA). The reads were aligned to the target gene reference in Mus musculus genome using Geneious software, 9-1-5. To detect the indels (insertions and deletions of nucleic acid sequence at the site of double-strand break), each mutation was evaluated carefully in order to exclude the ones that are caused by sequencing error or any off-target mutation. The variant frequencies (percentage to total) assigned to each read containing indels were summed up. i.e. indel percentage = total number of indel containing reads/ total number of reads. The minimum number of analyzed reads per sample was 70,000.
RNA sequencing for quantifying activator specificity. HEK293FT cells were cotransfected with 10 ng gRNA targeting MIAT locus, 200 ng Cas9 constructs, 100 ng MS2-P65-HSF1 (Addgene plasmid ID: 61423), and 25 ng transfection control. Total RNA was extracted 72 h post transfection using RNeasy Plus mini kit (Qiagen) and sent to UCLA TCGB core on dry ice. Ribosomal RNA depletion, and single read library preparation were performed at UCLA core followed by RNA sequencing using NextSeq500. Coverage was 14 million reads per sample. FASTQ files with single-ended 75 bp reads were then aligned to the human GRCh38 reference genome sequence (Ensembl release 90) with STAR 54 , and uniquelymapped read counts (an average of 14.8 million reads per sample) were obtained with Cufflink 55 . The read counts for each sample were then normalized for the library size to CPM (counts per million reads) with edgeR 56 . Custom R scripts were then used to generate plots.
Reporting summary. Further information on experimental design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The datasets generated and/or analyzed during the current study have been submitted to the Gene Expression Omnibus (GEO) repository and are accessible through GEO series accession number GSE125522 . Source data for this manuscript are available in the Source Data file. All other relevant data are available from the authors upon reasonable request.


Benchmarking GABA Quantification: A Ground Truth Data Set and Comparative Analysis of TARQUIN, LCModel, jMRUI and Gannet

Many tools exist for the quantification of GABA-edited magnetic resonance spectroscopy (MRS) data. Despite a recent consensus effort by the MRS community, literature comparing them is sparse but indicates a methodological bias. While invivo data sets can ascertain the level of agreement between tools, ground-truth is required to establish accuracy, and investigate the sources of discrepancy. We present a novel approach to benchmarking GABA quantification tools, using several series of phantom experiments with iterated GABA concentration. Each series presents a different set of background metabolites and environmental conditions allowing comparison of not only individual estimates, but the ability of tools to characterise changes in GABA across a range of potential confounds. The methodology of the phantom experiments is presented, as well as characterisation of the data. We also perform an initial comparative analysis of several common MRS quantification tools and in-house code, to illustrate utility of the dataset and the potential bias introduced by different quantification methods. The GABA-to-NAA ratios reported by each tool are compared to the ground-truth, and estimation accuracy is assessed by linear regression of this relationship. While the linearity of GABA-to-NAA gradients is generally captured by all of the tools, a large variation in the slope, offset and environmental stability of the gradient is observed. The primary driver of differences in linear combination modelling is the choice of basis function. However, tools employing a common basis and pre-processing still produce differences on the order of 4%. Less-strictly parametrised fitting approaches appear to improve the robustness of quantification, but accurate modelling of edit efficiency calculations is still necessary to avoid systematic offsets. In general, the level of variation suggests that comparisons across...

INTRODUCTION
Gamma-aminobutyric acid (GABA) is the primary inhibitory neurotransmitter in the central nervous system, maintaining the brain's excitation-inhibition balance. Its prominent role in neurotransmission and metabolism has led to extensive study and a plethora of applications for its detection. Disruption of GABAergic processes has been observed in Schizophrenic patients 1,2,3 and GABA receptor dysfunction has associations with epilepsy 4,5 . GABA processes are prominent in type I diabetes 6 and autism spectrum disorders 7 . GABA is also the subject of intense study by the psychology community, with GABA levels influencing impulsivity 8,9 , drug addiction 10 , anxiety disorders 11 and depression 12 . The most common method for non-invasive detection of GABA in-vivo is magnetic resonance spectroscopy (MRS). However, at clinically relevant field strengths the GABA signal is overlapped by more prevalent metabolites, necessitating other detection methods for reliable quantification.
While specialised techniques such as 2D J-spectroscopy 13 are used, the most common method is edited MRS 14,15,16,17 , in particular MEGA-PRESS 18,16,19,20 . MEGA-PRESS GABA-editing utilises two interleaved acquisitions -'edit on' and 'edit off' -which differ only in the frequency of applied editing pulses. The edit frequencies are chosen to take advantage of the J-evolution of the target metabolites, such that upon subtraction of the two acquisitions, only resonances within the edit bandwidth and those coupled to them remain. The widespread use of edited MRS has seen the development of a range of software packages aimed at processing and quantification of the spectra 21,22,23,24,25,26,27 . The tools process data by some combination of apodization, frequency and phase calibration, macromolecular baseline subtraction and residual water signal removal.
Each step may be realised by a diverse array of methods, providing many potential sources of variation for the quantified signal. This effect is compounded by further differences in the quantification procedure. The most common approach to quantification is the generation of basis spectra -via phantom scans or more commonly by quantum mechanical simulation -followed by the application of a fitting algorithm to decompose the MRS signal into weighted combinations of the basis functions. For edited MRS, with its simplified spectra, other fitting approaches -using Gaussian/Lorentzian modelling of resonances and numerical integration -become more feasible.
A recent large-scale effort has been made to establish consensus in the MRS community 28,29,30,31,32,33,34 . While recommendations were made around the quantification of MRS data -Choi et al. 29 reference basis set methods and Gaussian peak fitting for GABA quantification -there was no consensus on the best approach. While it has been previously argued that different methods, suitably applied, exhibit similar variation and clinical observations 35,36 , it is acknowledged in the literature that there is a procedural dependency in MRS analysis 37,38,39,40,41 , which provides a statistically significant variation in reported metabolite amplitudes or relative concentrations 42,43 . With the exact pre-processing, simulation and fitting parameters only partially reported in many research articles -an issue discussed by Lin et al. 30 -the methodological dependence of MRS quantification propagates, reducing the reproducibility of MRS studies. Despite this, there has been relatively little comparative study and benchmarking of the MRS analysis tools.
The goal of this study was to generate benchmarking data sets for the purpose of assessing the relative performance of GABA quantification tools and investigate the methodological dependence of GABA-edited data analysis. Previous work in this area has utilised either in-vivo data -for which there is usually no ground truth to evaluate the results -or simulated data sets, which inevitably cannot capture all sources of experimental variation. For this study, MRS benchmarking data sets were generated using several series of calibrated phantom experiments. While phantoms fail to capture all influences of the in-vivo environment, they do allow the inclusion of many experimental factors, while still enabling a ground truth comparison and ability to control the environmental characteristics expected to influence results. Four series of phantom data were experimentally acquired. In each case, the concentrations of various metabolites -creatine (Cr), N-acetyl-L-aspartic acid (NAA), glutamine (Gln) and glutamate (Glu) -are fixed at approximately in-vivo levels, while the concentration of GABA is iteratively increased. At each GABA concentration step, MEGA-PRESS spectra were acquired. The four series of experiments were designed to increase the complexity of the spectra, starting with pH and temperature calibrated solutions containing only NAA, Cr and GABA (E1), then adding Glu and Gln to the solutions to increase the complexity of the spectra (E3), and finally progressing to a series of gel phantoms which more closely mimic in-vivo spectra 44 (E4 1 ). To investigate the influence of the molecular environment, a series of non-neutral pH solutions of NAA, Cr and GABA (E2) were also acquired. The phantom data set has been designed with the benchmarking of new quantification tools in mind, for example in the use of convolutional neural networks for GABA quantification 45 . As such, the full data have been made available for others to benchmark their own analysis software, or optimise existing tools (DOI:10.21227/ak1d-3s20) 46 .
Furthermore, to demonstrate the utility of the dataset, MEGAPRESS spectra were quantified using a variety of state-of-the-art MRS analysis software to investigate the potential methodological bias that could be introduced. The GABA-to-NAA concentration ratio was estimated for each tool by applying appropriate signal attenuation corrections where needed. The reported ratios are plotted against the known ground-truth, and fit using linear regression. The GABA-to-NAA gradients are used to investigate the systematic offsets that can be introduced into reported concentration changes, and to inform future comparative benchmarking studies.

METHODS
Phantom preparation procedures, experimental setup and scan parameters, and details of the analysis are described.

Phantom preparation
The  Table 1 and representative images of the phantoms are shown in Figure A1.
The solution series were prepared by dissolving the required concentration of background metabolites in de-ionised water to create a base solution. 290 ml of the GABA-free base solution was transferred to a round bottom glass flask, filling it to the neck. To avoid dilution of the background metabolites, a small amount of this base solution was removed from the flask using a syringe and GABA was added to prepare a concentrated GABA solution (1 mg/ml). To vary the GABA concentration between scans, a small amount of the solution was removed from the flask with a syringe and replaced by the same volume of the concentrated GABA solution. The amount of solution replaced in each step varied depending on the desired GABA concentration change, but was generally between 0.5 ml and 2 ml. This procedure was repeated several times, iterating GABA to a relatively high concentration with the same scan protocols applied for each concentration step. This procedure limited experimental errors solely to the GABA concentration. The pH of the solution was monitored between scans and adjustments were made using a 36% hydrochloric acid (CAS-7647-01-0) solution and a 3.99% sodium hydroxide solution (CAS-1310-73-2) to maintain a pH of 7.2 ± 0.2. The amounts of HCl or NaOH added after the initial pH calibration of the base solution were marginal (≤ 1 ml for the 290 ml flask) and any dilution of the overall solution was deemed negligible. The E2 series was one exception to this procedure, where the pH was intentionally not adjusted, to investigate the effect this had on the various quantification procedures. The final pH of E2 was measured to be 3.0 ± 0.2.
For the gel series, the phantoms were made in advance of scanning. 800 ml of a base solution of Glu, Gln, NAA and Cr was prepared and divided into eight 100 ml portions. Different amounts of GABA were added to each solution and pH calibration was performed as above. Finally, 1 g of agar (CAS-9002-18-0) 2 was added as gelling agent. The mixtures were then heated to 90°C to 100°C while being stirred until the agar had fully dissolved. A small hole (< 2 mm diameter) was created in a spherical plastic mould using a plastic welding tool and the solutions were injected via a syringe. The arrangement was then allowed to cool over night. Once solidified, the gels were examined and the opening was sealed using a small amount of silicone sealant. Spherically-shaped phantoms were deemed the most suitable to reduce magnetic susceptibility-induced field inhomogeneity and minimise the spectral linewidth.
The accuracy of the scale used for mass measurements is 1 mg and the graduated cylinders used for the volume measurements have an accuracy of 1 ml. Assuming negligible errors in the manufacturer's stated molar masses, the uncertainty ∆c in the molar concentration c = m M V , where m is the mass in grams, M is the molar mass of the solute and V is the volume of the solvent, is ∆c c = ∆m m 2 + ∆V V 2 . Assuming (generously) an uncertainty of 2% for the mass measurements and 1% for the volume measurement, the uncertainty in the concentrations of base metabolites is estimated to be less than 3%. The uncertainty in the amount of solution added or removed from the flask in each step via 5 ml syringes is slightly higher and could be up to 10%, assuming an uncertainty of 0.1 ml per 1 ml to account for slight misreadings of the plunger position and tiny air pockets within the syringe.

Scan Protocols and Data Acquisition
All MR scans were conducted at Swansea University's Clinical Imaging Facility using a MAGNETOM Skyra 3 T (Siemens Healthcare GmBH, Erlangen, Germany). The scanner room was temperature controlled to 20 ±0 6°C. Between each scan, the phantom was removed to update concentrations, so the influence of heating is deemed negligible. Signal acquisition was done using the built-in spine coils, specifically the four channel spine coil element 'SP2'. The spine coils provided the platform for the most reproducible set-up and SP2 exhibited the highest SNR of the coil elements available. The phantom was aligned with this element and raised to isocenter, using a cardboard phantom holder to ensure maximum field homogeneity for the scans. Double-spin-echo field maps were acquired for each phantom to assess the homogeneity, and manual shimming and frequency calibration was performed to optimise the spectral width. Further manual calibration was performed to adjust the TX reference voltage and water-suppression pulses.
A 20 × 20 × 20mm 3 voxel at the isocenter of the magnet was selected for single voxel spectroscopy. The GABA-edited spectra were acquired using the MEGA-PRESS pulse sequence WIP859G 3 with CHESS water-suppression 47  For each experiment, the raw data acquired, the single-average coil-combined spectra and combined average edit-on, edit-off and difference spectra produced by the vendor-supplied spectroscopy software have been made available 46 . For the purpose of this paper, the combined channel and average dicom spectra were used. While raw data is generally preferable in order to appraise frequency and phase drift correction, the averaged dicom data was most universally supported by all the software packages evaluated, allowing a more encompassing cross-section of the methods.
All spectra used in this study were inspected for reconstruction issues, prior to their inclusion.

Spectral quality
Prior to quantification, the 2 ppm NAA and water peaks in each spectrum are fitted by Lorentzians, from which the position, full-width-halfmaximum (FWHM) and maximum signal amplitude of the resonances are obtained. The noise floor was estimated by computation of the standard deviation, σ, of the signal in the region 8−9 ppm, observed to be free of metabolite signal, then the signal-to-noise ratio (SNR) calculated: S 2σ . Water suppression factors are estimated by comparing the height and area of the water peak in the water unsuppressed and water-suppressed spectra.

Analysis methods
Analysis methods can be broadly separated into two categories: (i) basis set methods (TARQUIN 27 , jMRUI 22,26 , LCModel 24 ), which attempt to fit the data using a set of simulated or experimentally obtained basis spectra and (ii) peak fitting and integration tools (Gannet 21 , LWFIT -An in-house peak integration method used as a baseline in this study), which compute the area of peaks associated with various metabolites. As the absolute amplitudes are arbitrary and depend on the scanner design, transmitter calibration, coils used and the phantom itself. We focus on quantifying the relative contributions of metabolite signals, specifically the GABA-to-NAA ratio. The concentration ratio estimate, C, between a target metabolite, M , and reference, R, is given by: where S is the reported resonance amplitude, N is the number of contributing spins, R is the signal attenuation factor resulting from T1 and T2 relaxation, and E is the edit efficiency. For the purposes of our analysis, the difference between the relaxation properties of GABA and NAA are assumed negligible, as suggested by Choi et al. 29 . With this assumption, basis set methods report appropriately scaled metabolite amplitudes.
However, the peak fitting methods require corrections for N and E. Both Gannet and LWFIT quantify the 3 ppm GABA resonance and 2 ppm NAA singlet, resulting in N R N M = 3/2. By default, Gannet's edit efficiency correction factor is estimated as E R /E M = 0.5, and LWFIT's correction factor was calculated from simulation (Section 3.2). The concentration ratios reported by each tool were plotted vs. the known ground-truth ratio, and the relationship calculated via linear regression.
For the purposes of this study, the tools were employed as close to "out-of-the-box" as possible. Default settings were used with additional developer recommendations whenever these were readily available. While it is indeed possible to further optimise the parameters of each tool, our aim was to examine the potential methodological variation that could occur in GABA MRS studies. A brief description of each tool is provided below, and readers are directed to the referenced literature for further information.

Gannet
GABA-MRS analysis tool (Gannet) 21 is an open-source automated MATLAB-based analysis tool, specifically designed for automated processing of 3 T GABA MEGA-PRESS data. Gannet is able to process raw data, further reducing user dependence, performing channel combination, adding line broadening, frequency and phase corrections, outlier rejection and time averaging. Fitting is performed using non-linear least-squares algorithms (lsqcurvefit and nlinfi) with a five parameter Gaussian model to estimate the 3 ppm GABA signal in the difference spectrum, a six parameter Lorentzian model to estimate the 3 ppm Cr signal in the off spectrum and, if appropriate, a six parameter Gaussian-Lorentzian model to fit the unsuppressed water signal. Quantitative results are reported as integral ratios of GABA to Cr and a concentration relative to water, NAA or Glx (Glu and Gln combined). While Gannet attempts to remove user interaction from the analysis procedure, it allows user modification of the code and adjustment of the pre-initialisation script is encouraged depending on setup. For this paper Gannet version 3.0 was used and the pre-initialisation script was adjusted to account for room temperature phantom acquisition.

jMRUI
jMRUI 22,26 is proprietary software distributed under its own license terms and made freely available to registered users for non-commercial use.
It has a large suite of pre-processing, analysis and simulation options. For this study, two basis set methods were used to quantify the spectra.
The first, quantitation based on semi-parametric quantum estimation (QUEST) 48,49 , is a time-domain fitting tool, which includes a semi-parametric approach to handle spurious signals resulting from macro molecules and lipids. The second method, automated quantitation of short echo time MRS spectra (AQSES) 50 , uses a modified VARPRO variable projection NLLS algorithm 51 , which imposes prior knowledge in the form of upper and lower bounds on the non-linear parameters.
For this study, jMRUI version 6.0 beta was used. MEGA-PRESS basis sets were generated using NMR-SCOPE-B 52,53 by defining a single MEGA-PRESS pulse sequence with sequence parameters matching our experimental protocol. Two basis sets for GABA/Cr/NAA and GABA/Cr/-NAA/Glu/Gln were generated. Manual frequency and phase calibration and apodization were performed on the input spectra and fitting with both tools was applied using this common input. This allowed investigation of the influence of the fit algorithm itself, independent of pre-processing.

LCModel
LCModel 23,24 is a formerly commercial, Linux-based MRS analysis tool which recently moved to open-source. It can be used with its internal in-vitro basis set or any arbitrary basis set specified by the user. The baseline signal resulting from macromolecules is established using spline fits. Fitting is performed using a Gauss-Newton non-linear least-squares algorithm with Marquardt modification 54 with additional terms for T2, spectral shift and field inhomogeneities. For this paper, version 6.3-1L and the Purcell lab basis set 55 recommended by LCModel's creators were used.

LWFIT
Our in-house code 56 , LWFIT, written in MATLAB, aligns the edit-off and difference spectra so that the main NAA peak is precisely at 2.01 ppm and then calculates the areas for the NAA (1.91 ppm to 2.11 ppm, difference), Cr (2.90 ppm to 3.10 ppm, edit-off) and GABA (2.90 ppm to 3.12 ppm, difference) peaks, as well the two main peaks associated with Glu and Gln, GLX1 (2.25 ppm to 2.45 ppm, difference) and GLX2 (3.65 ppm to 3.85 ppm, difference), using the real frequency-domain spectra. This is done by numerical integration using piecewise-linear functions over the indicated fixed ppm ranges, selected to minimise contamination from other signals. Lorentzian, Gaussian and spline fitting and several baseline fitting and filtering methods were explored but initial testing indicated that numerical integration with minimal pre-processing yielded the most accurate estimation results of the methods considered. While more aggressive noise filters aesthetically improve the quality of the spectra, their application tends to exacerbate underestimation of GABA, especially for weak signals, suggesting that filtering eliminates some of the GABA signal present in the data 57 . For this work, only zero-filling (N = 4096 points) and a frequency shift to align the NAA peak in the spectrum to 2.01 ppm are performed.

Simulations
To further facilitate interpretation of the experimental results, quantum mechanical simulations were also performed. Simulation is fundamental in MRS quantification. Basis set fitting methods utilise simulated spectra to model the peak structure to which the experimental data are fit. While peak fitting methods are less tightly parametrised during the fit, they still must compensate for edit efficiency -the proportion of peak amplitude lost as a part of the editing process -in order to accurately quantify the metabolites. The simulations were used to investigate the influence of ideal RF-pulse simulation and GABA models on quantification, and to facilitate edit efficiency correction for LWFIT.
The experimental setup was emulated, modulating GABA concentration while maintaining a fixed concentration of other selected metabolites.
MEGA-PRESS simulations were performed using FID-A 61 , a simulation and data processing package for MRS. MEGA-PRESS difference spectra were generated, mirroring the sequence parameters of the phantom study with a finite-bandwidth editing pulse alternating between 1.9 ppm respectively. Simulated series were tested with all three GABA models to investigate any potential bias of our simulated data. For all simulated series, GABA concentrations were varied between 1 mM and 12 mM at 1 mM intervals and the combined spectra were analysed using our in-house code.

RESULTS AND ANALYSIS
We initially present the experimentally acquired spectra along with the corresponding calibration and quality measures. We then present the simulation and quantification results.

Calibration and characterisation of spectra
The stability of the experimental setup over the course of an entire series of experiments was assessed by considering the variation of core parameters such as transmitter voltage and frequency, shim gradients and water suppression settings. Table 2 shows that a small adjustment made to the experimental procedure for E3 and E4 to avoid table movements between successive scans significantly reduced adjustments of transmitter frequency, reference voltage and linear gradient offsets Gx, Gy, Gz, although all results are within acceptable limits.
The difference spectra for all four experimental series are shown in Figure 1. The expected pseudo-triplet structure, clearly visible for the wellcalibrated solutions E1 and E3, is still discernible in the gel phantom spectra E4, despite the broader linewidth. However, this structure is not present in E2 where we observe a broad, almost Gaussian resonance. Despite this, the highlighted area does appear to increase linearly with GABA concentration. Table 3 gives the linewidths of the water and NAA peaks obtained from Lorentzian peak fits, indicative values for the water suppression factors (WSF) and the SNR of the NAA peak. While the linewidths increase for the gel phantoms as expected, the width of the NAA peak in the difference spectrum remains around 4 Hz, suggesting good field homogeneity and frequency stability over the duration of the experiments. SNR over 200 for NAA in the solution series and over 100 for the gel phantoms are also well within the acceptable quality limits.

(a) Spectra for series E1 (b) Spectra for series E2
(c) Spectra for series E3 (d) Spectra for gel phantom series E4

FIGURE 1
Difference spectra aligned to 2.01 ppm NAA peak. The 3 ppm GABA peak is shaded according to LWFIT integration limits, and vertically offset spectra correspond to increasing concentrations of GABA for the phantoms.    Figure 2 shows the GABA-to-NAA amplitude ratio obtained by LWFIT vs. the actual imposed GABA-to-NAA ratio for three series of simulated experiments. Simulations using ideal excitation and refocusing pulses produce a gradient of 1.01 -within 1% of the actual ratio -an intercept of −0.05 and an R 2 value of 1.00 (up to 3 significant digits). However, for spectra simulated with finite-duration refocusing pulses with realistic pulse envelopes, the GABA-to-NAA ratio gradient is only 84% of its expected value for the NAA-GABA-Cre series and this is further reduced to 80% when Glu and Gln are added. Figure 3 shows the 3 ppm GABA peak for the edit-on, edit-off and difference spectra as a function of position on a 2D grid. The difference spectra near the boundary of the selected voxel indicate that the editing efficiency is significantly reduced. In an experimental setting, the editing efficiency will be further reduced by imperfect slice profiles of the excitation pulses and other factors such as field inhomogeneity. Assuming these issues lead to a further reduction of the editing efficiency by a factor of 0.84 and 0.80, respectively, the simulations  Edit efficiency is a crucial factor in edited MRS quantification. Basis set fitting methods should have this incorporated into the basis by appropriate simulation of the basis functions. Peak fitting methods, on the other hand, require an ad-hoc factor to compensate for the resulting reduction in peak amplitude. To ensure the best GABA estimates, full simulations with correct pulse envelopes are required to fully capture the J-modulated peak structure for basis set methods, and editing efficiency factor for peak fitting.

Quantification
The quantification results obtained by the various fitting approaches are summarised in Table 4. The results are reported as gradient, a, and intercept, b, with 95% confidence intervals and R 2 value of the linear fit. The concentration ratios are also plotted against the known ground-truth ratios in   Linear fit results for GABA-to-NAA ratios in all four experimental series. Reported are the gradient a and intercept b of the linear fit with 95% confidence intervals (CI) and the coefficient of determination R 2 (to two significant digits). Perfect quantification corresponds to a = 1, b = 0 and R 2 = 1. The percentage variation of reported amplitudes for NAA, whose actual concentration is constant, is also given. The large % variation of NAA amplitudes reported by TARQUIN for E1 is due to an unexplained rescaling of metabolite amplitudes between scans 10 and 11. The % variation of NAA for scans 1 -10 is 6.09%. There is no simple explanation for the large variation reported by TARQUIN for E4.  Figure 4b shows results for the intentionally miscalibrated solution series, E2. Unlike the other series, E2 did not exhibit the characteristic pseudo-triplet GABA resonance, and instead presented a single broad peak at 3 ppm. While a pH of 3.0 is beyond what is commonly observed invivo, this series was intended as a means to study the effect of chemical environmental changes upon quantification, a relevant consideration when investigating certain pathologies. Surprisingly, both TARQUIN and jMRUI give significantly more accurate GABA-to-NAA gradients, with a close to 1 and higher R 2 than for the pH-calibrated series. The TARQUIN change appears to be a consequence of its previous tendency to over-estimate GABA, while jMRUI's GABA gradient increase is likely driven by a reduced efficacy in NAA quantification. Peak fitting tools were expected to be more robust to environmental pH changes. With less tightly constrained models, the majority of signal is still captured by the tool, leaving only a systematic offset resulting from edit efficiency correction. We found that LWFIT and Gannet maintain a high R 2 , with LWFIT showing a reduced GABA-to-NAA gradient, likely a result of integration limits cutting off part of the broad resonance. Gannet, on the other hand, exhibits a modest increase, however there is significant overlap with the E1 confidence interval, so this finding may not be significant. Gannet's Gaussian model of GABA seems a suitable choice for this series, and the robustness of Gannet's estimations support this.  Table 2. Figure 4d shows the results for the gel phantom series with Glu and Gln. The fit quality for the gel series is generally lower, as expected due to higher FWHM and lower SNR for tissue-mimicking gels. The mean R 2 of E4 GABA-to-NAA gradient estimates across all tools is 0.65 for E4,

DISCUSSION
The main aim of this work was to develop an iterative phantom-based method for benchmarking GABA quantification software. The phantom data exhibit a high SNR with narrow metabolite line-widths, well above quality thresholds outlined in consensus work. The four experimental series provide unique challenges to quantification software, allowing elucidation of several confounds which would not be possible using in-vivo data, alone. E1 represents the simplest case with only GABA, NAA and Cr, E2 introduces environmental changes to pH, E3 adds Glx to investigate its potential to interfere in GABA quantification, and the gels of E4 present broader and lower SNR spectra. LWFIT -intended to serve as a baseline sanity check -reproduces the expected gradient to a high fidelity for E1, E3 and E4. LWFIT does, however, underestimate the GABA-to-NAA concentration gradient for E2, with its pH offset. The broad 3 ppm GABA resonance present in E2 is likely insufficiently captured by the integration limits used by LWFIT, furthermore, edit efficiency calculations are likely invalidated somewhat by the apparent structural changes in the spectra.
While all tools generally estimated a reasonable linearity in the GABA-to-NAA gradients, there is a systematic offset in these estimates. For TARQUIN, jMRUI and LCModel, these changes are primarily driven by the choice of basis set. This can be seen by comparing the two jMRUI methods -QUEST and AQSES -which share a common basis. Their mean estimates of GABA-to-NAA gradients differed by just 4% across all the experimental series. The discrepancy increases to around 106% when comparing TARQUIN and LCModel. The basis issue is further complicated by differences in the default GABA models used by the various tools. LCModel includes the Govindaraju and Kaiser model, while jMRUI uses its own modified version of the Govindaraju and Kaiser models. GABA spectra simulated using the Govindaraju, Kaiser and Near model, shown in Figure A2 exhibit clear differences in the structure of the 3 ppm GABA peak. Generally, the Kaiser and Near models appear to approximate the experimental spectra better than the older Govindaraju model, in agreement with observations in the literature. However, despite this, many tools still maintain the older Govidaraju GABA model by default. It is therefore not possible to fundamentally ascertain the performance of these tools without a systematic study with common basis sets and optimised settings. An interesting follow-up study would be to generate a common optimised basis set to be used in several basis set fitting methods. The phantom series can be further enriched by adding degrees of noise, linebroadening, experimentally acquired macromolecular baselines and other co-edited resonances to investigate their influence on quantification, providing many avenues of further study.
With the exception of TARQUIN, the tools present a tendency to underestimate GABA-to-NAA ratios to varying degrees. Simulations suggest realistic refocusing pulses reduce the editing efficiency and can lead to a reduction in reported GABA-to-NAA ratios of up to 20%. Including additional losses in experiments due to line broadening as a result of imperfect shimming, and excitation pulses with realistic slice profiles could therefore explain the underestimation of GABA-to-NAA gradients by many of the tools discussed. These losses can be taken into account by appropriate simulation during the creation of basis functions, or as a part of the edit efficiency correction for peak fitting tools. Table 5 contains a meta analysis of the reported GABA-to-NAA gradient values, specifically the mean gradient estimate and standard deviation of this value across all experimental series. The mean provides an indication of each tool's tendency to under or over estimate the changes in GABAto-NAA. Whereas the standard deviation, σ, provides an indication of the robustness of each tool's estimations to changes in the background conditions. A notable observation is that the peak fitting methods -Gannet and LWFIT -have far lower σ than all basis set methods. This suggests that the less strictly parametrised approaches might be preferable for edited spectroscopy analysis, where the landscape is greatly simplified compared to conventional MRS. This is in agreement with the observations of Marzola et al. 43 , who also reported the robustness of these methods in-vitro but note the poor performance of integration methods in-vivo. However, it should be noted that when series E2 is excluded, LCModel exhibits a similar results to Gannet, with a mean and standard deviation of 0.5583 and 0.0939, respectively.
The fit error is one aspect that was not approached in this study. While basis set methods generally report the Cramer-Rao lower-bound (CRLB) as the preferred measure of fit quality for a particular basis signal, this functionality is not present in Gannet and not valid for LWFIT, limiting comparability. Furthermore, the use of CRLB as a quality filter has been shown to cause bias in mean concentration estimates 66 . Fit residuals could provide some insight into fit accuracy, however these will be modulated by baseline approximation procedures, which differ substantially between tools, making comparisons potentially unfair. While efforts were made to ensure reasonable fit quality in this study, future comparative work will need to establish a means of common error reporting, using the same, sensibly justified metric. This is a fundamental consideration of quantification benchmarking.

CONCLUSION
The phantom benchmark dataset presented here provided a novel, interesting way to compare GABA MRS analysis software. The dataset is intended to form part of a larger systematic benchmarking of analysis tools, a necessary future step for the MRS community. While indeed preliminary, the initial comparison of software seems to indicate that a less-strictly parametrised peak fitting methods might be the best option for edited MRS, where the simplified spectral landscape does not necessitate complete modelling of the peak. This is, however, still a necessary step, as modelling of edit efficiency is important in reducing systematic offsets in this approach. The methodological dependence of the quantification results observed also suggests that care must be taken when comparing results across studies, where analysis pipelines may differ and under-reporting of fitting methods muddy the waters for study replication.


An empirical assessment of research practices across 163 clinical trials of tumor-bearing companion dogs

Comparative clinical trials of domestic dogs with spontaneously-occurring cancers are increasingly common. Canine cancers are likely more representative of human cancers than induced murine tumors. These trials could bridge murine models and human trials and better prioritize drug candidates. Such investigations also benefit veterinary patients. We aimed to evaluate the design and reporting practices of clinical trials containing ≥2 arms and involving tumor-bearing dogs. 163 trials containing 8552 animals were systematically retrieved from PubMed (searched 1/18/18). Data extracted included sample sizes, response criteria, study design, and outcome reporting. Low sample sizes were prevalent (median n = 33). The median detectable hazard ratio was 0.3 for overall survival and 0.06 for disease progression. Progressive disease thresholds for studies that did not adopt VCOG-RECIST guidelines varied in stringency. Additionally, there was significant underreporting across all Cochrane risk of bias categories. The proportion of studies with unclear reporting ranged from 44% (randomization) to 94% (selective reporting). 72% of studies also failed to define a primary outcome. The present study confirms previous findings that clinical trials in dogs need to be improved, particularly regarding low statistical power and underreporting of design and outcomes.

Regulatory status of treatment interventions.
For clinical trials evaluating treatment efficacy of small molecules (drugs) and biologics, we determined whether the anticancer treatments being evaluated had been licensed by the FDA for use in human patients prior to the veterinary trial. For this purpose, one investigator (YJT) searched Wikipedia and the Drugs@FDA database, as needed, to determine approval status as of December 2018; and to determine the length of time between approval and publication of the canine trial. When unlicensed, YJT searched PubMed and PubChem to find out if any human study using the same drug or biologic had been published before the publication of the respective trial in dogs. In addition, the reference list for included veterinary trials studying unlicensed interventions were searched to determine if they cited any human studies. When tested interventions in veterinary trials included combination of two or more drugs and/or biologic, we examined each one of them separately.
Sample size and power calculation. For each study, we manually extracted the sample size across the entire study, the sample size per American Cancer Society tumor type 20 tested per study, and the sample size per intervention arm. In addition, we manually extracted the number of events per arm for both survival and disease progression analyses. The number of events for survival analyses was defined as the number of deaths from any cause, including euthanasia, by the end of the study. For disease progression, this was defined as the number of dogs with disease progression at the end of the study. Only studies that defined progression as tumor growth by a set percentage were included in this analysis. Deaths without disease progression were not considered as a disease progression event. The detectable hazard ratio given the ratio of subjects per arm and number of events was calculated using the Schoenfeld formula 21 , assuming a Type I error rate of 0.05 and power of 0.8. A single assessor (YJT) performed data extraction and analysis.
Disease progression criteria extraction. To determine if the criteria for therapy response in clinical trials in tumor-bearing dogs had changed over time, we determined if the VCOG RECIST criteria was explicitly named as the criterion used, and extracted the minimum change in tumor size for partial response and progressive disease as defined by the authors in each study. Metrics for size as defined by the authors, such as volume or diameter, were also extracted. A single assessor (YJT) performed data extraction.
Cochrane risk of bias assessment. Each study underwent a Cochrane risk of bias assessment 22 . In brief, studies were judged for their risk of bias in six categories according to guidance provided for the Cochrane tool: random sequence generation, allocation concealment, blinded participants/personnel, blinded assessment of results, incomplete data, and selective reporting (i.e. trial registration). Studies were assessed to be at high risk of bias, low risk of bias, or to have unclear risk due to insufficient information. Owing to the subjective nature of this assessment, two assessors (YJT and RJC) performed risk of bias analyses independently. JPAI adjudicated any www.nature.com/scientificreports www.nature.com/scientificreports/ discrepancies between the two authors' assessments that could not be resolved. The consensus risk of bias assessment for all 163 studies agreed upon by the two assessors is presented in this manuscript.
Outcome reporting. Outcomes reported as primary outcomes within each paper, as well as outcomes listed in the abstract of each paper were assessed for statistical significance at the level defined by the authors (typically p < 0.05). Studies were assessed for the presence of a defined primary efficacy outcome, a reported abstract efficacy outcome, and ≥1 statistically significant primary or abstract efficacy outcomes. Studies were also assessed for the presence of significant outcomes showing that the experimental active intervention was harmful compared with standard of care or placebo. Studies that did not conduct statistical tests on primary or abstract outcomes were also identified. Untested outcomes were further characterized according to the authors' qualitative assessments of whether they were favorable, i.e. demonstrated an efficacy benefit versus standard of care or placebo control. A single assessor (YJT) extracted outcome data. Data cleaning. Raw data was manually extracted from publications before undergoing a data cleaning and analysis process using R to arrive at data in its presented form. The following R packages were used for cleaning, analysis, and plotting: readxl 23 , stringr 24 , tidyr 25 , splitstackshape 26 , plyr 27 , RISmed 28 , ggplot2 29 . Additionally, the RoB Summary function was used for creating risk of bias graphs 30 . Tableau was used for additional plotting.

Results
2317 publications were initially identified from Pubmed. Following exclusion of publications using the criteria detailed in the Methods, 163 papers, each containing a single trial with at least two treatment arms, remained eligible for analysis ( Fig. 1). Supplementary Appendix 1 lists the references of these trials. Some key characteristics of the eligible trials appear in Table 1.

Regulatory status of treatment interventions.
Of the 126 trials that evaluated small molecules (drugs) and biologics tested singly or in combinations, 52 tested drugs or biologics that had been licensed in humans as of December 2018. In all 52 trials, the licensing occurred at least five years prior to the publication of the trial conducted in dogs. The median time between human licensing and publication of the veterinary trial was 26.5 years (IQR 18-38.25). Of the 74 trials assessing drugs, biologics, or both, that had not been previously approved by the FDA for use in humans; 59 solely tested non-FDA licensed interventions, while 15 investigated a mix of licensed and unlicensed interventions. 7 of the 74 interventions were licensed in humans after trials in dogs. In 34 of the remaining 67, the unlicensed drugs or biologics had previously been used in human studies that were published before the publication of the respective veterinary trial. 37 trials tested interventions that had never been studied in humans and 3 trials investigated a mix of previously tested and untested interventions.
Sample sizes and power. The sample size distributions were heavily right skewed with the majority of studies having low sample sizes that severely limit their detectable hazard ratios (Fig. 2a-c). The medians for total study sample size, sample size by tumor type, and by intervention arm were 33(IQR 19.5-54), 11.5 (IQR 2-32), and 11 (IQR 5-24.5) respectively. The detectable hazard ratios (HRs) given the number of events reported in these clinical trials in dogs were extremely low as well (Fig. 2d,e). The median detectable HR was 0.3 (IQR 0.18-0.40) for overall survival, and 0.06 (IQR 0.004-0.24) for disease progression. Only 2 trials would have 80% power to detect a HR of 0.7 for overall survival and only 1 trial would have 80% power to detect a HR of 0.7 for disease progression. Among the 59 trials solely evaluating drugs or biologics not previously licensed in humans, medians for total study sample size, sample size by tumor type, and by intervention arm were 42 (IQR 23-70.5), 11 (IQR 2-35.5), and 12 (IQR 4-27) respectively.

Changes in response criteria over time.
In 2010 and 2013, the Veterinary Cooperative Oncology Group (VCOG) published consensus documents describing the canine response evaluation criteria for peripheral nodal lymphoma and for solid tumors (cRECIST v1.1) respectively 31,32 . These guidelines were based on the human RECIST guidelines, and introduced a standardized framework which require lesions to be measured with standardized cutoffs defined for progressive disease, partial response, complete response, and stable disease 33,34 . cRE-CIST defines partial response as >30% shrinkage in mean sum longest diameter for all target lesions compared with baseline; progressive disease is defined as >20% growth in mean sum longest diameter for all target lesions or appearance of new lesions 32 . www.nature.com/scientificreports www.nature.com/scientificreports/ Of the 112 trials that measured disease progression as an outcome, 94 clearly defined an objective measurement criterion. Of these, only 30 explicitly named cRECIST guidelines as the criterion used. The majority of the remaining trials (57/64) were published before the VCOG cRECIST guidelines became available. The 64 trials that did not cite cRECIST as the criterion measured progression as either the appearance of new lesions (27/64) or through changes in tumor size (37/64). Definitions of tumor size varied. The majority of these trials (29/38) defined tumor size as volume or product of diameters. Only 1 trial defined size as diameter. 7 of the 38 trials did not clearly define the metric used to assess tumor size.
Some trials that did not explicitly cite cRECIST guidelines differ from the cRECIST criteria in terms of stringency. All 27 trials that defined disease progression as appearance of new lesions, and the singular trial that defined progression according to changes in tumor diameter implemented thresholds that were equivalent to cRECIST. For the 29 trials that defined tumor size by volume, a 50% reduction in tumor volume is equivalent to a 30% reduction in diameter for partial response. Likewise, a 44% growth in volume is equivalent to a 20% growth in diameter for progressive disease. For partial response, the majority of these trials (28/29) used a 50% reduction in volume and hence have equivalent stringency to cRECIST. However, significant variation existed in the percent growth definition for progressive disease. 14 of the 29 trials used thresholds at 50%, representing a threshold that is less stringent. 7 trials used thresholds at 25%, representing more stringent thresholds. In addition, a further 7 used a more stringent criteria by considering any tumor shrinkage less than partial response as no response and did not make a distinction between no response and progressive disease.
Risk of bias assessment. We identified a high degree of unreported experimental details in the included veterinary clinical trials. This underreporting was prevalent in all categories of the Cochrane risk of bias tool, with the Selective Reporting, Blinded Assessment, and Blinded Participants/Personnel categories possessing the highest share at 94% (153/163), 77% (125/163), and 60% (98/163) not reported respectively (Fig. 3). None of the trials reported pre-registered protocols.

Discussion
Overall, our systematic review of 163 clinical trials enrolling tumor-bearing dogs covered a broad range of interventions. The large majority of them pertained to the evaluation of small molecules or biologics and most were interventions that had already been tested or approved for use in humans. However, a small number of trials addressed interventions without human precedent of use. Our assessment revealed three major insights about our population of eligible trials: low sample sizes limiting statistical power; increased leniency of response criteria over time; significant underreporting of outcomes and other important study design flaws.
Overall, we observed extremely low sample sizes in most studies, in line with a prior assessment of power and sample size in canine and feline trials 13 . This was true for total study size, size by intervention arm, as well as size broken down by tumor type. We observed that it was common for trials to evaluate the effect of an intervention on multiple tumor types within the same study. However, the low sample size per tumor type per study suggests that these trials may not have had sufficient power to detect differences in effect across multiple tumor types. The typical trial would be able to detect only unrealistically super-effective treatments both for survival and disease progression and few trials had sufficient power to detect very effective treatments. Indeed, based on the median HRs, current clinical trials in dogs were only sufficiently powered to detect a clinical benefit in 1 of the 18 oncology drugs approved by the FDA between 2000 and 2011 35 , if effect sizes in dogs were similar to those seen in humans. Current trials mostly give some preliminary data on how these drugs perform in dogs, but the estimated effect sizes carry very large uncertainty.
This may be acceptable for repurposing human drugs for use in dogs, since they typically already have extensive evidence from their previous use in humans; moreover, having some clinical trial evidence in dogs may be better than having no veterinary evidence at all, at least from a feasibility and safety standpoint. However, the above argument assumes considerable overlap in disease biology and physiology between the two species, which may not necessarily hold true in all cases. Furthermore, if the primary goal is to augment the drug development pipeline (i.e. test new drugs and biologics in veterinary clinical trials before testing them in humans), such very small sample sizes are not very informative. The sample sizes for interventions not previously licensed in humans were not observed to be appreciably larger.
The field would benefit from the conduct of larger trials using consortia and multi-center mechanisms, as has been done in human trials. This approach is also starting to be adopted for murine trials. The recently established Comparative Oncology Trials Consortium (COTC) is an encouraging step in the right direction 2,4,9 . However, all completed trials conducted by the COTC had a target enrollment of between 15-30, with the exception of one trial with a target enrollment of 80-90 9 . This indicates that there may be issues beyond awareness and collaboration that continue to limit sample size even in a multi-center setting, for instance, limited resources or difficulty in recruiting willing owners and eligible participants.
Of note, one might expect that the number of events and hence detectable HR for disease progression would be higher than that of overall survival. However, the results show the opposite. We believe that the observed data arose because there were many more studies reporting disease progression outcomes alone. We observed that a large proportion of the studies that listed only disease progression outcomes were smaller than the ones that listed both progression and survival outcomes. It is these small studies reporting only progression outcomes that led to a lower detectable HR for progression than survival. www.nature.com/scientificreports www.nature.com/scientificreports/ We observed that studies which did not adopt the VCOG cRECIST v1.0 guidelines adopted a variety of different thresholds for progressive disease. These thresholds were evenly split between being more and less stringent than VCOG cRECIST. In addition, we identified studies that failed to define the criterion or size metric used to measure disease progression.
In light of these observations, we largely applaud the widespread adoption of standards such as RECIST. Standardization makes comparative assessments across trials more straightforward. Beyond establishing a common set of thresholds, RECIST also defines a standardized methodology for measurement of lesions and tumors. For example, it includes clear guidance on methods that can and cannot be used for lesion measurement 31 . This enables the comparison of multiple studies and interventions, which was previously not possible due to heterogeneity of methods. However, it is crucial to recognize that these response thresholds are arbitrarily set with no inherent biological significance 34 . The level of tolerance for low efficacy may have to shift depending on a multitude of factors. This could include the efficacy of alternative therapies available, as well as the known adverse events associated with the investigational therapy. Candidate drugs with a lower efficacy are more likely to be deemed worthy of further investigation if no effective alternatives are available, or if they present a lower incidence of adverse events compared with the standard of care. Therefore, having continuous, detailed data on tumor size changes alongside RECIST-defined progression and response data would likely improve the ability to assess treatment efficacy with granularity, facilitating comparison with other drugs and potential interventions.
We observed a high degree of unreported experimental details as well as undefined outcomes in the clinical trials in dogs. In comparison, an analysis of 20920 human clinical trials conducted from 1986 to 2014 showed a similar proportion of unreported details for human trials for random sequence generation and allocation concealment; and a lower proportion of unreported details for the remaining categories 36 . The most recent human trials published between 2011 and 2014 showed a lower proportion of underreporting for all categories 36 . Clinical trials in dogs may lag behind human trials in this aspect, however, comparisons need to be tempered by the fact that even for human trials there is large variability in the risk of bias items across different disciplines.
Lack of methodological transparency is a considerable impediment for veterinarians hoping to draw trustworthy conclusions, and for investigators looking to prioritize interventions from veterinary trials into human studies. Moreover, the lack of reported details could be symptomatic of a more pernicious problem. Instead of being due to a simple lack of reporting, this might point to an outright failure to consider these sources of bias, or to define primary outcomes during planning. If these common sources of bias are frequently overlooked in this field, the internal validity of such veterinary clinical trial results might be at high risk. The relationship between lack of reporting and biased estimates of efficacy has been well explored. Studies that do not report experimental details, such as blinding, randomization, and allocation concealment, are consistently observed to yield larger estimates of efficacy. This phenomenon was observed in rodent models 37 , veterinary trials 17 , and human trials 38 . Likewise, the observation that primary outcomes are seldom defined begets questions about the prevalence of post hoc manipulation of outcomes and p-hacking practices 39,40 in companion animal trials. Beyond the implications for veterinary medicine, such suspicions could stand to prevent such trials from being seriously considered as a trustworthy link in the drug development pathway.
Not reporting methodological details does not necessarily mean that a trial did not use proper methodological safeguards. A study using self-volunteered information from authors found that human trials which failed to report some experimental details actually utilized these approaches in execution 41 . In one evaluation, 96% (52/54) of studies that failed to report concealment of allocation undertook allocation concealment measures in reality 41 . This proportion was 79% (64/81) for blinding of outcome assessors but only 20% for blinding of participants 41 (5/25). Nonetheless, this finding should not be presumed to hold equally true for studies in dogs. CONSORT reporting guidelines for human randomized controlled trials were first proposed in 1993 42 , the same year in which the Cochrane Collaboration was founded 43 . The field has had two decades to critique issues, design standards, and change mindsets. Though some effort has been made in recent years 13,14,16,19,44 , veterinary trials have not historically received the same scrutiny and emphasis as human trials, and as such there is no guarantee that studies which failed to report details have indeed implemented sufficient bias safeguards. Therefore, even with this qualification, the very presence of high degree of unreported details contributes to uncertainty about the quality of canine trials and should be addressed. Journals and editors could take the lead by requesting adoption of more informative reporting on major design aspects that can improve the quality of trials with translational benefit.
In human trials, reporting standards have increased the proportions of important experimental details in publications [45][46][47] . Other community-initiated policies, such as the mandatory registration of trial design initiated by the International Committee of Medical Journal Editors (ICMJE) 48 resulted in rapid increases in registration and improved compliance 49,50 . The aforementioned analysis of 20920 human clinical trials showed a steady reduction in the number of studies at unclear risk over time, across all categories in the Cochrane risk of bias assessment 36 .
Adoption of registration for veterinary clinical trials may also have a similar positive impact in this regard 12 . The American Veterinary Medical Association Animal Health Studies Database (AAHSD), which was launched in 2016, includes 345 entries of trials in dogs as of December 2018. However, unlike the ICMJE policy in human trials, registration is not required for veterinary clinical trials 12 . Furthermore, only 1 registered canine oncology trial had results available with a link to the publication. In addition, the AAHSD is only limited to studies conducted in the United States, Canada, and the United Kingdom 51 , further limiting the adoption of registration in veterinary medicine. Improved transparency may also allow the field to map itself better, enhance visibility of trials and provide stimuli for collaboration towards larger studies in the field.
Additionally, raw, individual patient level data was not always made available in the analyzed studies. None of the five biggest trials reported raw data [52][53][54][55][56] . Smaller trials were more varied, with some reporting individual patient level data [57][58][59] , while others did not 60,61 . The presence of such raw data might allow pooling of data across studies for analysis, which could be especially relevant in light of the low sample sizes observed.
www.nature.com/scientificreports www.nature.com/scientificreports/ Experience from human trials shows that often, optional data elements may be missing from registrations, most commonly involving primary and secondary outcomes as well as trial end dates 49,62 . In addition, differences in protocol details can be observed between the registry and published data in a substantial portion of trials 62 . Again, the existing experience from human trials may help to inform efforts in trials in dogs, especially regarding the effectiveness of solutions and ways in which these issues can be pre-emptively adverted.
Some limitations should be discussed. Clinical trials in domestic dogs are not typically reported in phases, unlike human trials 2 . As such, we were unable to distinguish exploratory, early stage trials from late stage trials in our HR analysis. The bulk of small trials with low detectable HRs could well be exploratory trials where statistical power is less critical than late stage trials meant to assess efficacy. Likewise, in the considerable proportion of trials evaluating licensed human therapies for use in dogs, statistical power may be less critical than in trials aiming to translate untested therapies for use in humans. This is because the existing body of efficacy data in humans for repurposed drugs could help address the uncertainty of estimates derived from small veterinary trials. Nonetheless, even with this caveat, the results still hold that the vast majority of trials in dogs have very limited statistical power.
In examining the population of companion dog trials evaluating drugs or biologics not previously approved for humans, we are unable to distinguish trials conducted with the aim of translating therapies for human use from those conducted for the licensing of veterinary-specific treatments. Nonetheless, both types of trials are being conducted by the same community of investigators regardless of primary objective, thus we believe that these insights remain broadly applicable.
Finally, we chose to focus on trials conducted on companion dogs with cancer due to the rapidly increasing prominence of this field. Indeed, all completed and planned trials by the COTC recruited dogs. However, more work should be conducted to determine if our findings are generalizable to other veterinary fields that focus on other animals or other disease types.
It would also be interesting to assess in future work the extent to which results of clinical studies in dogs are predictive of results in human trials. This would require ideally the evaluation of trials that were intended to precede human experimentation versus subsequent human trials. However, even with trials where the companion animal study has typically followed human use, it may be possible to examine whether animal and human trials of the same intervention and for the same cancer yield similar estimates of effect sizes. This will likely require the availability of more, larger trials in dogs with less uncertainty about effect size estimates in order to get more definitive evidence. Beyond prediction of efficacy in humans, other proposed applications of veterinary trials, such as the establishment of biomarkers or testing of drug delivery, should also be assessed empirically in the future, as a larger body of literature on these fronts becomes available. These will be crucial in helping to maximize the translational potential of such trials as part of the drug discovery pipeline.
In closing, our comprehensive assessment of clinical trials in tumor bearing companion dogs reveals numerous areas for improvement. These areas comprise low sample sizes, unclear reporting of response criteria and possible increased leniency over time, a high degree of underreporting of key protocol details outlined in the Cochrane risk of bias tool, and prevalent failure to define primary outcomes. This study builds upon the existing body of evidence examining similar pitfalls in veterinary studies 13,14,[16][17][18][19] , underscoring the importance of addressing these shortcomings for this research community. In particular, with trials on dogs starting to play a role in human drug discovery, this is an especially opportune time for increased awareness and improvements to the conduct of research.

Data Availability
Data and analysis scripts used in this study have been deposited on the Open Science Framework and are accessible at https://osf.io/5e37z/?view_only=e70e94fbe6824c9cac75c3e1e03a1c5b.