s2fieldsofstudy: Biology Chemistry

text:
Different duplex / quadruplex junctions determine the properties of antithrombin aptamers with mixed folding

Mixed duplex/quadruplex oligonucleotides have attracted great interest as therapeutic targets as well as effective biomedical aptamers. In the case of thrombin-binding aptamer (TBA), the addition of a duplex motif to the G-quadruplex module improves the aptamer resistance to biodegradation and the affinity for thrombin. In particular, the mixed oligonucleotide RE31 is significantly more effective than TBA in anticoagulation experiments and shows a slower disappearance rate in human plasma and blood. In the crystal structure of the complex with thrombin, RE31 adopts an elongated structure in which the duplex and quadruplex regions are perfectly stacked on top of each other, firmly connected by a wellstructured junction. The lock-and-key shape complementarity between the TT loops of the G-quadruplex and the protein exosite I gives rise to the basic interaction that stabilizes the complex. However, our data suggest that the duplex motif may have an active role in determining the greater anti-thrombin activity in biological fluids with respect to TBA. This work gives new information on mixed oligonucleotides and highlights the importance of structural data on duplex/quadruplex junctions, which appear to be varied, unpredictable, and fundamental in determining the aptamer functional properties.

INTRODUCTION
DNA structures, different from the canonical B-form double helix, have been known for a long time and extensively characterized. Several studies have also highlighted the remarkable shape variability of nucleic acids (1) and its close connection to the specific functions of these molecules. The alternative folds include various kind of double helices (2), triplex arrangements (3), four way Holliday junctions (4), i-motifs (5), and G-quadruplexes (6). The latter are particularly interesting because have been found in telomeric ends of chromosomes (7), in oncogenes (8) and in the promoter regions (9).
G-quadruplexes are formed by stacked G-tetrads, which are cyclic arrays of four guanines interacting through Hoogsteen hydrogen bonds. Depending on strand directionality, glycosidic bond angles, loop length and conformation, different G-quadruplex topologies have been observed (6).
Thanks to the numerous crystallographic (10) and NMR (11) structures of G-quadruplexes determined in the last decades, our knowledge of these peculiar motifs is considerably grown. Thus, in the last years the interest has moved towards the way these structures interface with other structural motifs (12), such as other quadruplexes (13)(14)(15) or duplexes (12,16,17). In particular, an intriguing issue concerns the way duplexes and quadruplexes, which are often characterized by different strand directionality, groove width and base orientation, can found accommodation in the same molecule. Recently, Phan et al. (17,18) have shown, through the NMR characterization of artificial constructs, several modes in which quadruplexes and duplexes can be connected (17). However, the molecules they investigated were designed to adopt a specific fold through the insertion of strand breakers or additional residues in key regions to promote the quadruplex-duplex transition. Much less is known about the fold of oligonucleotides in which these specific markers are not present. This is of particular interest in the case of aptamers, DNA (19) and RNA (20) strands able to bind with very high affinity and selectivity a large variety of molecules, ranging from small molecules, other oligonucleotides, peptides or large proteins (21). This ability mainly depends on the aptamer propensity to fold into several distinct three-dimensional shapes, which are difficult to establish a priori. Indeed, aptamers against a specific target are usually selected from random libraries of about 10 10 -10 20 sequences by SELEX (21). Using this procedure, aptamers against many pathologically relevant proteins (22), such as thrombin, nucleolin, HIV RNase H, prion and many others, have been selected (23)(24)(25)(26). The variety of folds these oligonucleotides may adopt includes also a duplex/quadruplex mixed structure. Based on their sequence, this very peculiar structure has been suggested to be adopted by an anti-RANK (human receptor activator of NF-kB) aptamer, the aptamer against the HIV RNaseH and several anti-thrombin aptamers (23), such as HD22 and RE31. HD22 aptamers (24) include a 27-and a 29-mer able to bind the exosite II of thrombin with an affinity about 100-200 times higher than that of the best known thrombin binding aptamer TBA, which recognizes exosite I (19). The 31-mer RE31 (25) is a new generation aptamer with increased affinity towards exosite I (26,27) and improved performances in biological fluids (28) with respect to TBA. The latter has been shown (29)(30)(31)(32) to adopt a unimolecular G-quadruplex structure that includes two Gtetrads, two TT loops and a TGT loop ( Figure 1A). Moreover, X-ray diffraction analysis has clearly shown that the binding motif of this molecule involves the two TT loops (29)(30)(31). HD22 and RE31 were supposed to have a TBA-like G-quadruplex ( Figure 1B and C) and a terminal duplex segment connected by two linkers three-and two-residue long for HD22 (24) and RE31 (26)(27)(28)32), respectively. The proposed global fold for the two aptamers is rather similar in contrast to their very different properties.
Some of us have recently solved the crystal structure of HD22 in complex with thrombin (33) demonstrating that it adopts a fold different from that previously proposed and depicted in Figure 1B. The structure embodies a duplex segment strictly connected to a highly distorted G-quadruplex (pseudo-quadruplex), in an unexpected perpendicular arrangement of the two regions. This peculiar steric combination favours the binding at the exosite II, with which both motifs interact (33).
Here, we present the crystal structure of the complex between thrombin and RE31, which clearly shows that the aptamer binds at the exosite I and adopts a markedly different fold with respect to HD22. The aptamer presents an elongated structure in which three base pairs mediate the transit from the duplex to the quadruplex region, determining a continuous base stacking that includes both regions. The results also provide a structural based interpretation of the improved binding and inhibitory properties of RE31 compared to TBA (26).
A stock solution of the aptamer was prepared at a concentration of about 1 mM in 100 mM potassium chloride, 10 mM potassium cacodylate pH 7.0. This solution was heated for 10 min at 85 • C and then slowly cooled down and stored at 20 • C overnight, in order to induce folding. The correct fold was verified by recording circular dichroic spectra in the 200-350 nm range. The complex with thrombin was prepared as previously described (29,30,34) by placing a 1.2fold molar excess of aptamer on a frozen sample of inhibited thrombin and leaving the sample for 3 h at 4 • C. Then, the sample was diluted and the buffer changed to 25 mM Kphosphate pH 7.1 and 100 mM KCl. The thrombin-RE31 complex was finally concentrated to about 7 mg/ml using 10 kDa-cutoff Centricon mini-concentrator and a refrigerated centrifuge.
Initial crystallization trials were set up using commercial kits and an automated crystallization workstation (Hamilton Robotics). In a couple of days, well-formed crystals grew in many different conditions. The best diffracting ones (Supplementary Figure S1) grew from drops consisting of 0.6 l thrombin-RE31 solution and 0.6 l reservoir solution (200 mM sodium malonate pH 7, 18% w/v PEG 3350 Da) that were equilibrated against 500 l reservoir solution.

Data collection, structure determination and refinement
Diffraction data were collected in-house on a Saturn944 CCD detector. The X-ray radiation used was Cu-K␣ radiation from a Rigaku Micromax 007 HF generator. After the addition of 25% (v/v) glycerol to the harvesting solution, crystals were flash-cooled at 100 K in supercooled N 2 gas (Oxford Cryosystems) and maintained at 100 K during the data collection. Diffraction data were processed using HKL2000 (35). Crystals of thrombin-RE31 complex belong to the orthorhombic space group P2 1 2 1 2 and diffract up to 2.99Å resolution. Matthews' coefficient calculations suggested the presence of a 1:1 complex in the asymmetric unit and a solvent content of 45%. Detailed statistics on data collection are reported in Supplementary Table S1.
The structure was solved using the program Phaser (36). The coordinates of thrombin and the two G-planes of TBA, which were derived from the structure of the thrombin-TBA complex (PDB code: 4DII) (30), were utilized as two independent search models. A clear solution was obtained with a Z-score of 14.3 and a log-likelihood gain (LLG) of 848.
The starting model was refined using Refmac5 (37) and Phenix (38) programs. Each run was alternated with manual model building using WinCoot (39). Analysis of electron density maps, calculated with (F o -F c ) and (2F o -F c ) coefficients, allowed the rebuilding of the aptamer missing residues, with the only exception of the two terminal 5 and 3 residues, the positioning of a potassium ion in the quadruplex core and of several water molecules.
The final model has R-factor/Rfree values of 0.166/0.231. Statistics and parameters of the refinement are given in Supplementary Table S2.
The drawings were prepared with Pymol (http://pymol. org) and Chimera (40). The coordinates of the structure have been deposited in the Protein Data Bank (Code 5CMX).

Structural analysis
3DNA-dssr (41) was used to calculate local and overall geometric parameters of the aptamer. Superpose program from CCP4 package (42) was used to calculate root mean square deviations. Features of the thrombin-RE31 interface were calculated using Cocomaps server (43), whereas contacts between the two molecules, as well as packing interactions between the aptamer and symmetry related thrombin molecules, were found by using 3DNA-snap (41) and Pisa (44) programs. All the results were verified by visual inspection of the structure with WinCoot (39).

Overall crystal structure
Best crystals of the complex between thrombin and RE31 in the presence of K + , diffract X-rays up to a resolution of 2.99Å and belong to orthorhombic space group P2 1 2 1 2. Initial model obtained by molecular replacement showed that the aptamer G-planes are located close to thrombin exosite I. Electron density maps calculated from initial phases allowed unambiguous model building of the whole DNA molecule (Figure 2A and B), with the only exception of the terminal residues (Gua1 and Cyt31). For what concerns thrombin, the heavy chain (residues 16-245) and the light chain (residues 1B-14K) were well defined in the electron density maps, with the exception of the ␥ -autolysis loop (residues 146-150) of the heavy chain. The crystal structure was refined to Rfactor/Rfree values of 0.166/0.231, respectively. Detailed statistics of the refinement process are reported in Supplementary Table S2.

Global aptamer fold
Based on sequence analysis and spectrophotometric studies in solution, a mixed duplex-quadruplex structure was proposed for RE31 aptamer ( Figure 1C) (25,28,32); crystallographic analysis shows that RE31 in complex with thrombin indeed adopts this mixed fold.
Moreover, good quality electron density maps ( Figure  2B) allowed a very detailed description of the aptamer fold. Five Watson-Crick base pairs, starting with the Thy2-Ade30 pair, define a regular B-type double helix. Two additional non-Watson-Crick base pairs formed by residues belonging to linkers and loops extend further the duplex region. The quadruplex motif involves the fifteen TBAresidues arranged in the well-known chair-like fold (30), with two G-tetrads (residues 10, 13, 19, 22 for tetrad I, and 9, 14, 18, 23 for tetrad II), which coordinate a potassium ion, two TT loops and a TGT loop. The latter is in close proximity to the duplex segment, and, in particular, is directly involved in the transition from the duplex to the quadruplex segment (see the following section). After superposition of the quadruplex core of RE31 and TBA, excluding the TGT loop, the root mean square deviation (RMSD) between all the atoms of the two G-quadruplexes is only 0.89Å, indicating that the overall fold of the quadruplex remains basically unaltered. The quadruplex domain of RE31 keeps the alternating syn-anti conformation for each adjacent pair of guanines along the strand and the cyclic array of Hoogsteen-like hydrogen bonds in each tetrad.
The two tetrads are essentially planar, with RMSD values calculated with respect to the average G-quartet planes of 0.37 and 0.42Å for tetrads I and II, respectively, in very good agreement with values obtained for TBA (0.35 and 0.39Å, for TBA in complex with thrombin in the presence of K + ions, PDB code 4DII). The similarity between TBA and the TBA-portion of RE31 extends to position of bases belonging to TT loops. Indeed, thymines 12 and 21 stack on the G-tetrad I and are connected to each other through a water molecule, similarly to the Thy4-Thy13 pair of TBA. Minimal differences are observed for Gua18, whose position is fixed by the stacking geometry of Thy17 and the other core guanines (see below).
Careful analysis of the duplex region (Supplementary Table S3) highlights the regularity of the structure, despite the presence of the quadruplex, likely due to the gradual transition from one motif to the other.

Duplex/quadruplex junction
On the basis of the sequence, the duplex and the quadruplex motifs of RE31 were supposed (25,28,32) to be connected by two 2-residue linkers (TA and GG respectively) ( Figure  1C). This suggestion is not validated by the crystal structure. Indeed the TGT loop from the quadruplex region is involved in the transition between the two domains by forming base pairs with residues from the linkers. Overall, three base pairs constitute the duplex/quadruplex junction (Fig-ure 3). At the end of the double helix, the first two residues from the linkers (Thy7 and Gua25) form a Wobble pair. In the second pair (Ade8 and Thy17), the bases are linked through a reverse Hoogsteen hydrogen bond network. It is worth to stress that Thy17 is part of the TGT loop. The last loose pair presents a non-canonical G-G interaction between Gua24 and Gua16; the latter is the guanine of the TGT loop. Detailed hydrogen-bonded interactions among residues of the junction are reported in Table 1 and depicted in Supplementary Figure S2. It should be stressed that the first base pair, Thy7-Gua25 is planar and well integrated in the duplex segment; whereas the last one, Gua16-Gua24, barely defines a plane. However, the arrangement of these six residues is such that a continuous stacking from the duplex segment to the quadruplex one is observed ( Table 2). In particular, Thy7, Thy17 and Gua16 (junction) are aligned with Gua6 (duplex) and Gua14 (quadruplex) and Gua25, Ade8 and Gua24 (junction) are aligned with Cyt26 (duplex) and Gua23 (quadruplex) (Figure 4). The junction pairs are also crucial to mediate the change of groove width from 17 A (measured for the P-P distance) for the minor groove of the duplex to 19Å for the wide groove of the quadruplex.

Thrombin-RE31 interaction
In the asymmetric unit RE31 binds exosite I of thrombin through its quadruplex motif in a manner very similar to TBA (30). Detailed interactions between thrombin and RE31 are reported in Table 3. Comparing these data with those previously obtained for TBA (30), it is clear that  the differences in the mode of binding of the two molecules are negligible, despite the doubled size of RE31. Hydrogen bonds, as well as hydrophobic interactions, with the target protein are very well conserved ( Figure 5). It should be noted that none of the residues belonging to the duplex or the junction region of RE31 is involved in the stabilization of the binary complex.

Crystal packing
In addition to the aptamer-thrombin interactions described above, the duplex segment of RE31 is in contact with three symmetry related thrombin molecules (Supplemen-tary Figure S3). In particular, the phosphate group of Thy28 strongly interacts with two arginine side chains of the molecule at −x-1/2, y-1/2, −z-1: salt bridges are formed between OP1 and OP2 of Thy28 and guanidinium groups of Arg126 and Arg233, belonging to thrombin exosite II. Less strong interactions are found between RE31 and the thrombin molecule at −x-1/2, y-1/2, −z-2. In this case the protein caps the end of the duplex with residues Pro60B, Pro60C, Trp96 and Arg97 burying an interface area of 158Å 2 . Finally, Thy15, one of the three residues not involved in intramolecular contacts (the other two, Thy11 and Thy20, bind exosite I), occupies a patch on surface of the molecule at −x-1, −y, z, defined by residues Phe34,    Lys36, Gln38, Leu65, Arg67, Tyr76, Lys81, Ile82, Met84. Thy15 also contacts Thy20 of the corresponding aptamer.
Although the interface areas of the duplex contacts are appreciable, they are much less extended than that relative to the well-characterized TT loops-exosite I interface (586.2 A 2 ).

DISCUSSION
The structural investigation on oligonucleotides able to adopt a mixed duplex/quadruplex fold is still understated, despite the great interest these molecules assume in different biological contexts, including the case of several biomedical aptamers (23).
The present crystallographic analysis shows that RE31 adopts a mixed duplex/quadruplex fold, as it was previously suggested using a combination of spectrophotometric studies in solution and sequence analogies with well characterized molecules such as the 15-mer TBA with which RE31 share the sequence of the core residues. Indeed, these fifteen residues adopt the antiparallel G-quadruplex structure observed for TBA (30). However, in RE31 the TGT loop is crucially involved in the transit from the duplex to the quadruplex region with the result that a continuous base stacking runs along the whole molecule: the duplex gradually converts into the quadruplex organization and the helical axes of the two regions maintain almost the same orientation. This rod-like shape prevents a simultaneous interac-tion of both the duplex and the quadruplex region with the same thrombin molecule, and the binding does not involve the duplex, but only the quadruplex domain, as observed for TBA-like aptamers (29)(30)(31). This behaviour is markedly different from that of HD22, the other duplex/quadruplex thrombin-binding aptamer, whose structure has been elucidated by X-ray crystallography (33). In HD22 aptamer, the first guanine of the quadruplex is directly linked to the last residue of the duplex and the transition is mediated by a G-fork between the first guanine of the quadruplex and a linker residue. This kind of junction determines the formation of a distorted G-quadruplex structure and an almost orthogonal orientation of the helical axes of the two motifs. In this bend conformation both the duplex and the quadruplex can interact with exosite II of thrombin, which has a significantly larger surface than exosite I, the binding site of RE31 and TBA-like nucleotides.
In the present crystal, the duplex of RE31 establishes contacts with three different regions of the protein surface of symmetry related thrombin molecules. Although these contacts involve a smaller number of interactions with respect to the aptamer/protein binding site described above, they may contribute to the greater protein affinity measured for RE31 with respect to TBA (27). On the other hand, the compact structure of the aptamer, and in particular the recruitment of the TGT loop in the junction, may increase the resistance to nucleases that primarily attack single-stranded oligonucleotide fragments (quadruplex loops) (45,46) and explain the greater inhibitory effect of RE31 compared to that provided by TBA. Furthermore, it can be surmised that the presence of the duplex phosphodiester bonds, whose rupture is not likely to affect the functionality of the aptamer, may also contribute to a decrease of the nuclease efficiency.
As a final comment, it is worth to stress that RE31 and HD22, which have the duplex motif at the 5 and 3 ends of a quadruplex core, exemplify a different structural case with respect to the constructs designed and characterized by Phan and co-workers (17,18,47,48). In the latter case, the duplex region is inserted in a loop of well-structured G-quadruplexes, where it produces negligible effects on the quadruplex architecture. The RE31 model shows that the link of an antiparallel B-form duplex to an antiparallel Gquadruplex may occur with an almost gradual transition from one structural motif to the other, allowing the regularity of both domains to be preserved through the formation of a continuous base stacking. On the other hand, for HD22 the preservation of the duplex structural features, severely affects the regularity of the quadruplex (33).
In conclusion, the present thrombin-RE31 structure represents the first atomic resolution model of a hybrid duplexquadruplex aptamer interacting with protein exosite I. This structure, in comparison with the different architecture of HD22, proves that structural characterization of duplex/quadruplex and quadruplex/duplex molecules is a fascinating topic worth to be deeply analysed. Indeed the sterically accessible space for residues that connect the two motifs allows varied and unpredictable conformations of the link and appears to play an important role in determining functional properties of these oligonucleotides.

SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.


Naked-Eye Fingerprinting of Single Nucleotide Polymorphisms on Psoriasis Patients

We report a low-cost test, based on gold nanoparticles, for the colorimetric (naked-eye) fingerprinting of a panel of single nucleotide polymorphisms (SNPs), relevant for the personalized therapy of psoriasis. Such pharmacogenomic tests are not routinely performed on psoriasis patients, due to the high cost of standard technologies. We demonstrated high sensitivity and specificity of our colorimetric test by validating it on a cohort of 30 patients, through double-blind comparison with two state-of-the-art instrumental techniques, namely reverse dot blotting and sequencing, finding 100% agreement. This test offers high parallelization capabilities and can be easily generalized to other SNPs of clinical relevance, finding broad utility in diagnostics and pharmacogenomics. Page 2 of 23 Nanoscale N an os ca le A cc ep te d M an us cr ip t Pu bl is he d on 2 8 A pr il 20 16 . D ow nl oa de d by N an ya ng T ec hn ol og ic al U ni ve rs ity o n 07 /0 5/ 20 16 0 9: 46 :0 2. View Article Online DOI: 10.1039/C6NR02200F

Sequencing of the human genome has shown that a significant part of the genetic variability among individuals is due to single nucleotide polymorphisms (SNPs). 1 Nowadays, a bulk of knowledge is being accumulated on the involvement of an increasing number of SNPs in constitutional susceptibility to diseases and predictability of the response to drugs. 2,3 These relationships are the focus of pharmacogenomics, the genetic pillar at the basis of personalized medicine. Current technologies for SNP diagnostics are mainly based on sequencing, real-time PCR, and hybridization assays, which however have limited clinical applicability, owing to their significant costs and instrument requirements, thus hampering the implementation of personalized therapies. 4 Psoriasis is a chronic inflammatory skin disease affecting 2% of the world population. 5 It is a heterogeneous disease, both in its clinical symptoms and in the underlying genetics. In recent years, Genome-Wide Association Studies (GWAS) have demonstrated the role of chromosomal regions, called PSORiasis Susceptibility (PSORS), in the risk of developing psoriasis. Moreover, several lines of evidence support the hypothesis that patients having a particular psoriasis-associated SNP can respond to treatment better than others. 6 In particular, the HLA-C*06:02 (Human Leucocyte Antigen-C*06:02), located in the PSORS1 region on chromosome 6, was the first identified allele associated with psoriasis susceptibility 7,8 and correlates with disease severity. 9 HLA-C*06:02 positive patients showed a better response to the biological drug ustekinumab. 10 Moreover, several SNPs in many other genes, mainly involved in the immune response, are associated with psoriasis. [11][12][13] Notably, patients harboring certain mutations in the TNFAIP3 (Tumor Necrosis Factor-Alpha-Induced Protein 3) gene, respond better to anti-TNFα biological therapy. 14,15 In particular, rs610604 G allele in the TNFAIP3 gene 16 has been shown to be associated with a better response to several anti-TNFα drugs (etanercept, infliximab, and adalimumab). 15 Currently, HLA genotype and TNFAIP3 SNPs are not routinely analyzed in psoriasis patients, due to economic constraints and relative novelty of the understanding of their pharmacogenomics implications. However, recent reports demonstrated that a four SNPs tagging is sufficient to identify both the presence and zygosity of the psoriasis susceptibility allele HLA-C*06:02, thus possibly substituting the more demanding whole gene analysis. 17,18

Results
In this study, we show the low-cost detection of SNPs in the HLA locus and in the TNFAIP3 gene, aimed at identifying patients carrying the HLA-C*06:02 allele and the rs610604 G allele of TNFAIP3. We exploited a rapid and instrument-free colorimetric method for SNP detection, recently developed in combination with isothermal DNA amplification. 19 We achieved the colorimetric genotypization of 30 real patients' samples, and validated the results through double blind comparison with two standard techniques, namely reverse dot blotting and sequencing, demonstrating 100% accuracy.
In order to identify the C*06:02 allele, we exploited the 4 tagging SNPs typing (rs3130457, rs6904246, rs7745906 and rs10484554) 12,17 recently reported by Nikamo et al. 18 The fifth SNP analyzed in this test was the above mentioned rs610604, falling in the TNFAIP3 gene. A total of 30 psoriasis patients were included in the study. 28 patients had been diagnosed with psoriasis vulgaris and 2 patients with palmoplantar psoriasis. Genomic DNA from patients' blood samples was amplified by PCR with primers specific for the regions of the HLA locus and TNFAIP3 gene harboring the SNPs of interest (Table S1 †).
Raw PCR samples were directly used for the colorimetric detection, according to the scheme depicted in Fig. 1. Briefly, streptavidinated paramagnetic microparticles were used to capture biotinylated PCR amplicons, which include the target SNPs to be tested. The PCR amplicons were denatured and subsequently hybridized with two different probes, the first one specific for the SNP (discriminating probe) and the second one, linked to AuNPs, specific for a fixed portion of the discriminating probe, providing the colorimetric detection (detection probe). After incubation with the streptavidinated magnetic microparticles, magnetic washing was performed to denature the amplicons and obtain single stranded DNA targets for the subsequent hybridization with the probes. For each SNP, two aliquots of the amplicon were hybridized with the discriminating probes specific for either allele. The universal AuNP probes were then added to the reactions. These consist of 40 nm AuNPs, functionalized with polyT oligonucleotides, which recognize a polyA tail in the discriminating probes.
Details about the synthesis, functionalization, and characterization of AuNP probes are provided in Materials and methods. Finally, samples were subjected to magnetic washing and the colorimetric result was assessed visually. In particular, red color ("ON" result), due to hybridized AuNPs, corresponds to a perfect match between the discriminating probe and the sample, indicating the detection of the specific allele. Conversely, yellow color ("OFF" result), due to pristine magnetic microparticles alone, indicated mismatch among the sample and probe, i.e., the absence of the allele tested for (Fig. 1B). Finally, as a further control, the negative PCR sample was hybridized with an equimolar mixture of the two discriminating probes. The resulting yellow color confirmed the absence of PCR contaminations, without the need for electrophoretic analysis. For TNFAIP3, homozygous GG, homozygous TT and heterozygous GT patients could be directly discerned from the colorimetric results. For HLA locus, the presence or absence of the C*06:02 allele was assessed by interpreting the results shown in Fig. 2. Indeed, the presence of a particular HLA allele can be identified by different combinations of SNPs, or haplotypes. 17 Colorimetric test results are shown in Fig. 3. Each line of the plate corresponds to a patient, while each column corresponds to a discriminating probe, i.e., to an allele variant. The red/yellow experimental results corresponded well with the ON/OFF response of the test, and can be interpreted clearly and without bias by the naked-eye. Furthermore, the double test performed on each SNP allows distinguishing homozygous and heterozygous samples. For instance, according to the legend shown in Fig. 1B, the first patient is heterozygous for the first SNP (rs10484554), T allele homozygous for the second SNP (rs3130457), A allele homozygous for the third SNP (rs6904246), and so forth. The attribution of the C*06:02 allele could be easily inferred by visual comparison with a printed legend (Fig. 2). This latter describes all the possible combinations of the SNP status, corresponding to the absence or presence of the C*06:02 allele, or to heterozygosity for C*06:02. For each SNP, each of the two possible alleles may be present (red), absent (yellow) or either of the two (grey, indifferent). For instance, the first patient is C*06:02 negative (because it Fig. 1 Strategy for the colorimetric SNP fingerprinting. (A) PCR reactions were directly analyzed by the colorimetric assay, using two discriminating probes, each specific for an allele variant ("a" or "b"). The PCR product hybridizes with the matching probe only (in this example, the one specific for allele a), capturing the AuNP probes on the microparticle surface, giving the sample a red color (left). In the case of mismatch (allele b is not present in this example), the sample remained yellow (right). The PCR negative control was analyzed with a mixture of the two discriminating probes and the resulting yellow color confirmed the absence of any amplified product (not shown in the scheme). (B) Legend for the interpretation of the results. Red color indicates "ON" result, i.e., the presence of that particular SNP variant. Yellow color indicates "OFF". If both the samples containing the two discriminating probes turn red, the patient is heterozygous for that particular SNP. (C) Spot control of PCR specificity by polyacrylamide gel electrophoresis. matches the status of the 3 rd row in Fig. 2). Moreover, several different combinations of SNPs identify the presence of the C*06:02 allele in homozygosity (Fig. 2, last line). Indeed, any combination of SNPs rs6904246 and rs7745906 (homozygous for any of the alleles, or either one or both heterozygous) is compatible with C*06:02 homozygosity, while only SNP rs10484554 and SNP rs3130457 should necessarily be T and C, respectively. 18 Similarly, other combinations of SNPs indicate heterozygosity or the absence of C*06:02.
Notably, after standard PCR amplification, the duration of the colorimetric assay is ca. 30 min and requires no instrumentation. Moreover, all reagents, as well as the final colored samples, resulting from the presence or absence of hybridized AuNPs (Fig. 3, lower panel), are stable for months, allowing easy recording/archiving with no time-constraints.
All samples were also analyzed with two different state-ofthe-art technologies, i.e. direct sequencing and reverse dot blotting (Fig. S2 †). The primers used for direct sequencing are listed in Table S2. † All these studies were performed in doubleblind. Interestingly, a comparison of our colorimetric test with standard techniques showed 100% agreement (Table 1), both when using the predictive method based on reverse dot blotting and when using direct sequencing, demonstrating an insignificant risk of producing erroneous results (both false positives and false negatives). The detailed results of the SNP genotypization by colorimetric test, sequencing, and reverse dot blot test, for each patient, are provided in Table S3. † Fig. 2 Interpretation of the four SNPs tagging, 18 to determine C*06:02 presence and zygosity. Each column corresponds to a SNP variant, as indicated on top of the scheme. Each row corresponds to a possible combination of SNPs that can be found in a patient. Red color: "ON" result; yellow: "OFF"; grey: indifferent result (both presence and absence of the SNP variant return the same interpretation of the HLA-C*06:02 status). On the right, the interpretation of the HLA-C*06:02 status, for each SNP combination, is shown. Fig. 3 Results of the colorimetric test applied on 30 patients. Each row corresponds to a different patient, while each column to a different discriminating probe. For clarity, the negative control of the PCR reaction (always yellow), was not included in the picture (each patient's data row has been taken as a separate picture; pictures were then collated together for clarity of presentation). Bottom panel: representative SEM image of magnetic microparticles decorated with hybridized AuNPs (white dots).

Discussion
The availability of a robust, rapid, and cost-effective tool for SNP fingerprinting can enhance pre-treatment screening and promote personalized medicine. 4 Recently developed biological drugs for psoriasis have a better safety profile compared to conventional systemic drugs, but a cost up to ten times higher. 11 Given the high variability of the therapeutic response among different patients, 11 it would be crucial to predict the efficacy of a drug before treatment, as this could reduce unnecessary treatments and shrink the treatment-to-cure window, with great benefit to patients. However, treatments are most often administered on a trial-and-error basis, exposing the patient to avoidable side effects from ineffective drugs and delays in symptoms relief, which can be frustrating, especially for the severe forms of the disease, associated with serious and even potentially life-threatening comorbidities.
Gold standard methods to genotype HLA alleles rely either on sequencing, or on PCR with sequence-specific-primers (SSP-PCR), [20][21][22] or on array hybridization tests. In particular, a recently developed method, called the INNO-LiPA HLA-C kit, can detect human leukocyte antigen C alleles at the allele group level (first two digits of the HLA nomenclature, 23 e.g. HLA-C*06:xx), in a strip format, with one amplification and one detection line. This method does not directly detect the specific allele (second two digits in HLA nomenclature, e.g. HLA-C*06:02), but it can provide an estimate of it, based on statistical data, through the dedicated LiRAS software. Moreover, Olerup et al. have developed a SSP-PCR procedure, 20,21 which requires two steps (and two different kits) to identify first the serological group, then the allele within that group. The PCR SSP results carry a large amount of information, as several allelic variants can be identified with the same test, and have to be carefully interpreted according to reference tables provided by the manufacturer. On the other hand, typing of tagging SNPs 17 allows for simple interpretation of the results. The HLA-C*06:02 status can be assessed easily, based on a combination of four tagging SNPs 18 (Fig. 2). Our test offers a low-cost and instrument-free alternative to standard techniques for SNP genotyping. Several AuNP-based colorimetric methods have been reported in the literature for the detection of nucleic acid targets, and also for single nucleotide polymorphisms (reviewed, for instance, in ref. 24,25). Many of them exploit AuNPs aggregation and their consequent colorshift. [26][27][28][29] When comparing simplified, colorimetric SNP detection assays with state-of-the-art instrumental techniques, the main issues are specificity and sensitivity. 24 Although several reports have developed different strategies to deal with these issues, as a matter of fact, a great majority of commercial diagnostic tests for SNP genotyping currently rely on instrumental assays, mirroring the drawbacks that molecular strategies face when transferred to real-world samples. Compared to other hybridization-based methods, our test is very robust and simple to operate, and, importantly, displays zero background. This is due to a smart design that combines the use of AuNPs as a signal transducer and cooperative hybridization of probe oligonucleotides. 30 Indeed, beyond the sensitive colorimetric detection, AuNPs provide enhanced specificity in mismatch discrimination, through their direct effect on the melting properties of DNA. 31 Cooperative hybridization of adjacent probes, in turn, further enhances the stability of the binding of perfectly matched targets, over that of mismatched targets. In addition, the yellow-to-red colorimetric result proved to be bias-free, making the test appropriate for the development of an integrated point-of-care device, which can be operated and read even by untrained users, or by automated image analysis software (as shown, for instance, in Fig. S3 †).
In this work, we have reported and validated a novel colorimetric test for the genotypization of a panel of SNPs, relevant in the pharmacogenomics of psoriasis. Given its diagnostic accuracy, simplicity and affordability, this test may find great utility in promoting personalized medicine of psoriasis patients. Furthermore, the use of universal AuNP detection probes, together with the great parallelization potential given by the reduced number of processing steps, simplifies the application of this assay to the colorimetric fingerprinting of other clinically relevant SNP panels.

Materials and methods
Synthesis and functionalization of 40 nm gold nanoparticles 40 nm citrate capped gold nanoparticles were synthesized as previously reported. 32 AuNPs were conjugated with thiolated DNA probes consisting of a sequence of thymines (5′ T(30)-(O-CH2-CH2)3-SH 3′), purchased from IDT DNA, by standard gold-thiol chemistry. Functionalized AuNPs were centrifuged and washed four times with 0.3 M NaCl in 10 mM phosphate buffer, 0.01% SDS, in order to remove the excess unbound oligonucleotides. These conjugates DNA-AuNPs were stored at 4°C until their use. Their concentration was determined by UV-vis spectrophotometry. In order to quantify the density of functionalization, functionalized AuNPs were digested overnight at 40°C with 10 mM dithiothreitol (DTT) in order to displace the DNA from the AuNP surface. The digested suspension was centrifuged at 13 400 rpm for 15 min and the supernatant was analyzed using the Quant-iT OliGreen ssDNA kit (Invitrogen). The number of oligonucleotides bound on the surface of a single AuNP was 430 ± 60. Both the as-synthesized AuNPs and AuNP probes were thoroughly characterized by UVvis, Dynamic Light Scattering (DLS), and Transmission Electron Microscopy (TEM) (Fig. S1 †).

View Article Online
Clinical samples and extraction of genomic DNA A total of 30 psoriasis patients (21 females and 9 males, average age 47 years), under treatment at the department of Dermatology of Vito Fazzi Hospital in Lecce (Italy), were included in the study. Of these, 28 patients were affected by psoriasis vulgaris, with psoriatic lesions in specified body regions (upper and lower extremities, trunk, scalp, and nails), 2 patients by palmoplantar psoriasis. All patients provided written informed consent prior to any study-related procedures. In 4 patients the skin disorder was associated with psoriatic arthropathy (according to Moll and Wright criteria for PsA), diagnosed by imaging studies (enthesitis and tenosynovitis). PASI (Psoriasis Area and Severity Index) was used to assess psoriasis disease severity and the clinical efficacy of treatment. DNA samples for the genetic study were extracted from blood samples from all patients by using a NucliSENS® easyMAG® automated nucleic acid extraction system (BioMérieux, Durham, NC). The samples were then anonymized prior to testing. The concentration and purity of the extracted DNA were determined by measuring the absorbance at 260 nm by using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Inc.). The average concentration of genomic DNA in extracted samples was ≤5 ng µl −1 . Only 1 µl of this material was used for the amplification and yielded enough material for the subsequent hybridization step.

Live subject statement
All experiments were performed in compliance with the relevant laws and institutional guidelines. The study was approved by the ethical committee of the Vito Fazzi Hospital. Written informed consent was obtained from patients for DNA testing.

Polymerase chain reaction (PCR) amplification and gel electrophoresis
Genomic DNA was amplified by traditional PCR to obtain approximately 100 bp fragments containing SNPs rs3130457, rs6904246, and rs7745906 for the determination of HLA C*06:02 allele and SNP rs610604 in the TNFAIP3 gene. The PCR reaction was set up in a 50 μL mixture containing 100 pg μL −1 of genomic DNA, 100 nM FW primer, 100 nM biotinylated REV primer, 2.5 mM MgCl 2 , 1× PCR reaction buffer, 200 μM dNTPs (Sigma-Aldrich) and 0.05 units per μL of Taq DNA Polymerase (Sigma-Aldrich). Sequences of all primers used are listed in Table S1. † DNA amplification was performed using the following thermal cycling conditions: 95°C for 5 min, followed by 34 cycles of 94°C for 1 min, 55°C for 1 min, 72°C for 1 min, with a final extension at 72°C for 5 min. For the amplification of the region of the genome containing SNP rs10484554, which is a GC-rich region, the reaction mixture contained 200 pg μL −1 of genomic DNA, 1.875 mM MgCl 2 , 200 μM CleanAmp™ dNTP (Sigma-Aldrich) and 3% dimethyl sulfoxide (DMSO). In this latter case, thermal cycling conditions were: 95°C for 5 min, followed by 38 cycles of 94°C for 1 min, 59°C for 30 s, 72°C for 1 min, with a final extension at 72°C for 5 min. PCR products from the first sample were checked by polyacrylamide gel electrophoresis, in order to assess the specificity of the primer sets and exclude off-target amplifications. Notably, the step of gel electrophoresis was not needed for subsequent samples, as these were directly analyzed with the colorimetric method, which also allows checking for possible interfering cross-contaminations of PCR reactions. The PCR products were separated on an 18% denaturing polyacrylamide gel, run in Tris-Borate-EDTA (TBE) buffer, they were then stained with SYBR® Gold (Thermo Fisher Scientific) and imaged on a Typhoon scanner (GE Healthcare).
Colorimetric discrimination of single nucleotide polymorphisms (SNPs) Colorimetric discrimination of the SNPs was performed according to a previously reported protocol 19 with some optimizations. Briefly, each amplification product related to each SNP was analyzed in duplicate, in order to have a nakedeye discrimination of the results after hybridization with the two discriminating probes (Integrated DNA Technologies, Inc.) matching the two possible allele variants (Table S1 †). For each duplicate, a quantity of 15 μL of the biotinylated PCR product was incubated for 5 min at room temperature with 5 μL of Dynabeads M-280 Streptavidin (Invitrogen) paramagnetic microparticles, preventively washed twice with hybridization buffer (HB) (1× PBS pH 7.4, 5% w/v PEG 600) and resuspended in an equal volume of HB. The microparticles coated with biotinylated PCR products were incubated with 50 μL of 0.15 M NaOH for 5 min to denature the double-stranded amplification product and magnetically separate the interfering non-biotinylated strand. The suspension was washed with 100 μL of 0.15 M NaOH and subsequently with 100 μL of HB, and finally it was resuspended in 20 μL of HB. The two samples were incubated for 5 min at room temperature with 10 μL of a 10 μM solution of one of the two discriminating probes. Then, 300 fmol of AuNPs conjugated with detection probes was added and kept at room temperature for 10 min. After a magnetic wash with 100 μL of HB at room temperature, the mixtures were incubated at 52°C for 8 min and then resuspended in 25 μL of HB for the readout of the colorimetric result.
Determination of the HLA C*06:02 allele HLA-C*06:02 genotyping was determined with the SNP tagging method described by Nikamo

Reverse dot blotting
The HLA typing tests were performed by reverse hybridization with the INNO-LiPA HLA-C kit (Fujirebio, Gent, Belgium). The test was performed using the Auto-LiPA automated test system, and the results were analysed with the LiRAS software (Fujirebio, Gent, Belgium). Briefly, biotinylated amplified targets, were chemically denatured, and then hybridized with membrane-immobilized probes. Afterwards, stringent washing removed mismatched hybridization. Then, streptavidin-conjugated alkaline phosphatase was added, followed by a chromogenic substrate for this enzyme. A purple/brown coloration of a probe band indicated perfect match of that specific probe with the target.

Automatic analysis of the dataset
Automatic analysis of the dataset was performed on a photograph taken with a smartphone, exploiting established image processing tools, in two steps. First the contour of each well was detected as a circle using the circular Hough transform 33 and, subsequently, a clustering step based on the internal color of each circle was applied to automatically group the data. In the second step, the RGB pixel values were converted to the YUV color space, because it is more similar to the average human perception. Then, the k-means clustering algorithm 34 was used to separate the wells, selecting the average YUV color within the wells as a descriptive feature. k-Means is an iterative data-partitioning algorithm that assigns n observations to exactly one of the k clusters, where k is chosen before the algorithm starts. In the case of this dataset, n is 300, equal to sample cardinality (i.e. the number of wells) and k is 2. A warning label (blue), which requires a human control, was used in the case of rare light red-orange points wrongly clustered as yellow (because the RGB orange value is nearer to the yellow level than the red one). The k-mean algorithm has the advantage of finding two groups without any prior knowledge of the RGB level of yellow and red (so it can work with any two colors under any light conditions).


extfieldsofstudy: Biology Medicine

text:
Identification of a novel deletion in the MMAA gene in two Iranian siblings with vitamin B12-responsive methylmalonic acidemia

Background Adenosylcobalamin (vitamin B12) is a coenzyme required for the activity of methylmalonyl-CoA mutase. Defects in this enzyme are a cause of methylmalonic acidemia (MMA). Methylmalonic acidemia, cblA type, is an inborn error of vitamin B12 metabolism that occurs due to mutations in the MMAA gene. MMAA encodes the enzyme which is involved in translocation of cobalamin into the mitochondria. Methods One family with two MMA-affected children, one unaffected child, and their parents were studied. The two affected children were diagnosed by urine organic acid analysis using gas chromatography-mass spectrometry. MMAA was analyzed by PCR and sequencing of its coding region. Results A homozygous deletion in exon 4 of MMAA, c.674delA, was found in both affected children. This deletion causes a nucleotide frame shift resulting in a change from asparagine to methionine at amino acid 225 (p.N225M) and a truncated protein which loses the ArgK conserved domain site. mRNA expression analysis of MMAA confirmed these results. Conclusion We demonstrate that the deletion in exon 4 of the MMAA gene (c.674 delA) is a pathogenic allele via a nucleotide frame shift resulting in a stop codon and termination of protein synthesis 38 nucleotides (12 amino acids) downstream of the deletion.

adenosylcobalamin (AdoCbl) and methylcobalamin (MeCbl). AdoCbl is required for activity of the mitochondrial enzyme methylmalonyl-CoA mutase (MCM, EC 5.3.99.2), which catalyzes the conversion of methylmalonyl-CoA to succinyl-CoA during catabolism of branched-chain amino acids and odd-chain fatty acids. Disorders that impair synthesis of AdoCbl cause accumulation of methylmalonic acid in urine and blood due to decreased activity of methylmalonyl-CoA mutase, which leads to methylmalonic acidemia (MMA) [2].
Individuals affected with this disorder present with lethargy, failure to thrive, recurrent vomiting, dehydration, respiratory distress, and hypotonia in the first year of life [2].
The gene responsible for cblA has been identified as MMAA [3]. MMAA is located on chromosome 4q31.1-2, and consists of seven exons that encode a protein containing 418 amino acids. Two different enzymatic functions have been identified for the MMAA gene product: a role in vitamin B12 transport into the mitochondria, and the conservation or re-activation of MCM [4]. In this study, we report the clinical, biochemical, and genetic defects and structural analysis of MMAA in one Iranian family with isolated MMA.We demonstrate that this novel deletion in exon 4 of MMAA is a pathogenic allele that results in a premature stop codon and subsequent truncated protein and the loss of the important binding sites of the protein due to loss of the ArgK (PF03308) conserved domain.

Patients and MMA diagnosis
One family with two MMA-affected children, one unaffected child, and their parents were studied. A daughter who died five days after birth was not assessed. The two affected children were diagnosed after presenting with clinical symptoms and high levels of methylmalonic acid by urine organic acid analysis. All procedures followed were in accordance with the ethical standards of the committee on human experimentation of MUMS (Mashhad University of Medical Science). Informed consent was obtained from all subjects for being included in the study.

Structural and functional analysis
The novel mutation (p.N225M) was considered for the purpose of structural modeling in order to completely understand the vast structural damages caused by truncation of the protein after the occurrence of the deduced novel mutation in comparison with the normal structure. Although the protein loses a large part of its structure, this structural modeling would help to better understand the vast dysfunctioning of the mutated protein by visually analyzing it. To accomplish it, the template most similar to the deduced normal protein with 95.7 % sequence similarity (Swiss Model Template ID (SMTL ID): 2www.1) was chosen and used to build the tertiary structure of the normal and mutated proteins by using the online Swiss-Prot server for automated modeling [6]. The 2www.1 structure has been solved to 2.64 Å using the X-ray diffraction method [7]. The structures were then considered for the purpose of the energy minimization task using ZMM software. ZMM software uses amber-all-atom force field with the cut-off distance of 10 Å to minimize the conformational energy in the space of generalized coordinates [8]. Low energy conformation was reached after 100 sequential tasks of energy minimization failed to actually enhance the quality of the assembled structures [9].
The resulting models were finally analyzed in terms of essential accuracy aspects such as G-Factor, bond length and bond angles by the help of PROCHECK [10] and WHAT-IF [11] online programs for of which all the results satisfied the important criteria. Moreover, most of the residues were inside of the favorable regions of the Ramachandran map.
In addition to revealing the structural defects, the effects of this novel mutation (p.N225M) on the function of the deduced protein was further investigated by considering its mutational role on the function of the conserved ArgK (PF03308) domain which starts at LUE p.101 and ends at THR p.384.

Clinical symptoms and biochemical findings
This novel deletion was identified in two children: a 4-year-old girl with metabolic problems, and her 1.5-year-old brother. Both children had signs and symptoms of MMA. The family pedigree is shown in Fig. 1a. The clinical symptoms and biochemical findings of the two affected children are presented below: The girl, who was a product of a consanguineous marriage, was born at term following a normal pregnancy. She was admitted to the hospital at 2 months of age with a history of persistent vomiting and progressive lethargy. Pyruvate was 1.49 mg/dL (normal, 0.3-0.7), ammonia was 197.60 mmol/L (normal, 64-107), lactate was 42.1 mg/dL (normal, 0 to 20), venous pH was 7.48 (normal, 7.35-7.45), and bicarbonate was 21.9 mmol/L (normal, 22-28). Tandem mass spectrometry analysis showed significant elevation of propionylcarnitine (13 μmol/L), which may be indicative of MMA or propionic acidemia (PA). Urine organic acid analysis indicated abnormally high excretion of methylmalonic acids (1250 mmol/mol creatinine). She was diagnosed with MMA and treated with protein restriction, daily intramuscular hydroxycobalamin injections (1.0 mg), Shohl's solution (4.5 ml/day), benzoate tablets (250 mg/day), L-carnitine (5 ml/ day), and metronidazole syrup (4.5 ml/day). At 4 months she presented with lethargy, vomiting, hepatomegaly, hypotonia, anorexia, developmental delay, and failure to thrive. At 10 months she had begun to achieve sufficient weight and showed developmental and hypotonia improvement with protein restriction and L-carnitine and hydroxycobalamin injections. In the third year of treatment, her urinary methylmalonic acid excretion was 12 mmol/mol creatinine, but development and movement were delayed compared with normal individuals.
The boy's birth history and development in the first year of life were normal and no symptoms of MMA were observed. He was hospitalized at 1.5 years of age for persistent vomiting and infection. Glucose was 250 mg/dL (normal, 70-105), venous pH was 7.07 (normal, 7.35-7.45), bicarbonate was 6.5 mmol/L (normal, 22-28), blood urea was 55 mg/dL, and creatinine was 0.5 mg/dL. He was treated for MMA with L-carnitine, sodium benzoate, and vitamin B12. Dialysis was initiated, but he died on the second day of hospitalization.

Novel MMAA deletion and mRNA expression analysis
The family pedigree and sequence chromatograms for exon 4 of MMAA from the five family members are shown in Fig. 1. A homozygous deletion in exon 4 of MMAA, c.674delA, was found in both the affected son and daughter. An adenine deletion was detected at position c.674 in exon 4. This deletion predicted a nucleotide frame shift and amino acid change from asparagine to methionine at position 225 (p.N225M), resulting in a premature TAA stop codon at position 238 and truncation of the protein 12 amino acids downstream of the methionine produced by the frame shift. We confirmed the mutation's homozygosity through analyses of the parents' and unaffected son's samples. We found one heterozygous variant in exon 4 of MMAA, c.674delA. This mutation has been recorded in the GenBank database with accession number KR026958. These results, from two affected and one unaffected siblings and their parents, confirmed the diagnosis of cblA-deficient MMA. MMAA mRNA expression analysis identified the expected truncated message. The mRNA expression analysis and sequencing chromatogram of expression analysis are shown in Fig. 2. The premature TAA stop codon 12 amino acids downstream of the methionine produced by the frame shift are shown in this figure. Mutational effects on the structural defects and dysfunctioning of the protein As discussed above and shown in Fig. 3, the novel mutation (p.N225M) is the underlying cause of a big frame shift and the loss of a large part of the deduced novel protein due to the early truncation which results in conformational instability. The modeling of the deduced novel mutation also reveals the structural importance of the cutting regions in dysfunctioning of the MMAA novel mutated protein. The MMAA protein contains the ArgK (PF03308) conserved domain which starts at p.101 and ends at p.384. Therefore a very large part of the protein lost (from p.225 to the end) contains the elements of this conserved domain site. The function of this domain is not completely known yet. However, it may have GTPase and ATPase activity and may also contain important binding sites of the MMAA protein due to its roles and conservation history. Moreover, the MMAA protein is a cobalamin cofactor which is considered to have a confident (vastly reported in the literature: data collected from www.string-db.org at String: 9606.ENSP00000281317) interaction with MUT protein. The MUT gene has a B12binding site which can bind two different types of cobalamin cofactors, methylcobalamin and adenosylcobalamin, like the one described here [12][13][14]. It can thus be assumed that the loss of a large part of the MMAA protein (involving the conserved domain ArgK (PF03308) due to the novel mutation p.N225M) and the alteration of twelve amino acids (LYS p.226 to THR p.237) which happened due to the frame shift caused by c.674delA may drastically affect the binding activity of these two proteins. This effect can be considered as either removing the complete interaction sites or reducing the binding parts between them in such a way that the resulting proteins may fail to interact completely. Therefore in both situations the resulting interaction would not be completely functional.

Discussion
In 1968, Rosenberg and Lindblad described infants with severe metabolic ketoacidosis and developmental delay with abnormally high levels of methylmalonate in their blood and urine. Treatment of these infants with AdoCbl reduced both the severity of their symptoms and methylmalonate excretion [15,16]. Further studies indicated that these infants expressed normal methylmalonyl CoA mutase, but synthesis of AdoCbl was impaired [17,18]. Kaye et al. (1974), reported two patients with MMA who expressed normal methylmalonyl CoA mutase, but were defective in AdoCbl synthesis and unresponsive to cyanocobalamin treatment. Biochemical and genetic analyses showed that cblB was also defective in these patients [19].
Since then, more than 45 different disease-causing mutations have been identified in MMAA, which include 30 missense or nonsense mutations, four splicing mutations, five small deletions, five small insertions, and one gross deletion. All the mutations are available on the Human Gene Mutation Database (HGMD) [20][21][22][23][24]. It is also very interesting to know that the MMAA protein has a low complexity region which starts at p.213 and ends at p.224, an amino acid before the reported novel mutation (p.N225M). Low complexity regions are prone to mutation as described previously [25] and in our case; the mutation occurred immediately after the end of this region. This low complexity region is demonstrated in both normal and novel mutated (p.N225M) MMAA protein in Fig. 3 above. Molecular analysis has not been performed on MMAA in the Iranian population. Early detection of MMA is critical; the girl in this family who was diagnosed at 2 months has survived, while the boy, who was diagnosed at 1.5 years, died soon after. Early diagnosis can prevent metabolic crises and mental retardation. Identification of mutations in MMA patients aids prenatal diagnoses and determination of the parents' and relatives' carrier status.


Lethal gene drive selects inbreeding

Gene drive has come of age with advances in genetic engineering. An unresolved question is whether gene drive systems will achieve population suppression or will instead succumb to resistance in the target population. Here is it shown that the evolution of inbreeding can be favored in response to a gene drive system, substantially reducing its efficacy.

INTRODUCTION
The proposed use of selfish genes to suppress or extinguish populations is at least half a century old [1][2][3], but the feasibility of actually engineering selfish genes is new. There has thus been much excitement about the possibility of using these approaches to eradicate disease vectors, balanced by concerns about the possibility of unforeseen harm. Perhaps the most tangible approach is one outlined by [4], of creating a homing endonuclease gene (HEG) that inserts itself into an essential gene. Under the idealized assumptions of 100% segregation distortion in gametes of heterozygotes (germ line only), normal heterozygote viability and fertility but homozygote lethality, such a selfish gene is expected to evolve to such an extreme as to cause a 50% reduction in population fecundity if the original research article segregation distortion is limited to one sex [5,6]. A segregation distortion that operates in both sexes can evolve to fixation and death of all progeny, ensuring extinction [4,6,7]. A 50% reduction in population fecundity might seem to be so inconsequential as to eliminate it from further consideration. However, the 50% reduction applies to a single HEG, the combination of multiple HEGs at different locations in the genome ('stacking') has been proposed as a way of suppressing population fecundity to arbitrary levels [4].
The HEG need not work as completely as expected. Extreme levels of population suppression from the HEG are sensitive to even minor variations in parameter values [8,9]. More importantly, an HEG that harms population fitness will select resistance mechanisms [10,11]. Since HEGs target specific DNA sequences, the most obvious form of resistance to the HEG is a change in the target sequence so that the HEG will no longer duplicate itself in heterozygotes [4]. Resistance could also take the form of interfering with endonuclease expression or functionality. The problem of target sequence evolution has been countered with the suggestion of deploying multiple HEGs simultaneously [4], but other resistance mechanisms would not obviously be thwarted by that approach.
Here I address another possible evolutionary mechanism that may interfere with the spread and long term maintenance of a lethal HEG: evolution of inbreeding. It is appreciated that inbreeding reduces the population impact of 'lethal' HEGs and other gene drive systems [3,4,12,13]. What is not clear is whether inbreeding is actually favored and how much it rescues mean fitness once a lethal HEG has invaded the population. Although a fixed level of inbreeding should reduce the incidence of the recessively lethal HEG, an allele that increases the level of inbreeding will itself suffer increased loss from any excess inviable progeny that it creates, perhaps selecting against inbreeding and even favoring increased outcrossing. It will in fact be shown here that inbreeding does evolve under some conditions, but the extent to which population fitness recovers depends heavily on the magnitude of inbreeding depression. Furthermore, the level of inbreeding that evolves is often not the level that would maximize population fitness were inbreeding imposed on the population.
The models analyzed here incorporate two major simplistic assumptions: the evolution of inbreeding is treated as the evolution of self-fertilization in an infinite population of simultaneous hermaphrodites with no spatial structure, and inbreeding depression is treated as a static quantity. Analysis of these models seems a reasonable first step in deciding whether the problem justifies inclusion of greater reality.

THE MODELS
Accommodating inbreeding as selfing All models assume a population of diploid hermaphrodites capable of a mix of outcrossing and self-fertilization (selfing). Each individual produces a constant amount of sperm and of eggs. Eggs can be fertilized either by sperm chosen randomly from an outcross pool or by self sperm, from the individual who produced the ova.
The models assume that selfed offspring have a fitness lower than that of outcrossed offspring ( < 1), with inbreeding depression parameter d ¼ 1 À s: In this initial study, inbreeding depression is invariant throughout the evolutionary process. In real systems, inbreeding depression is often partially purged upon extended inbreeding, but allowing inbreeding depression to be static is a reasonable starting point and, if anything, provides a conservative measure of the vulnerability of gene drive systems to be suppressed by inbreeding (see Discussion section).

Genotypes and phenotypes
The models assume two unlinked loci, each with two alleles (A/a, D/d). An individual's level of selfing is controlled by its genotype at the A/a locus independently of the genotype at the other locus ( Table 1). The D/d locus affects viability and experiences gametic drive. Specifically D is a recessive lethal that enjoys complete segregation distortion in heterozygotes-Dd produces all D gametes in spermatogenesis (and also in ovogenesis in some models). For most analyses, Dd is considered to be of normal viability and fertility (the drive does not operate in somatic tissues so those remain heterozygous and function normally), some analyses instead impart a slight fitness disadvantage to the heterozygote. The DD genotype of both sexes dies at conception, so viable genotypes at this locus are limited to dd and Dd.

Modifications of reproductive output
The net reproductive output of a genotype that produces selfed offspring or experiences zygote loss can be modeled in different ways, each of which may be observed in nature. For example, there are two extremes that span the possibilities for the impact of selfing on sperm contributed to the outcrossed pool [14,15]. At one extreme, sperm contributed by an individual to the outcross pool may be reduced in proportion to the fraction of ova selfed by that individual (sperm are 'discounted'). Alternatively, the outcross contribution may be unaffected (sperm not discounted) because few sperm are actually used for selfing compared with the huge number released for outcrossing. Both cases will be evaluated here.
In species with post-zygotic investment in offspring, reproductive compensation may also operate for lost zygotes, whereby the total number of viable offspring released from the maternal parent is larger than expected from the fraction of viable zygotes [15,16]. Here, the models are limited to maternal offspring production in direct proportion to viable zygotes.

Four models
The foregoing sections identified two important biological variations that affect the evolutionary dynamics: (a) Drive operates just in male gamete production or in both male and female gamete production (b) The contribution of sperm to the outcrossed pool is discounted in proportion to the fraction of ova selfed, or alternatively, selfing has no effect on the contribution to the outcross pool.
By considering these variables in all combinations, there are four models.
Interest is in whether selfing evolves specifically in response to the presence of D and the load from DD inviability. Appropriate analyses are thus limited to parameter values in which selfing would not evolve if D was absent. It is now well known that the evolution of selfing is sensitive to both the magnitude of inbreeding depression and whether selfed sperm are discounted [15,17,18]. When selfed sperm are not discounted from the outcross pool, there is an extra benefit to male function. Consequently, selfing is intrinsically beneficial at low values of inbreeding depression and is favored until inbreeding depression exceeds 0.5, the appropriate range for our problem is thus > 0.5. If instead selfed sperm are discounted, then selfing is favored only if inbred offspring are more fit than outcrossed offspring, so the appropriate range for our problem includes even small values of inbreeding depression, > 0.
The models will be analyzed for gene frequency evolution and mean fitness in the population (measured as viable offspring), equations are given in the Appendix. Applied interest in these evolutionary processes is primarily for population control and the potential for extinction [4,[19][20][21]. Yet, as has long been realized in applications of the sterile insect technique (used to suppress target pest populations), the impact of a particular intervention on adult population size depends heavily on ecology, which is often species-specific. Thus the analyses here omit any translation between mean fitness and population suppression.

Expectations
The evolution of gene drive has different consequences depending on whether drive operates in one sex (males, here) or both sexes. A drive allele with 100% distortion in sperm and ova of the same Dd heterozygote can evolve to the extreme that all offspring in the population are inviable DD [7]. This outcome ensures population extinction. Drive operating only in one sex can only evolve to the point that half the offspring are DD, half are Dd [5,6] with a 50% reduction in population fitness (fecundity). A minimal expectation is thus that selfing should be favored for higher levels of inbreeding depression when drive operates in both sexes than when it operates just in one sex. For example, one might anticipate that selfing would not evolve in response to male-drive when inbreeding depression exceeds AADd s 2 0.5 ( < 0.5) because mean fitness in the absence of selfing is 0.5.

Evolutionary properties of single cases
When accounting for sex differences in segregation distortion, sperm discounting, inbreeding depression levels and initial conditions, there are many combinations of models to consider. To enable an easy comparison across this variety, a summary analysis of cases will be presented below. Before progressing to this summary, a few cases will be explored in detail.
The first case will assume (i) gene drive in males only, (ii) no sperm discounting (no sperm lost to selfing), and (iii) inbreeding depression of = 0.51, just barely large enough to prevent the evolution of selfing in the absence of gene drive.

Effect of selfing on mean fitness
As a first step, mean fitness was calculated under the assumption that all individuals experienced the same level of selfing. In the absence of the HEG (allele D), mean fitness declines linearly with inbreeding, with a slope of À (blue line, Fig. 1A). With the HEG present and the system allowed to equilibrate, equilibrium fitness rises from 0.5 in the total absence of selfing to a maximum of 0.66 at a selfing level of 0.67 (yellow curve, Fig. 1B). This relationship appears linear up to this point. The HEG is maintained in populations with selfing rates less up to and including 0.67 but is lost in populations with levels of selfing of 0.68 and higher. With the HEG absent, the HEG-selfing mean fitness coincides with the inbreeding-fitness line and descends with yet higher imposed levels of selfing.
Genetic variation in selfing: invasion and equilibrium fitness The preceding analyses did not allow for the evolution of selfing, merely imposing a fixed level on the entire population. It might be inferred that, because mean fitness increases with selfing in the presence of the HEG (up to a point), selfing will evolve. But this evolution need not happen, or it may happen in unintuitive ways. An easy first step in evaluating evolution is to consider whether a selfing allele will invade a fully outcrossing population with the HEG at equilibrium. Equations for male drive, no sperm discounting (Appendix) were linearized for a rare allele A that increased selfing in a fully outcrossed population. Linearization assumed the equilibrium at which the HEG had gone to fixation with aadd absent, aaDd fixed.) For the genotype frequency vector ðAadd; AAdd; AaDd; AADdÞ 0 , the linearized transition matrix is with leading eigenvalue (for the parameter range of interest) Recall that s ¼ 1 À d. A combination of analytical and numerical investigations of (2) reveals that selfing cannot evolve if does not exceed 0.25 (if ! 0.75). Regions of the (s 1 , s 2 ) space where eigenvalue (2) exceeds 1.0 are plotted in blue in Fig. 2A for the case of = 0.51.
In addition to the invasion analysis, numerical analyses of the full equations were conducted, providing mean fitness and average selfing rates at equilibrium. Evolved mean fitness values are given in Fig. 2B, the diameter of each blue dot is proportional to the excess of mean fitness above the baseline 0.5. It is evident from both the invasion analysis and the equilibrium analysis that evolution of selfing requires that the A allele enacts high levels of selfing in either or both Aa and AA genotypes. Small levels of selfing do not invade, and as seen in Fig. 2B, the highest equilibrium fitnesses occur with s 2 at or near 1.0. Selfing does not evolve in this system incrementally. Fig. 2C gives the average level of selfing at equilibrium, with the diameter of the dot proportional to the selfing rate above 0.01. The highest levels of selfing occur in the middle of the parameter space, but these do not attain the highest fitness levels (although fitness levels in this region are not much below maximum).
It is easily inferred from Fig. 2B and C that (i) invasion of the selfing allele reduces the frequency of the drive allele B from its level with pure outbreeding, and (ii) invasion of the selfing allele does not usually lead to its fixation (because the average selfing rate in the population is often less than the value of s 2 ). The former is easy to understand and plausibly follows from the analysis in Fig. 1B. The lack of fixation of selfing alleles presumably stems from a balance between the negative effect of inbreeding depression accompanying selfing and the benefit of selfing purging the drive allele, but an  intuitive explanation for the specific outcome in individual cases seems out of reach.
To gain further insight to the interlocus dynamics, a correlation of the number of selfing alleles with the number of drive alleles within a zygote was calculated across all viable zygotes. Across many trials, the correlation at equilibrium was invariably negative, typically in the range of À0.2 to À0.3 when mean fitness was 0.6 or higher. This result matches the intuition offered above: the selfing allele experiences a diminished association with the drive allele presumably by purging it (see below for an exception).
For comparison, the same model was analyzed when the inbreeding depression parameter was set to = 0.7. Note that this value is near the expected boundary of 0.75, above which selfing cannot invade. The invasion region of the (s 1 , s 2 ) space and equilibrium mean fitnesses are shown in Fig. 3. Average selfing rates are not shown because the values are too small to be visible using the same scale as in Fig. 2C.
A similar pattern was evident for male drive with sperm discounting (not shown). Invasion of the selfing allele A occurred for a much wider range of (s 1 , s 2 ) values for = 0.01 than for = 0.45, again, values near (s 1 ¼ 0; s 2 ¼ 0) did not invade. However, for both models of 2-sex drive, invasion occurred for even the smallest selfing values tested (s 1 ¼ 0:01; s 2 ¼ 0:01).
The male-drive models have another interesting property regarding the evolution of selfing. At equilibrium with 100% segregation distortion and all individuals Dd, a mutant A allele will always arise in the Dd background, yielding a AaDd genotype. Henceforth, there will be no way for A to become paired with the dd genotype because all outcross and self pollen it receives are DÀA is permanently trapped with D. Complete coupling of A with D ensures that the increased selfing imposes a load on A from the formation of DD homozygotes. Allele A cannot escape this load and is selected against relative to allele a. This outcome is, no doubt, a pathology of the strict parameter values that can be imposed in the model but would not likely apply in nature: The trapping of the selfing allele disappears if segregation distortion is less than complete, or if the starting population is not 100% Dd, whence Aadd is then formed and selfing can perpetuate it (and create AAdd). Nonetheless, the result reinforces intuition about the different processes affecting the selfing allele.
Broad patterns in the joint evolution of selfing and gene drive The four classes of models were analyzed for a total of four initial conditions across 120 different (s 1 , s 2 ) values (Fig. 4). For each model, two different inbreeding depression values () were chosen within the feasible range, panels (E) and (F) are from the analyses in Figs. 1, 2 and 3.
Each panel of Fig. 4 combines the two types of fitness curves of Fig. 1, again using blue for fitness in the absence of the HEG and (dashed) yellow when the HEG is present. The slopes and intercepts of the lines vary among panels, but the patterns are the same whereby the yellow curve rises to intersect the blue, and they then decline together at higher selfing levels. (In panel F, the yellow curve actually declines across the entire spectrum of selfing levels due to the high level of inbreeding depression.) For 2-sex drive, the yellow curves intersect the origin, reflecting the death of all individuals at equilibrium in a fully outcrossing population.
The numerical studies of Figs. 2 and 3 reveal a plethora of evolved mean fitness and average selfing rates across the spectrum of possible s 1 and s 2 selfing rates. Instead of presenting similar comprehensive data for each of the 6 other panels in Fig. 4, the panels show merely the highest mean fitness values and associated mean selfing levels obtained across the set of (s 1 , s 2 ) values and initial conditions tested (black triangles). Where similarly high fitnesses were obtained across a range of average selfing rates, the triangles include those outcomes. The fitness values included were typically limited to within 0.02 of the observed maximum, with the goal of showing the span of average selfing rates giving rise to similarly high fitnesses.
There are several properties of these evolutions that defy intuition. The black triangles often do not coincide with maximum fitness on the yellow curve, and in a few cases they even lie moderately off the yellow curves. In Fig. 4E, the selfing alleles that attain maximal fitness are largely recessive, whereas those attaining maximal fitness in F are largely dominant, the two cases differ only in . Furthermore, it's not even clear in Fig. F why selfing is favored, as mean fitness declines when selfing is imposed on the population (yellow curve). In some cases, the highest fitness spanned a considerable range in average selfing rates (e.g., Fig. 4C and D), whereas in other cases, they were limited to a tiny range [e.g., in both Fig. 4B and Fig. 4F, the three triangles fall on top of each other]. A more comprehensive analytical approach will be required to resolve these puzzles.
A summary of the outcomes is given in Fig. 5. This is a plot of the maximum fitness observed versus the fitness of inbred offspring () in each of the eight panels of Fig. 4. Each of the eight points is coded to indicate one-sex drive versus two-sex drive ('1' versus '2') and whether sperm are discounted (triangles) or not (circles).
From the two preceding figures, several points are noteworthy.
(i) Given appropriate selfing rates for the Aa and AA genotypes, selfing evolved and increased mean fitness at least slightly above the fitness evolved with the HEG alone. Evolved fitness was equal to the fitness of selfed offspring () for the sperm discounting models but somewhat above for the models with no sperm discounting. (ii) Mean fitness did not necessarily attain the maximum that could be obtained if selfing was imposed on the population. (iii) In the male-drive models, selfing never evolved if the fitness of selfed offspring () was too low. Invasion analyses indicated that needed to exceed 0.5 in the sperm discounted model and 0.25 in models where sperm were not discounted. (iv) The drive allele (D) was lost only when a selfing rate of 1 could evolve (except for cases of apparent floating point errorsee the following section).
Given that DD is lethal, the assumption that Dd is fully viable becomes questionable: fitnesses of recessive lethal heterozygotes are typically slightly below maximal [22]. The analyses of Fig. 4 were, therefore, conducted again but assigning a viability factor of 0.98 to all Dd genotypes. (A Dd produced by selfing had fitness 0.98.) Although quantitative effects of this fitness adjustment were observed, they were slight, and a parallel figure to that of 4 was effectively indistinguishable (not shown).
A caution: the practical loss of drive for some initial conditions In some models, certain initial conditions led to trajectories with vanishingly small frequencies of the D allele (e.g. 10 -50 or less). When such low frequencies occur, the deterministic retention of D likely becomes irrelevant-it would be lost in any real population. The possibility of stochastic loss of D means that any implementation of a gene drive system should be analyzed stochastically for initial conditions that might be used.
For illustration, consider the model and value in Fig. 4A for s 1 = s 2 = 0.9. The presumed deterministic trajectory is shown in Fig. 6 for a specific set of initial frequencies. At first, D increases rapidly, in turn selecting increases in A. The rise in A and consequent increase in overall selfing drives D down to very low levels, less than 10 -30 for over 250 generations. However, A declines during this period, allowing D A plot of average selfing rates corresponding to that of Fig. 2C is not shown because it appears empty-the maximum observed selfing rate across the space is 0.028, too small to be visible as a point  The blue and yellow curves represent mean fitnesses at equilibrium for a uniform selfing rate (with and without the HEG), as in Fig. 1. The black triangles represent equilibria when selfing was allowed to evolve in the presence of the HEG, showing outcomes at or near the highest mean fitness observed across hundreds of trials with different selfing genotypes (selfing rate in the graph is to rebound, leading again to its rise and to further, damped cycles with A. The end point is an equilibrium with both D and A. D would be prone to extinction during these nadirs, and indeed, D was actually lost in some numerical trials (with slightly different initial frequencies than used in Fig. 6), until the floating point precision of the trial was increased.
(That such losses were due to floating point error was independently confirmed by M. Edgington, personal communication.) Anomalous loss of D was also seen in some numerical trials of other models. Again it seemed that the losses were not deterministic but instead resulted from inadequate floating point precision to maintain the infinitesimally low frequencies of D.

Extended evolution of selfing
The models assumed a single, bi-allelic locus for selfing with an initial state of outcrossing (s 0 = 0). It is possible that models allowing an extended evolution of selfing would lead to higher selfing rates with a higher mean fitness. Toward this end, further analysis was carried out of some trials in Fig. 4 in which the A allele had fixed and left the population with an intermediate level of selfing (D also present). The fixed A endpoint was then used to seed trials of the 2-locus model in which the evolved selfing rate was the baseline upon which further variation could act. Fig. 4B, D and E were evaluated: when a wide spectrum of selfing rates was tested against the new baseline, the opportunity for extended evolution often led to yet higher selfing rates and higher fitness above the starting state. Yet the extended evolution   did not exceed the highest mean fitness found in the collective set of trails in Fig. 4 (black triangles), although extended evolution often achieved the same highest fitness, or nearly so.
Models restricted to two alleles at one selfing locus may also not provide full insight to long term evolutionary dynamics. A 3-locus model was created to investigate this possibility (two loci-A/a and B/b-affected selfing, the D/d locus encoded the drive).
The case of Fig. 4H was evaluated, as the mean fitness evolved under the 2-locus model was substantially below the maximum possible in the presence of the HEG. The initial conditions specified selfing rates of 0, 0.9 and 0.9 for the aadd, Aadd and AAdd genotypes, respectively, as led to the modest invasion of selfing in Fig. 4H. The B allele forced increases in the lowest level of selfing but maintained 0.9 as the highest level, thus evolution of B would have forced the baseline level of selfing above 0. None of those parameter combinations or the others tested resulted in an increase of B, suggesting that a multi-locus selfing system does not easily evolve to the level of selfing that would maximize mean fitness in the presence of the HEG.
Taken together, these results suggest that the trials for Fig. 4 captured the highest fitnesses capable of evolving under selfing (given a fixed value of inbreeding depression). There are clearly constraints on which selfing alleles are favored (e.g. Fig. 2), but the full suite of alleles considered may have captured the maximal fitness evolution via selfing.

DISCUSSION
There is much justified excitement about the possibility of employing gene drive systems to limit wild populations of undesirable species [9,12,19,20,23]. In particular, homing endonuclease genes (HEGs) can now be developed that are recessive lethals but experience close to complete segregation distortion in heterozygotes [24]. Such HEGs can theoretically spread to fixation or to the point that the entire viable population is heterozygous [5][6][7], with a major reduction in mean fitness, possibly even extinction. Evolution of resistance to the HEG, or to its effects, becomes highly relevant in understanding the possible limitations of these engineered systems.
This study indicates that the final frequency of a recessive lethal enjoying complete segregation distortion can be reduced by the evolution of inbreeding. There is a consequent increase in mean fitness above that which evolves in the absence of selfing. One important result is that the fitness mitigation achieved by selfing is limited largely by the magnitude of inbreeding depression. In models assuming sperm discounting, selfing enabled mean fitness to avoid the low value expected from an uncontrolled HEG, but mean fitness at equilibrium was the same as that of a selfed offspring (s ¼ 1 À d).
Recovery was somewhat higher in the models with no sperm discounting (Fig. 5).
Although the deleterious population consequences of the HEG could be partly mitigated by the evolution of selfing, it was also true-in maledrive models-that selfing was favored only if the selfing allele enacted a sufficiently large degree of selfing. Alleles with low levels of selfing could not necessarily invade when alleles with high levels could. A conjecture is that an allele needs to attain high levels of selfing to purge the drive allele (which is a recessive lethal) and thereby escape the load that is otherwise associated with inbreeding. The deterministic trajectory of the drive allele may traverse very low frequencies, values so low that allelic extinction would be expected in real populations. Indeed, some numerical trials appeared to lose the drive allele simply from inadequate floating point precision. Gene drive releases will thus require extensive analyses to identify initial conditions that are least prone to unwanted loss of the drive allele, although those trajectories will also depend on the presumably unknown genetics of inbreeding variants in the target population.
Precedents in the evolution of mating systems It has been long appreciated that suppressors of recessive lethal meiotic drive systems are strongly favored [10,11]. However, the mechanisms of suppression may often be intricately tied to the molecular bases of the distorter [25]. The point of this paper is to suggest a type of suppression-inbreeding-whose origin and evolution is not sensitive to the molecular bases of the distortion. A resistance mechanism that is robust to the underlying mechanisms of distortion may thereby operate in many different systems.
The suggestion from this study is not without precedent, however. Evolution of inbreeding in response to a lethal HEG is one of now several examples of mating/genetic system evolution in response to a selfish allele (or alternatively, an altruistic allele). In some cases, it has merely been suggested that inbreeding limits the impact of a selfish element [e.g., 3,4]. In one study, however, inbreeding was shown to evolve as a response to an evolving altruism [26]: the allele for increased inbreeding became coupled with the allele for altruism, with inbreeding reinforcing the benefit of altruism by increasing the number of altruists interacting in family units. A few studies have shown that the number of sires-polyandry or monogamyevolves in response to the presence of a selfish element [27,28] or in response to altruism [29], the number of sires affects relatedness within families and thus affects the competition between selfish and non-selfish alleles. In a system of yet higher complexity, Lande and Wilkinson [30] found that a sex-linked segregation distorter could be partially or wholly suppressed by female preference of a sex-linked male trait if that trait was tightly linked to the segregation distorter. Brandvain and Coop [31] found that segregation distorters in female meiosis can select changes in the female recombination rate.

LIMITATIONS OF THE MODELS
Other forms of inbreeding With many proposed applications of gene drive systems to disease vectors and crop pests [4,20,21], the evolution of inbreeding as a possible 'escape' from a lethal HEG is especially relevant. The results here strictly apply to hermaphrodites (mostly plants), and even more narrowly to plants lacking self-incompatibility systems. There may well eventually be direct applications of these models in weed control, and the frequently observed evolution of selfing in plants [32] is a caution that such HEG implementations may have limited impact.
Yet the results of this study are also interesting for their implications to species with males and females, as in insects, where the most immediate uses of lethal HEGs and other genetic methods of control are entertained. For those species, selfing is not possible, so inbreeding would need to involve sib mating or some other localized mating structure. An obvious extension of the work here is thus to consider the evolution of inbreeding in those species. It might be anticipated that the results here will generalize to other forms of inbreeding, but further work is needed to go beyond mere speculation. Indeed, one potentially critical difference is that selfing more easily achieves high inbreeding coefficients than does sib mating. As suggested by R.
Lande (personal communication), the fact that selection favored selfing alleles only if they enacted moderate to high rates of selfing (e.g. Figs. 2 and 3) may indicate that the evolution of sib mating will face greater constraints than does the evolution of selfing. Such a result would raise hope that the evolution of inbreeding is less likely to thwart an HEG in insects than in selfing species. However, the restrictions against evolution of low selfing rates were observed only for the male-drive models, so it is not clear whether the evolution sib mating would be constrained with 2-sex drive.

Evolution of inbreeding depression
The models assumed that inbreeding depression () was fixed throughout the evolutionary process. In contrast, inbreeding depression is known to evolve and can be at least partly overcome when inbreeding is imposed on formerly outcrossing populations [33]. In the results here, the magnitude of inbreeding depression played a critical role in invasion and ultimate recovery of population fitness (e.g. Fig. 5). Any purging of inbreeding depression would not necessarily be relevant to invasion conditions, but it would improve the recovery after invasion beyond that seen in models with fixed . An important extension of these models is thus to incorporate the evolution of inbreeding depression with the evolution of inbreeding.
Recent work points toward ways in which dynamic inbreeding depression might be accommodated. A consensus is emerging that inbreeding depression is due largely to deleterious mutations, but there are two important classes of deleterious mutations that contribute and have different consequences for purging: weakly deleterious mutations with large additive effects and strongly deleterious mutations such as lethals whose effects are largely recessive [14,33,34]. The weakly deleterious mutations are abundant, whereas the strongly deleterious mutations are much less common. In the short term, purging is chiefly from loss of the strongly deleterious class [14,34]. A reasonable conjecture, therefore, is that results based on fixed inbreeding depression levels would accrue to inbreeding depression stemming from the weakly deleterious component, with the recessive lethal component being rapidly purged.
Another possible extension of the model is to expand the biological ramifications of inbreeding depression to include reproductive compensation for inviable zygotes. In species with post-zygotic parental investment in offspring, reproductive compensation boosts offspring number by replacing inviable genotypes with viable ones [16,15]. The ramifications of this could be to facilitate the evolution of drive by supplementing Dd offspring to replace inviable DD.

Implications for population suppression
The models analyzed here are of population genetics evolution. Yet interest in applying gene drive systems is often for population control-to limit numbers of adults. Unfortunately, across diverse ecological settings, there is no straightforward relationship between mean fitness and adult population size (the one exception being that a mean fitness of 0 ensures extinction).
Of particular relevance to population suppression by a recessive lethal HEG is that the lethality will typically operate at the zygotic stage. A reduction in zygotes will be especially prone to be overcome by many types of density-dependent regulation, whereby the number of adults is much less reduced than is the number of zygotes [35]. In the long history of work on the sterile insect technique-the goal of which is to suppress pest populations by artificially introducing sterile males or females-species-specific ecology often underlies the difference between success and failure [36][37][38][39]. It is thus to be expected that the impact of gene drive systems to suppress a target species will also be ecology dependent, unless of course it can destroy all progeny.
The models here merely indicate that selfing (and possibly other forms of inbreeding) will evolve if there is genetic variation for inbreeding and if inbreeding depression is not too severe. At least for mosquitoes-perhaps the most obvious target species for gene drive control efforts-there is hope that both constraints may obstruct the rapid evolution of inbreeding. There is a long history of inbreeding mosquitoes to establish laboratory stocks, they not uncommonly show strong inbreeding depression, and they also exhibit high heterozygosity despite inbreeding-suggestive of balanced lethals that will persist indefinitely as a load from inbreeding [40]. As to whether sib mating is even feasible in the absence of inbreeding depression, a survey of mosquito mating systems reveals that a swarming mating system is widespread in groups of diseasecarrying mosquitoes, suggesting that they are neither predisposed to sib mating nor to adopting a life style that could easily convert to sib mating [41]. Pupal mating is known in some distantly related groups, however.

Future
As noted above, an important direction motivated by the present study is extending the models to the organisms most likely to be targeted by lethal gene drive systems: mosquitoes and other insects, which have separate sexes. Such extensions would include sib mating instead of selfing, relevant ecology, and population structure. This level of modeling is obviously challenging, but models of gene drive spread with realistic ecology have already been developed [42][43][44], so adding the evolution of sib mating should be within reach.
The empirical feasibility of evolving high levels of inbreeding during assault by a recessive lethal HEG remains to be seen. We have essentially no field experience with gene drive systems. There is, however, over half a century of experience with various forms of sterile insect applications [45][46][47]. Sterile insect techniques almost universally rely on the release of lab-reared insects that, when mated with wild insects, cause death or sterility of the progeny [23]. The assault from sterilizing, lab-reared insects should also favor inbreeding as one of several mechanisms that avoid matings that produce sterile progeny. Yet resistance via assortative mating in wild populations subjected to the sterile insect technique has rarely been reported, despite many applications and successes (summarized in 13]. In this comparison, it may be important that the sterile insect technique relies on inundation of the wild population with lab-reared insects, the wild individuals becoming increasingly overwhelmed with sterilityinducing matings as the population declines. Gene drive systems may have the opposite effect, encouraging de facto consanguinity as the population density declines. Nonetheless, there are many reasons to be hopeful that gene drive systems will be able to achieve long-standing population control in at least some species.

references
Here, D i,j represents a diploid adult, with i A alleles and j D alleles. Diploids may have 0, 1 or 2 A alleles but only 0 or 1 D alleles because D is a recessive lethal. M k,n (F k,n ) represents outcrossed male (female) gametes with k A alleles and n D alleles. S i,j represents diploid offspring by selfing, with similar subscripting as for D i,j .
Recursion equations for the gametes and selfed offspring are specific to each model, as follows. It is assumed that all ova are fertilized, and sperm frequencies are normalized. Although the body of the paper restricts the analyses to complete drive, these equations use an unsubscripted D to represent the fraction of gametes of a D heterozygote carrying the D allele, facilitating the derivation. Selfing rates are s 0 for aa, s 1 for Aa, and s 2 for AA, regardless of the other locus.