成人做爰A片免费播放乱码,亚洲日韩中文字幕无码一区,国产亚洲精品久久久久秋,巨大黑人video

銷售熱線

13818239648
主營產(chǎn)品:試劑,耗材,血清,細(xì)胞,抗體
  • 技術(shù)文章ARTICLE

    您當(dāng)前的位置:首頁 > 技術(shù)文章 > Bridge RNAs direct programmable recombination of target and donor DNA

    Bridge RNAs direct programmable recombination of target and donor DNA

    發(fā)布時(shí)間: 2024-07-01  點(diǎn)擊次數(shù): 856次

    Bridge RNAs direct programmable recombination of target and donor DNA

    Nature volume 630pages984–993 (2024)Cite this article

    • 56k Accesses

    • 1041 Altmetric

    • Metricsdetails


    Abstract

    Genomic rearrangements, encompassing mutational changes in the genome such as insertions, deletions or inversions, are essential for genetic diversity. These rearrangements are typically orchestrated by enzymes that are involved in fundamental DNA repair processes, such as homologous recombination, or in the transposition of foreign genetic material by viruses and mobile genetic elements1,2. Here we report that IS110 insertion sequences, a family of minimal and autonomous mobile genetic elements, express a structured non-coding RNA that binds specifically to their encoded recombinase. This bridge RNA contains two internal loops encoding nucleotide stretches that base-pair with the target DNA and the donor DNA, which is the IS110 element itself. We demonstrate that the target-binding and donor-binding loops can be independently reprogrammed to direct sequence-specific recombination between two DNA molecules. This modularity enables the insertion of DNA into genomic target sites, as well as programmable DNA excision and inversion. The IS110 bridge recombination system expands the diversity of nucleic-acid-guided systems beyond CRISPR and RNA interference, offering a unified mechanism for the three fundamental DNA rearrangements—insertion, excision and inversion—that are required for genome design.

    Main

    Evolution has dedicated a vast number of enzymes to the task of rearranging and diversifying the genome. This process enables the emergence and functional specialization of new genes, the development of immunity3 and the opportunistic spread of viruses and mobile genetic elements (MGEs)1,2. MGEs are abundant throughout all domains of life and often mobilize through a transposase, integrase, homing endonuclease or recombinase. These enzymes typically recognize DNA through protein–DNA contacts and can be broadly classified by their target sequence specificity, which ranges from site-specific (for example, Cre and Bxb1 recombinases)4,5 to semi-random (for example, Tn5 and PiggyBac transposases)6,7.

    Insertion sequence (IS) elements are among the most minimal autonomous MGEs, and are found abundantly across bacteria and archaea. Many characterized IS elements use a self-encoded transposase that recognizes terminal inverted repeats (TIRs) through protein–DNA interactions8. IS elements have been categorized into approximately 28 families on the basis of their homology, architecture and transposition mechanisms, but they can be broadly grouped by the conserved catalytic residues of their encoded transposases. These include DDE, DEDD and HUH transposases, and, less frequently, serine or tyrosine transposases8.

    IS110 family elements are cut-and-paste MGEs that scarlessly excise themselves from the genome and generate a circular form as part of their transposition mechanism9,10. Given what is known about this mechanism and life cycle, IS110 transposases are more accurately described as recombinases. Although circular intermediates are found in other IS families, IS110 is the only family that uses a DEDD catalytic motif in its recombinase. The N-terminal DEDD domains of IS110 recombinases share homology with RuvC Holliday junction resolvases, suggesting that they have a unique mechanism of action compared with other IS elements. IS110 elements typically lack TIRs and appear to integrate in a sequence-specific manner, often targeting repetitive elements in microbial genomes11. Although the mechanism of DNA recognition and recombination for IS110 elements remains unclear, previous studies have suggested that the non-coding ends of the element flanking the recombinase ORF regulate recombinase expression12,13.

    Here we show that the IS110 circular form drives the expression of a non-coding RNA (ncRNA) with two distinct binding loops that separately recognize the IS110 DNA donor and its genomic insertion target site. By bridging the donor and target DNA molecules through direct base-pairing interactions, the bispecific bridge RNA facilitates DNA recombination by the IS110 recombinase. Each binding loop of the bridge RNA can be independently reprogrammed to bind and recombine diverse DNA sequences. We further show that this modularity enables a generalizable mechanism for DNA rearrangement through sequence-specific insertion, inversion and excision.

    IS621 recombinase binds to a ncRNA

    IS110 elements encode recombinases that are around 300–460 amino acids (aa) in length and have an N-terminal DEDD RuvC-like domain (Pfam ID: PF01548) and a C-terminal domain with a highly conserved serine residue8,14 (Pfam ID: PF02371) (Fig. 1a and Extended Data Fig. 1a,b). They use this recombinase to scarlessly excise out of their genomic context, yielding a double-stranded DNA (dsDNA) circular form that is inserted into specific genomic target sequences such as repetitive extragenic palindromic (REP) elements9,12,15,16 (Fig. 1b and Supplementary Table 1). Recombination of the circular form and the target centres around a short core sequence, and the intervening sequences between the cores and the recombinase coding sequence (CDS) are defined as the left (LE) and right (RE) non-coding ends. IS110 recombinases are highly diverse and widespread in prokaryotes, but only a small subset have been catalogued by curated databases or functionally characterized (Fig. 1c).

    Fig. 1: IS110 mobile genetic elements express a ncRNA that is bound by its encoded recombinase.
    figure 1

    a, Schematic representation of the IS110 recombinase protein sequence. b, Schematic representation of the structure and life cycle of an IS110 element. Core sequences are depicted as green diamonds, the genomic target site is shown in blue and the non-coding ends are orange. Sequences are from IS621. c, A midpoint-rooted phylogenetic tree constructed from 1,054 IS110 recombinase sequences. d, Distribution of non-coding end lengths across eight IS families. The maximum of the LE and RE lengths is plotted for each family. Box plots show median (centre line), interquartile range (IQR) (box edges) and 1.5 × IQR (whiskers). Outliers not shown. n?=?268 for IS110; n?=?18–184 for other families (Extended Data Fig. 2). e, Small RNA-seq coverage plot of the concatenated non-coding ends of IS621 and five related orthologues expressed from their endogenous promoter in E. coli. Top, sequence logo of the conservation of the σ70 promoter motif. TSS, transcription start site. f, MST of a fluorescently labelled IS621 recombinase with either WT or scrambled ncRNA to measure the equilibrium dissociation constant (KD). Mean?±?s.d. of three technical replicates. g, Consensus secondary structure of ncRNAs constructed from 103 IS110 LE sequences.

    Full size image

    We found that IS110s have the longest median non-coding end lengths, with a relatively narrow distribution, compared with other IS families (Fig. 1d, Extended Data Fig. 2). Upon excision, the circular form of the element reconstitutes a promoter across the core sequence of the concatenated RE–LE far upstream of the recombinase CDS12,13 (Fig. 1b), which suggests that a ncRNA could be expressed from this region. Previous reports have shown that the non-coding ends of IS200 and IS605 family elements are transcribed into RNAs that resemble CRISPR RNAs to guide endonuclease activity17,18, and small RNAs have been thought to modulate recombinase expression for the IS110 family member ISPpu9 (ref. 19).

    To investigate the potential presence of an IS110-encoded ncRNA, we focused on the IS110 family member IS621, which is native to some strains of Escherichia coli, and five closely related orthologues (Supplementary Table 2). Small RNA sequencing (RNA-seq) of E. coli containing a plasmid that encodes the concatenated RE–LE sequences of the predicted circular forms revealed a continuous peak spanning around 177?bp of the LE, starting from the predicted endogenous σ70-like promoter (Fig. 1e).

    Next, we measured the affinity of an in-vitro-transcribed 177-nucleotide (nt) ncRNA from IS621 and its purified cognate recombinase using microscale thermophoresis (MST). We found that the IS621 recombinase binds to the LE-encoded ncRNA, but not to a scrambled 177-nt RNA control, with high affinity (dissociation constant (KD)?=?2.1?±?0.2?nM) (Fig. 1f). Our data indicate that IS110 element excision reconstitutes a promoter to drive the expression of a ncRNA that specifically binds to its recombinase enzyme, suggesting that the ncRNA might have a role in recombination.

    ncRNA covaries with target and donor DNA

    We evaluated the ncRNA consensus secondary structure across 103 diverse orthologues, and revealed a 5′ stem-loop followed by two additional stem-loops with prominent internal loops (Fig. 1g and Extended Data Fig. 3a,b). The first internal loop has relatively low sequence conservation across orthologues, whereas the second is much more conserved (Extended Data Fig. 3c).

    We next asked whether the ncRNA might assist the recombinase in recognizing the target site or the donor DNA (that is, the IS110 element itself). To assess this, we systematically reconstructed the insertion sites and circular forms of thousands of IS110s (Fig. 2a). An iterative search using a custom structural covariance model of the IS621 ncRNA enabled the prediction of thousands of ncRNA orthologues encoded within LEs20 (Methods). We first created a paired alignment of IS110 ncRNAs with their respective target and donor sequences. To assess the possibility of base-pairing between the predicted ncRNAs and their target and donor sequences, we then performed a covariation analysis across 2,201 donor–ncRNA pairs and 5,511 target–ncRNA pairs. We overlaid a base-pairing concordance analysis to identify stretches of the ncRNA that might bind to either the top or the bottom strand of the target or donor DNA21 (Supplementary Data 1). Nucleotide sequence covariation would indicate evolutionary pressure to conserve base-pairing interactions between ncRNA positions and target or donor positions.

    Fig. 2: Identification of IS621 bridge RNA binding loops with sequence-specific recognition of target and donor DNA.
    figure 2

    a, Schematic of the computational approach to assess the base-pairing potential between the IS110 ncRNA and its cognate genomic target site or donor sequence. Covariation analysis between target–ncRNA or donor–ncRNA pairs yields a matrix in which diagonal stretches of red signal indicate ncRNA complementarity to the bottom strand of the DNA and blue stretches indicate complementarity to the top strand. b, Nucleotide covariation and base-pairing potential between the ncRNA and the target (left) and donor (right) sequences across 5,511 ncRNA–target pairs and 2,201 ncRNA–donor pairs. The IS621 ncRNA sequence is shown across the x?axis, along with dot-bracket notation predictions of the secondary structure. Covariation scores are coloured according to strand complementarity, with ?1 (blue) representing high covariation and a bias toward top-strand base-pairing, and 1 (red) representing high covariation and a bias toward bottom-strand base-pairing. Regions of notable covariation signal indicating base-pairing for IS621 are boxed. Complementary nucleotides within covarying regions are highlighted in bold. c, Schematic of the in vitro recombination (IVR) reaction with IS621. d,e, Gel electrophoresis of the IVR LD–RT PCR product (d) or LT–RD PCR product (e). Results are representative of three technical replicates. Rec, recombinase. f, Binding of target and donor DNA sequences by an IS621 RNP containing fluorescently labelled recombinase and ncRNA, using MST. Mean?±?s.d. of three technical replicates. g, Schematic of the IS621 bridge RNA. The target-binding loop contains the LTG and RTG (blue), and the donor-binding loop contains the LDG and RDG (orange). h, Base-pairing model of the IS621 bridge RNA with cognate target and donor DNA.

    Full size image

    This combined analysis clearly indicated potential base-pairing between the two internal loops of the ncRNA and the target and donor DNA sequences, respectively (Fig. 2b and Extended Data Fig. 4a,b). Projecting this covariation pattern onto the canonical IS621 sequence and ncRNA secondary structure, we inferred that the first internal loop might base-pair with the target DNA, whereas the second internal loop might base-pair with the donor DNA. The 5′ side of each loop seems to base-pair with the bottom strand of the target or donor with a stretch of eight or nine nucleotides, whereas the 3′ side of each loop seems to base-pair with the top strand of the target or donor using four to seven nucleotides (Fig. 2b). The strong covariation and base-pairing signal suggest that ncRNA base-pairing with target and donor DNA is a conserved mechanism across diverse IS110 orthologues.

    IS621 ncRNA bridges target and donor DNA

    Previous attempts to study IS110 activity have been successful only in IS110 host organisms, with no reports of successful in vitro reconstitution9,12,15. We reasoned that the ncRNA could be the missing component required for recombination. To test this, we combined in-vitro-transcribed ncRNA with purified IS621 recombinase and dsDNA oligonucleotides containing target and donor DNA sequences to assess in vitro recombination. Strikingly, we found that the ncRNA is necessary for in vitro recombination, and that the four components (ncRNA, recombinase, target DNA and donor DNA) are sufficient to produce the expected recombination product (Fig. 2c–e and Supplementary Fig. 1). MST also revealed that the recombinase–ncRNA ribonucleoprotein (RNP) complex binds to wild-type (WT) target and donor dsDNA oligos (target KD?=?13?±?6?nM; donor KD?=?77?±?3?nM), but not to non-complementary DNA molecules (Fig. 2f). Together, these findings indicate that the ncRNA bound by the IS621 recombinase enables sequence-specific binding to both target and donor DNA molecules to facilitate recombination.

    We named this ncRNA ‘bridge RNA’, on the basis of its bispecific role in bridging the target and donor DNA molecules for recombination. We refer to the two internal loops of the bridge RNA as the target-binding loop and the donor-binding loop (Fig. 2g). The target-binding loop comprises two key regions that base-pair with the top and bottom strands of the target DNA, respectively: the left target guide (LTG) base-pairs with the left side of the bottom strand of the target DNA (left target; LT), whereas the right target guide (RTG) base-pairs with the right top strand of the target DNA (right target; RT). The donor-binding loop has an analogous architecture, in which a left donor guide (LDG) base-pairs with the bottom strand of the left donor (LD) and a right donor guide (RDG) base-pairs with the top strand of the right donor (RD) (Fig. 2h). Of note, the core dinucleotide is included in each of the base-pairing interactions (LTG–LT, RTG–RT, LDG–LD and RDG–RD), which results in an overlap between the right top and left bottom strand pairings.

    To lend further support to our hypothesis that the bridge RNA target-binding loop guides the selection of the genomic target sequence, we analysed insertion loci across diverse IS110 orthologues. Binning natural IS110s by sequence similarity of their LTG and RTG, we created a consensus genomic target site motif for each LTG, RTG pair. The target motif was highly concordant with the target-binding loop sequences of the bridge RNA (LTG and RTG), with zero to two mismatches in most cases (Fig. 2b, Extended Data Fig. 5a and Supplementary Data 2). Our covariation data further indicated that the RTG of some IS110 orthologues is longer than the RTG for IS621 (Fig. 2b and Extended Data Fig. 5b,c). We also observed evidence of a distinct base-pairing pattern between the RDG and the RD, in which a stretch of nine bridge RNA nucleotides base-pairs discontiguously with a stretch of seven donor DNA bases (Fig. 2b,h).

    Programmable target site selection

    The base-pairing mechanism of target and donor recognition by the bridge RNA suggests programmability. To assess this, we set up a two-plasmid recombination reporter system in E. coli: pTarget encodes the IS621 recombinase, a 50-bp target site and a promoter, and pDonor encodes the RE–LE donor sequence containing the bridge RNA and a promoter-less gfp. Recombination of pDonor into pTarget would place gfp downstream of the promoter, with successful recombination events detected using flow cytometry (Fig. 3a). Using the WT IS621 donor and target sequences, we detected the expression of GFP and confirmed the expected recombination product using nanopore sequencing (Fig. 3b). Substituting conserved catalytic residues with alanine (Extended Data Fig. 1a,b) abolished recombination, as did substituting pDonor with a version lacking the RE–LE (and therefore lacking the bridge RNA) (Fig. 3b).

    Fig. 3: The IS621 target site is reprogrammable and is specified by the bridge RNA.
    figure 3

    a, Schematic representation of the plasmid recombination assay with bridge RNA in cisb, GFP fluorescence of E. coli after DNA recombination of the plasmid reporter system using catalytic variants of the IS621 recombinase. Plots are representative of three replicates. c, Schematic of reprogrammed target and bridge RNA target-binding loop sequences. d, GFP mean fluorescence intensity (MFI) of E. coli after plasmid recombination using the indicated reprogrammed bridge RNA target-binding loop and target sequences (WT and T1–T7). Bold bases highlight differences relative to the WT target sequence. Mean?±?s.d. of three biological replicates. e, Schematic of bridge RNA expression in transf, Comparison of recombination efficiency with bridge RNA expressed in cis and in trans. Mean?±?s.d. of three biological replicates.

    Full size image

    We next selected seven target sequences (T1–T7) and designed reprogrammed bridge RNAs with matching target-binding loops (Fig. 3c). These T1–T7 reprogrammed bridge RNAs abrogated recombination with the WT target while enabling high rates of recombination (13.8–59.5% of all cells) with each cognate target sequence (Fig. 3d, Extended Data Fig. 6a and Supplementary Fig. 2). We next asked whether the bridge RNA could be expressed in trans rather than within the RE–LE context. We truncated the RE–LE (298?bp) to a 22-bp donor around the core dinucleotide, which eliminated the ?35 box of the natural σ70 promoter (Fig. 3e). This variant of pDonor did not support recombination into the T5 target plasmid (Fig. 3f) until we supplied the full-length T5 bridge-RNA-encoding sequence in a distinct site on pDonor under the control of a synthetic promoter. The in trans bridge RNA increased the total GFP fluorescence signal by nearly twofold compared with the same bridge RNA expressed from the native RE–LE promoter (Fig. 3e,f). Together, these results indicate that the bridge RNA target-binding loop can be reprogrammed to direct target site specificity for DNA recombination in E. coli.

    To comprehensively assess the mismatch tolerance and reprogramming rules of bridge RNAs, we designed an E. coli selection screen that links thousands of barcoded pairs of DNA targets and bridge RNAs on a single plasmid. Successful recombination with a WT donor plasmid induces a kanamycin resistance cassette (KanR) for survival (Fig. 4a). Using this approach, we first confirmed that base-pairing between the bridge RNA and both strands of the CT target core sequence was strongly preferred, in line with the high conservation of the CT core sequence in both the target and the donor (Fig. 4b and Extended Data Figs. 5d,e and 6b).

    Fig. 4: High-throughput characterization of IS621 target specificity shows flexible programmability.
    figure 4

    a, Schematic representation of the target specificity screen. Successful recombination enables the survival of E. coli through the expression of a kanamycin resistance cassette (KanR). The target sequence and bridge RNA are separated by a 12-nt barcode (BC). NGS, next-generation sequencing. b, Mismatch tolerance of the core dinucleotide. Core-binding nucleotides of the target-binding loop are summarized by IUPAC codes, including D (not C) and V (not U). Average counts per million (CPM) of two biological replicates. Box plots show median (centre line), IQR (box edges) and 1.5 × IQR (whiskers). c, Mismatch tolerance between non-core sequences of the target and target-binding loop. Average CPM of two biological replicates. Box plots show median (centre line), IQR (box edges) and 1.5 × IQR (whiskers). d, Mismatch tolerance between target and target-binding loop, as indicated by the percentage of total detected recombinants carrying each nucleotide at each position. Average of two biological replicates. e, Nucleotide enrichment among the top 20% most efficient matched pairs of targets and target-binding loops. f, Schematic of the genome insertion assay in E. colig, Genome-wide mapping of insertions mediated by the WT IS621 bridge RNA. The percentage of total reads mapped to each insertion site is depicted and binned by the number of differences from the intended sites as measured by Levenshtein distance. Average of two biological replicates. h, Target site preference of IS621. Sequence logos depict the target site motifs among natural (top, Methods) and experimentally observed (bottom, Fig. 4g) IS621 target sites. i, Genomic specificity profile of four reprogrammed bridge RNAs. Two biological replicates.

    Full size image

    Next, we varied the nine non-core positions of the target and the corresponding positions of the LTG and RTG to assess single and double mismatch tolerance at each position. We observed a strong preference for perfect matches across all nine positions of the target-binding loop and target, and a high degree of reprogramming flexibility at all positions (Fig. 4b–d, Extended Data Fig. 6b,c, Supplementary Table 3 and Supplementary Fig. 3). As expected, double mismatches were even less tolerated than were single mismatches, with bias for certain combinations of mismatch positions (Fig. 4c and Extended Data Fig. 6d). Overall, we show that the target-binding loop is broadly programmable at each position, with a low mismatch tolerance (Fig. 4e).

    Programmable insertion in the E. coli genome

    To evaluate the genomic site selection and specificity of WT IS621, we measured the insertion of a replication-incompetent plasmid (4.85?kb) bearing the 22-bp WT donor sequence into the E. coli genome using the WT IS621 bridge RNA and recombinase (Fig. 4f). After selection, we mapped insertions genome-wide and observed 173 unique insertion sites, with 144 of these insertions occurring within the REP elements that are known16 to be targeted by IS621 (Fig. 4g, Supplementary Table 3 and Supplementary Fig. 4). Of all insertion sites, 74.5% (129 sites) matched the naturally observed target sequence (ATCAGGCCTAC), and two more sites exactly matched the specificity encoded by the target-binding loop (ATCGGGCCTAC); together, these accounted for 96.21% of all detected insertions (Extended Data Fig. 7a–c). Our assay therefore recapitulated the specificity of IS621 elements found in nature, including tolerance for a mismatch at position 4 of the target site (Fig. 4h). Structural analysis of the IS621 recombination complex indicates that this mismatch results in a non-canonical rG:dT base pair, which could explain the high frequency of insertions into these target sites22.

    Further scrutiny of the insertion sites revealed that four of the ten most frequently targeted sites were flanked on the 3′ end of the RT sequence by 5′-GCA-3′—complementary to the 5′-UGC-3′ that occurs immediately 5′ of the RTG in the WT bridge RNA (Fig. 2b and Extended Data Fig. 7a–d). This suggested to us the potential of an extended base-pairing interaction beyond the predicted RTG–RT for IS621 (7?bp instead of 4?bp), which was supported by the observation that some IS110 orthologues naturally encode longer RTGs (Fig. 2b and Extended Data Fig. 5b,c).

    To investigate genome-wide insertion specificity, we reprogrammed bridge RNAs to target sequences found only once in the E. coli genome. We tested four distinct genomic sites with two bridge RNAs for each: one containing a short 4-bp RTG (IS621 RTG) and one with a long 7-bp RTG (Extended RTG) to directly assess the effect of RTG–RT base-pairing length on specificity. In each case, we found that the expected genomic target site was the most frequently targeted, representing between 51.6% and 94.0% of all detected insertions (Fig. 4i). Off-target insertions were also observed, with individual off-target sites each representing between 0.11% and 31.16% of insertions across all bridge RNAs, with the more frequently detected off-targets typically carrying one or two mismatches with the expected target (Extended Data Fig. 7e).

    The extended RTG improved the specificity of insertion into the on-target site from an average of 69.4% (range 51.2–89.4%) to an average of 84.9% (range 65.4–94.0%). It also resulted in markedly fewer insertions into off-target sites for bridge RNA 2 and bridge RNA 3, eliminating 18 out of 45 and 14 out of 25 off-target sites, respectively (Fig. 4i). Notably, some off-target sites seemed to indicate tolerance for insertions in the target sequence, whereas some low-frequency insertions seemed to more closely resemble the 11-bp WT donor sequence, rather than the programmed target (Extended Data Fig. 7e,f). Of the 117 genomic off-target insertion sites detected across the 8 experiments, 102 (87.2%) had the expected CT dinucleotide core, 56 (47.9%) closely resembled the target or donor sequence (Levenshtein distance?<?3) and the remaining sites were enriched for long k-mer matches to the target or donor sequence (Extended Data Fig. 7g), suggesting that most or all of the detected off-target insertions were bridge-RNA-dependent. In addition to off-target insertions, genomic deletions and inversions between experimentally observed insertion sites were detected in rare cases (allele frequency?<?0.05) (Supplementary Note 1). Altogether, these experiments provide evidence of the robust capability of IS621 to specifically insert multi-kilobase cargos into the genome, and offer further insights into the mechanism of recombination.

    Programming the donor specificity of bridge RNAs

    Among IS621 elements, the donor sequence is more highly conserved than the genomic target sequence, which suggests that the donor-binding loop may be less readily reprogrammed than the target-binding loop (Extended Data Fig. 5d,e). To assess this, we designed a donor specificity screen in which we varied the 7-bp LD and 2?bp of the RD flanking the core dinucleotide, all within the context of a full-length RE–LE expressing the bridge RNA in cis. Successful recombination with the T5 target plasmid would induce KanR expression (Fig. 5a). Analysis of thousands of donor and donor-binding loop pairs revealed that the donor sequence can be fully reprogrammed (Fig. 5b). Similar to the interaction between the target and the target-binding loop of the bridge RNA, LD–LDG mismatches and RD–RDG mismatches were generally poorly tolerated (Fig. 5c, Supplementary Fig. 5 and Supplementary Table 4). Position 7 of the LD was an exception, exhibiting a strong bias against cytosine and therefore appearing to be more mismatch tolerant than other positions (Fig. 5d,e).

    Fig. 5: Bridge RNA donor recoding enables fully programmable insertion, inversion and excision.
    figure 5

    a, Schematic representation of the donor specificity screen. A unique molecular identifier (UMI) identifies each paired donor and donor-binding loop. b, Reprogrammability of donor sequences by the number of nucleotide differences from the WT donor. WT donor abundance is indicated by the dashed line. Average CPM of two biological replicates. Box plots show median (centre line), IQR (box edges) and 1.5 × IQR (whiskers). c, Mismatch tolerance between non-core sequences of the donor-binding loop and donor. Average CPM of two biological replicates. Box plots show median (centre line), IQR (box edges) and 1.5 × IQR (whiskers). d, Mismatch tolerance between bridge RNA donor-binding loop and donor by position, as measured by the percentage of total detected recombinants with each indicated mismatch. Average of two biological replicates. e, Nucleotide enrichment among the top 20% most efficient matched pairs of donors and donor-binding loops. f, Schematic representation of the paired reprogramming of the donor and the donor-binding loop. g, Specific recombination using reprogrammed donor and donor-binding loop sequences. Donor sequences are listed on the left, and the bridge RNA is reprogrammed to base-pair with the indicated sequence. Bold bases highlight differences relative to the WT donor sequence. Mean?±?s.d. of three biological replicates. h, Schematic representation of the programmable excision assay. i, Schematic representation of the programmable inversion assay. j, Efficient programmable excision of DNA. Pairs of donor and target are denoted. k, Efficient programmable inversion of DNA. Pairs of donor and target are denoted. In j,k, negative control (NC) expresses the reporter and recombinase but no bridge RNA; and data are MFI?±?s.d. of three biological replicates.

    Full size image

    In these experiments, the core dinucleotide (CT) was held constant, which could limit the sequence space of potential target and donor sites. To address this, we modified the cores of target T5 and the WT donor, along with their associated bridge RNA positions in both loops, from CT to AT, GT or TT (Extended Data Fig. 8a,b). Although non-CT cores were generally less efficient, efficiency was improved by extending the length of RTG–RT base-pairing from 4?bp to 7?bp, informed by our previous results on RTG extension (Fig. 4i and Extended Data Fig. 8c,d).

    Next, we investigated the ability of the bridge RNA to combinatorially control the recognition of target and donor sequences simultaneously. Using our in trans GFP reporter assay, in which the target-binding loop of the bridge RNA recognizes target T5 (Fig. 3e), we reprogrammed the donor sequence and the donor-binding loop of the bridge RNA to one of nine distinct donor sequences (D1–D9) with varying levels of divergence from the WT donor (Fig. 5f). D1–D9 reprogrammed donor-binding loops supported robust recombination with their cognate donor sequences (26.9–95.0% of all cells) but not with the WT donor (Fig. 5g and Extended Data Fig. 8e). Together, these data show that the bridge RNA allows modular reprogramming of both target and donor DNA recognition.

    Programmable DNA rearrangements

    In addition to their use for DNA insertion, recombinases such as Cre have been routinely used for the excision or inversion of DNA sequences. Typically, such approaches require pre-installation of the loxP recognition sites in the appropriate arrangement, with two sites oriented in the same direction resulting in excision, and sites oriented in opposite directions resulting in inversion. Given our understanding of the IS621 insertion mechanism, as well as the reported existence of invertase homologues of IS110s14,23, we hypothesized that IS621 recombinases could mediate programmable excision and inversion.

    We first generated GFP reporter systems for both excision and inversion (Fig. 5h,i and Extended Data Fig. 9a–c). Testing the same four pairs of donor and target recognition sites in both reporters, we showed that both excision and inversion occur robustly and in a programmable manner (32.2–98.9% and 4.54–98.2% of all cells, respectively) (Fig. 5j,k and Extended Data Fig. 9d,e). Overall, the ability of IS110 recombinases and their bridge RNAs to insert, excise and invert DNA in a programmable and site-specific manner enables remarkable control over multiple types of DNA rearrangements with a single unified system.

    Diverse IS110s encode bridge RNAs

    Finally, we investigated whether the bridge RNA is a general feature of the IS110 family. The IS110 family is divided into two groups: IS110 (which includes IS621) and IS1111. IS1111 elements also encode DEDD recombinases, but have been categorized into a separate group on the basis of the presence of sub-terminal inverted repeat sequences (STIRs) that range in length from 7 to 17 bp8,10,13. We examined our covariation analysis of IS110 group termini and identified a short 2–3-bp STIR pattern that flanks the programmable donor sequence, suggesting an evolutionary relationship with the longer STIRs of IS1111 elements (Extended Data Fig. 10a,b). Amongst all IS110 and IS1111 elements annotated in the ISfinder database, we found that IS1111 elements have much longer REs than LEs—in contrast to the IS110 subgroup, in which the LE is significantly longer than the RE (Fig. 6a).

    Fig. 6: IS110 subfamilies encode distinct and diverse bridge RNA secondary structures in different non-coding end sequences.
    figure 6

    a, Non-coding end length distribution for IS110 and IS1111 group elements. Box plots show median (centre line), IQR (box edges) and 1.5 × IQR (whiskers). b, Location of predicted bridge RNA for IS110 and IS1111 group elements. c, Phylogenetic tree of the 274 IS110 recombinases catalogued by ISfinder. d, Bridge RNA consensus structures from six diverse IS110 elements. Secondary structures are shown with internal loops coloured according to the sequence that they complement: target (blue), donor (orange) or core (green).

    Full size image

    Using RNA structural covariance models, we predicted a bridge RNA in 85.7% of IS110s and 93.0% of IS1111s (Fig. 6b). The vast majority of IS110 group members appeared to encode a bridge RNA within the LE, whereas IS1111 group members appeared to encode a bridge RNA within the RE. This is consistent with a previous report that correlated target site preference with sequence conservation in the RE of IS1111 elements and, on this basis, speculated that an RNA might be involved in target site selection24. Notably, the location of the bridge RNA closely predicted the phylogenetic relationship between IS110 and IS1111, which strongly suggests that these two groups emerged from a common ancestor in which the bridge RNA translocated between the ends of the element and the length of the STIR was modified (Fig. 6c and Supplementary Table 5).

    We predicted bridge RNA structures and manually inspected the loops of six diverse IS110 and IS1111 elements for evidence of complementarity with their cognate target and donor sequences. This analysis yielded diverse structures with clear evidence of a base-pairing pattern (8–14?nt) between internal bridge RNA loops and DNA targets and donors (Fig. 6d and Extended Data Fig. 10c). Of note, in many IS1111 orthologues, the predicted bridge RNA has potential donor-binding nucleotides in a multi-loop structure rather than the simple internal loop observed for IS621 and other members of the IS110 group. Altogether, we conclude that the IS110 family encodes diverse predicted bridge RNAs that direct sequence-specific and programmable recombination between target and donor sequences.

    Discussion

    Non-coding RNA molecules that specify a nucleic acid target are central to both prokaryotic and eukaryotic life. Nucleic acid guides are a widely used mechanism in fundamental biological processes; for example, the tRNA anticodons that govern ribosomal translation; small interfering RNAs and microRNAs of RNA interference; CRISPR RNAs of CRISPR–Cas immunity; and small nucleolar RNAs (snoRNAs) for gene regulation. The bridge RNA that we discovered in this work is the first example, to our knowledge, of a bispecific guide molecule that encodes modular regions of specificity for both the target and the donor DNA, coordinating these two DNA sequences in close proximity to catalyse efficient recombination. Bridge RNAs encode all of this complex molecular logic in a remarkably compact (around 150–250-nt) sequence along with their single effector recombinase (around 300–460-aa) partner.

    IS110 targeting is achieved using internal binding loops that are reminiscent of tRNA hairpin loops or snoRNA internal loops, distinct from the terminal binding sequences of CRISPR–Cas or Argonaute guide RNAs. Each RNA loop encodes segments that base-pair with staggered regions of the top and bottom strand of each cognate DNA binding partner, in contrast to the single-strand base-pairing mechanisms of known RNA-guided systems. Furthermore, the RNA-guided self-recognition of the IS110 element in donor form illustrates a previously unobserved mechanism of DNA mobility.

    Mobile genetic elements have been shaped throughout evolution to insert, excise, invert, duplicate and otherwise rearrange DNA molecules. Bridge RNAs enable IS110 recombinases to exploit the inherent logic of RNA–DNA base-pairing, directly bypassing the complex target site recognition codes of other known transposases and recombinases, which depend on extensive protein–DNA or short single-stranded DNA–DNA interactions that offer much less opportunity for straightforward programmability25,26,27,28. The IS110 family is evolutionarily diverse and widespread across prokaryotes, providing a rich landscape for further functional insights. In our initial survey of diverse IS110 orthologues, we uncovered a variety of bridge RNA structures and lengths, suggesting that there is considerable mechanistic diversity both between and within each of the two major IS110 and IS1111 subfamilies.

    Our accompanying cryo-electron microscopy analysis of the IS621 recombinase in complex with bridge RNA, target DNA and donor DNA, captured in several stages of the recombination reaction, is copublished with this paper22. Together, our two studies detail the unique mode of dual-strand recognition of the target and donor DNA through programmable base-pairing interactions with the bridge RNA. The synaptic complex structures illustrate how two recombinase dimers associate with the target-binding loop and the donor-binding loop of bridge RNAs, coming together to form an adaptable recombination complex (ARC) with composite subunit-spanning active sites when both target and donor DNA are engaged by the ARC system. This elegant licensing mechanism enables nicking and exchange of the top strands between the donor and target, resulting in a Holliday junction intermediate that is resolved by the cleavage of the bottom strands. Together, our genetic, mechanistic, computational and structural characterization of the bridge recombination system lays the foundation for protein and RNA engineering efforts to improve and optimize its capabilities.

    Guide RNAs are underpinning a technological revolution in programmable biology29,30,31,32,33,34,35. The direct enzymatic activity of stand-alone, naturally occurring programmable RNA-guided proteins has been notably limited to the endonuclease function30,36. Successive generations of programmable nucleases and nickases have advanced the prevailing genome-editing method from the original homology-based capture of a DNA donor37 to the targeted stimulation of donor insertion, all of which require a complex interplay with endogenous DNA repair processes31,34,38,39,40. Functional diversification of these systems beyond nucleic acid binding or cleavage has generally required the recruitment or fusion of additional effector proteins, resulting in increasingly large and intricate engineered genome-editing fusions41,42. The IS110 bridge system, in contrast, uses a single and compact RNA-guided recombinase that is necessary and sufficient for direct DNA recombination (Fig. 2d,e). Modular reprogramming of target and donor recognition by the bispecific bridge RNA uniquely enables the three fundamental DNA rearrangements of insertion, excision and inversion for manipulating large-scale DNA sequences and overall genome organization. With further exploration and development, we expect that the bridge recombination mechanism will spur a third generation of programmable RNA-guided tools beyond RNA interference- and CRISPR-based mechanisms to enable a new frontier of genome design.

    Methods

    Development of metagenomic and genomic sequence database

    A custom sequence database of bacterial isolate and metagenomic sequences was constructed by aggregating publicly available sequence databases, including NCBI, UHGG43, JGI IMG44, the Gut Phage Database45, the Human Gastrointestinal Bacteria Genome Collection46, MGnify47, Youngblut et al. animal gut metagenomes48, MG-RAST49 and Tara Oceans samples50. The final sequence database included 37,067 metagenomes, 274,880 bacterial and archaeal metagenome-assembled genomes, 855,228 bacterial and archaeal isolate genomes and 185,140 predicted viral genomes.

    Analysis of conserved residues in IS110 protein sequences

    Genomic sequences were annotated using Prodigal51 to identify coding sequences. All unique protein sequences were then combined into a single FASTA file and clustered at 30% sequence identity using mmseqs2 (ref. 52). Two Pfam domains DEDD_Tnp_IS110 (PF01548) and Transposase_20 (PF02371) were used to search against these clustered representative proteins using the hmmsearch tool in the hmmer package53. DEDD_Tnp_IS110 was used to identify the RuvC-like domain, and Transposase_20 was used to identify the Tnp domain. All members of the matched 30% identity clusters were then extracted, and the same IS110 Pfam domain significance thresholds were applied to filter these candidates. Next, only proteins that met E?<?1?×?10?3 for both domains were retained. Next, RuvC-like domains were only retained if they were between 125 and 175 aa in length, and Tnp domains were only retained if they were between 60 and 110 aa in length. Any sequences with ambiguous residues were removed. Protein domains were then clustered at 90% using mmseqs (‘easy-cluster --cluster-reassign -c 0.8 --min-seq-id 0.9 --cov-mode 0’). Cluster representatives were then aligned using hmmalign (‘--trim --amino’)53. Alignment columns with more than 50% gaps were removed, and the alignments were visualized using ggseqlogo in R54.

    Phylogenetic analysis of IS110 transposases

    A phylogenetic analysis of IS110 transposases was also performed. Full-length IS110 proteins were clustered at 90% identity using the mmseqs2 easy-cluster algorithm (‘--cluster-reassign -c 0.85 --min-seq-id 0.9 --cov-mode 0’)52. Next, using the identified 90% protein sequence clusters, a representative from each cluster was selected that was closest to the 80th percentile in total length. This resulted in a curated set of 90% identity cluster representatives. Next, 90% identity cluster representatives were clustered at 30% identity across 70% of the aligned sequences using the mmseqs2 easy-cluster algorithm (‘--cluster-reassign -c 0.70 --min-seq-id 0.30 --threads 96 --cov-mode 0’). This resulted in 1,686 30% identity cluster representatives. RuvC-like and Tnp-like domains were extracted from these proteins using the corresponding Pfam pHMM models and hmmsearch53. These extracted domains were then individually aligned using hmmalign (‘--amino --trim’) and concatenated into a paired alignment. All pairwise percentage identity values were calculated for this alignment, and redundant sequences were removed using a 60% identity cut-off, resulting in 1,054 aligned sequences. A phylogenetic tree was then constructed using iqtree2 v.2.1.4-beta, with all default parameters55, midpoint rooted and visualized in R with ggtree56. Additional metadata about each sequence was mapped onto the tree, including host kingdom and phylum, ISfinder group and notable orthologues.

    Curated ISfinder transposases were analysed separately to produce another phylogenetic tree. IS110 transposase sequences were extracted from the database available through the prokka software package57. Only IS110 transposases of more than 250 aa were retained. Protein sequences were then clustered using mmseqs2 (‘easy-cluster -c 0.5 --min-seq-id 0.9 --threads 8 --cov-mode 0’)52. Cluster representatives were then aligned using mafft-ginsi (‘--maxiterate 1000’)58. Alignment columns with more than 50% gaps were removed. A phylogenetic tree was then constructed using iqtree2 v.2.1.4-beta with all default parameters55.

    Analysis of LE and RE lengths across IS110 elements

    Sequence coordinate information about individual IS elements was collected through the ISfinder web portal59. This included information about the total length of each IS element, as well as the start and end coordinates of the recombinase CDS. The LE non-coding length was calculated from the CDS coordinates for each IS110 element as the distance between the 5′ terminus and the start of the CDS, and the RE non-coding length was calculated as the distance between the end of the CDS and the 3′ terminus. Tn3 family elements were excluded owing to highly variable passenger gene content.

    Predicting IS110 element boundaries

    To identify the boundaries of each element, an initial search was conducted using comparative genomics to identify putative pre-insertion and post-insertion examples within the metagenomic sequence database. IS110 protein candidates were clustered at 30% identity using mmseqs2 (ref. 52), and within each cluster all relevant genomic loci were identified. Nucleotide sequences were then extracted from the database by adding 1,000 base pairs to the 5′ and 3′ ends of the IS110 CDS, and extracting the complete intervening sequence. These IS110 loci were then separated into ‘batches’ on the basis of 90% identity protein clusters. These batches were then searched against up to 40 metagenomic or isolate samples in the custom database, prioritizing samples that already contained related recombinases. Putative pre-insertion sites were identified if the distal ends of the loci aligned by BLAST to a contiguous sequence60, but the IS110 CDS did not. Precise boundaries of the IS110 element were then predicted using a modified method similar to what was implemented by the previously published tool MGEfinder61. Core sequences were identified as repeated sequences near the end of the predicted element. Next, an iterative BLAST search was used to extend IS110 element boundary predictions beyond those that could be detected by identifying pre-insertion sites. IS110 elements were searched using BLAST against all IS110 loci. Hits were retained only if both ends of the element aligned, and if the core was concordant between query and target. This then generated a new set of IS110 elements and their boundaries, which were recycled as query sequences, and the search was repeated for another iteration. This repeated for 36 iterations before convergence (no new IS110 elements were found). The combined set of IS110 boundaries were kept for further analysis.

    Identification of bridge RNA consensus structures

    A pipeline was developed to identify conserved RNA structures in the sequences immediately flanking the recombinase CDS. First, the IS621 protein sequence was searched against the complete IS110 database for orthologues using blastp (‘-max_target_seqs 1000000 -evalue 1e-6’). Only hits that were at least 30% identical at the amino acid level with 80% of both sequences covered by the alignment were retained. Up to 2,000 unique proteins were then selected in order of descending percentage amino acid identity. Flanking sequences for the corresponding proteins were then retrieved from the database, with flanking sequences defined as a 5′ flank of up to 255?bp (including 50?bp of 5′ CDS) and a 3′ flank of up to 170?bp (including 50?bp of the 3′ CDS). These flanks were then further filtered to exclude sequences that were more than 35 bases shorter than the target flank lengths. Sequences were filtered to exclude those with ambiguous nucleotides. Protein sequences were then clustered using mmseqs2 easy-linclust with a minimum percentage nucleotide identity cut-off of 95% across 80% of the aligned sequences, and one set of flanks for each representative was retained. Flanking sequences were then clustered at 90% nucleotide identity across 80% of the aligned sequences, and only one representative flanking sequence pair per cluster was retained. Then, up to 200 sequences were selected in order of decreasing percentage identity shared between the IS621 protein sequence and their corresponding orthologue protein sequence. The remaining sequences were then individually analysed for secondary RNA structures using linearfold62. Sequences were then aligned to each other using the mafft-xinsi (IS621 orthologue sequences) or mafft-qinsi (all other ISfinder elements) alignment algorithms and parameter --maxiterate 1000 (ref. 58). Alignment columns with more than 50% gaps were removed. The conserved RNA secondary structure was then projected onto the alignment, and manually inspected to nominate bridge RNA boundaries. This region was exported as a separate sequence alignment file, and a consensus RNA secondary structure was predicted using ConsAlifold63. This structure was then visualized using R2R64. This same pipeline was used to analyse hundreds of other IS110 elements, resulting in diverse predicted secondary structures. For visualization purposes, consensus secondary structures with minimally structured terminal ends were trimmed to the primary structured sequence. These consensus structures were converted into covariance models using Infernal20, and these were then searched across thousands of sequences to identify putative bridge RNAs20.

    Nucleotide covariation analysis to identify bridge RNA guide sequences

    To identify programmable guide sequences in the bridge RNA of the IS621 element, the following approach was taken. First, the IS621 protein sequence was searched against our collection of IS110 recombinase proteins with predicted element boundaries using blastp. Next, only alignments that met a cut-off of 20% amino acid identity across 90% of both sequences were retained. Next, a covariance model of the bridge RNA secondary and primary sequence was used to identify homologues of the bridge RNA sequence in the non-coding ends of these orthologous sequences20. Fifty nucleotide target and donor sequences were extracted centred around the core. For elements with multiple predicted boundaries, boundaries with a CT dinucleotide core were prioritized. Next, elements that were identified at earlier iterations in our boundary search were prioritized. Next, elements that were similar in length to the known IS621 sequence element were prioritized. Only one element per unique locus was retained. Alignments were further filtered to remove redundant examples by clustering targets or donors and bridge RNA sequences at 95% identity, taking one representative per pair and then taking at most 20 examples for each 95% identity bridge RNA cluster. Predicted bridge RNA sequences were then aligned using the cmalign tool in the Infernal package20. Two paired alignments were then generated that contained concatenated target and bridge RNA sequences, and concatenated donor and bridge RNA sequences. These alignments were then further filtered to remove all columns that contained gaps in the IS621 bridge RNA sequence. These alignments were then analysed using CCMpred (‘-n 100’) to identify covarying nucleotides between targets or donors and bridge RNA sequences65. These covariation scores were normalized by min-max normalization and multiplied by the sign of the column-permuted base-pairing concordance score (see next paragraph), with +1 corresponding to bottom-strand base-pairing and ?1 corresponding to top-strand base-pairing. The signal was visualized as a heat map and interactions were identified within the two internal loops of the bridge RNA, leading to the proposed model for bridge RNA target or donor recognition. The same covariation analysis was performed on the donor alone, leading to the identification of short STIR sequences for IS110 elements.

    A separate analysis was performed on the same paired alignment used in the covariation analysis to determine whether certain pairs of nucleotides were biased toward base-pairing. The observed concordance was first calculated for each pair of columns as:

    ??????=absmax(??=1??CheckEqual(??????,??????),??=1??CheckComplementary(??????,??????))??,

    where C is the concordance score, i refers to the first column (or position), j refers to the second column, n refers to the total number of rows (sequences) in the alignment, ski refers to the nucleotide in bridge RNA sequence k at position i and tkj refers to the nucleotide in target (or donor) sequence k at position j. absmax(a,b) is a function that returns the value with the largest absolute magnitude, CheckEqual(a,b) is a function that returns one when a?=?b and 0 otherwise and CheckComplementary(a,b) is a function that returns ?1 if a and b are complementary nucleotides and 0 otherwise. All positions in which the nucleotide is a gap in either sequence are ignored and discounted from n. All observed values of Cij are then compared with two different null distributions of Cij scores. The first is generated by randomly permuting the rows of the bridge RNA alignment 1,000 times and recalculating C for each permutation, and the second is generated by randomly permuting the columns of the bridge RNA alignment 1,000 times and recalculating C. The mean and standard deviation of these permuted C distributions are then used to convert the observed C scores into z-scores, and positive and negative values are then separately min-max normalized to maintain the ?1 to 1 scale. The sign of this score is then used to project base-pairing information onto the covariation scores as generated by CCMpred.

    Small RNA-seq of IS110 bridge RNAs

    BL21(DE3) E. coli were transformed with plasmids bearing a concatenated RE–LE sequence and plated on an LB agar plate with appropriate antibiotics. A single colony was picked and grown in terrific broth (TB) to an optical density (OD) of 0.5. RNA isolation was performed using the Direct-zol RNA Miniprep kit (Zymo Research). RNA was prepared for small RNA-seq according to the following protocol. In brief, no more than 5?µg total RNA was treated with DNase I (NEB) for 30?min at 37?°C then purified using the RNA Clean & Concentrator-5 kit. Ribosomal RNA was depleted from samples using the Ribo-Zero Plus rRNA Depletion kit (Illumina) and purified using the RNA Clean & Concentrator-5 kit. Depleted RNA was treated with T4 PNK for six hours at 37?°C, supplementing with T4 PNK and ATP after six hours for one additional hour. RNA was purified using the RNA Clean & Concentrator-5 kit and subsequently treated with RNA 5′ polyphosphatase (Lucigen) for 30?min at 37?°C. RNA was purified with the RNA Clean & Concentrator-5 kit, and the concentration was measured by NanoDrop. NGS libraries were prepared using the NEBNext Multiplex Small RNA Library Prep Kit (NEB) according to the manufacturer’s protocol. Libraries were sequenced on an Illumina MiSeq using a 2×150 Reagent kit (v.2).

    Analysis of small RNA-seq data

    Demultiplexed fastq files were cleaned and merged using BBtools (bbduk and bbmerge), respectively66. Merged fastq files were aligned to the RE–LE-bearing plasmid using bwa-mem67. Small RNA-seq coverage was normalized according to the maximum read depth observed for each orthologue across the entire RE–LE plasmid.

    In vitro transcription of bridge RNAs

    In vitro transcription was performed on a linear DNA template using the HighScribe T7 High Yield RNA Synthesis Kit (New England Biolabs) as per the manufacturer’s instructions. The DNA template was prepared by cloning into a pUC19 backbone and the plasmid was linearized using the SapI restriction enzyme (NEB) and purified using DNA Clean & Concentrate (Zymogen). After in vitro transcription, RNA was purified using the Monarch RNA Cleanup kit. Where necessary, bridge RNA was further purified by denaturing polyacrylamide gel electrophoresis, extracted from the gel using UV shadowing and recovered by ethanol precipitation.

    IS621 protein preparation

    The IS621 recombinase gene was human codon optimized and cloned into a modified pFastBac expression vector (Addgene, 30115), which includes an N-terminal His6-tag, a TwinStrep-tag and a tobacco etch virus (TEV) protease cleavage site. To express IS621 recombinase protein Sf9 cells (ATCC, CRL-1711) were cultured in Sf-900 III SFM medium (Thermo Fisher Scientific) supplemented with 10?µg?μl?1 gentamicin and 5% heat-inactivated fetal bovine serum (Gibco). For baculovirus production, recombinant bacmids were first generated by transforming MAX Efficiency DH10Bac competent cells (Thermo Fisher Scientific) with the pFastBac construct. Site-specific recombination between pFastBac and the baculovirus shuttle vector was then confirmed by PCR and Sanger sequencing. For large-scale protein expression, a high-titre P1 recombinant (pFastBac) baculovirus stock was used; cells were infected with pFastBac baculovirus at a multiplicity of infection of 5 plaque-forming units per cell at a cell density of 3?×?106 cells per ml and grown in suspension culture at 28?°C. Cells were collected by centrifugation (300g, 15?min, 4?°C) 48?h after infection and lysed by sonication in buffer containing 20?mM Tris-HCl, pH 7.5, 1?M NaCl, 2?mM MgCl2, 1?mM dithiothreitol (DTT), 10% glycerol and 2% Triton-X, supplemented with cOmplete EDTA-free mini protease inhibitor cocktail (Roche). Then the lysate was clarified by ultracentrifugation at 45,000g and filtered through a 0.45-µm PVDF syringe filter (Millipore Sigma). The supernatant was applied to a 5-ml Strep-Tactin Superflow high-capacity FPLC column (IBA Lifesciences) and washed with 20 column volumes of wash buffer containing 20?mM Tris-HCl, pH 7.5, 0.5?M NaCl, 2?mM MgCl2, 1?mM DTT and 10% glycerol, and the protein was eluted with wash buffer containing 80?mM biotin. Eluted protein was concentrated using a 10-kDa molecular weight cut-off (MWCO) ultracentrifugal concentrator (Millipore Sigma) at 4?°C and the His-TwinStrep-tag was cleaved using TEV protease (NEB) at 37?°C for 4?h. His-TwinStrep-tag cleaved protein was then applied to a 5?ml HisTrapFF Crude immobilized metal affinity column (Cytiva) equilibrated with wash buffer containing 20?mM Tris-HCl, pH 7.5, 0.5?M NaCl, 2?mM MgCl2, 1?mM dithiothreitol (DTT) and 10% glycerol. Wash fractions expected to contain His-TwinStrep-tag-removed IS621 recombinase protein were collected and bound protein was eluted using wash buffer containing 0.5?M imidazole. Notably, IS621 recombinase remained bound to the HisTrapFF column despite His-TwinStrep-tag removal and eluted in the presence of high imidazole. Finally, elution fractions containing recombinant protein were concentrated using a 10-kDa-MWCO ultracentrifugal concentrator (Millipore Sigma) and buffer exchanged during centrifugation into size-exclusion chromatography (SEC) buffer containing 20?mM Tris-HCl, pH 7.5, 0.5?M NaCl, 2?mM MgCl2, 1?mM DTT and 10% glycerol. SEC was performed using a Superdex 200 Increase 10/300 GL column (Cytiva) to further purify the protein, and the peak fractions were collected, concentrated as described above and stored at ?80?°C until use.

    Microscale thermophoresis (MST)

    MST was performed using a Monolith NT.115Pico series instrument (NanoTemper Technologies). IS621 recombinase was labelled for MST using the RED-MALEIMIDE 2nd Generation cysteine reactive kit (NanoTemper Technologies) as per the manufacturer’s instructions. Labelled protein was eluted in a buffer containing 20?mM Tris-HCl, 500?mM NaCl, 5?mM MgCl2, 1?mM DTT and 0.01% Tween20, pH 7.5. To determine the affinity of recombinase for RNA, 20?nM recombinase was incubated with a dilution series (2,500–0.076?nM) of purified LE-encoded ncRNA or a scrambled RNA of equivalent length. MST was performed at 37?°C using premium capillaries (NanoTemper Technologies) at 30% LED excitation and medium MST power. Data were analysed using the NanoTemper MO.affinity analysis (v.3.0.5) software package and raw data were plotted on GraphPad Prism (v.10.2.0) for visualization. The binding affinities of the IS621 RNP for donor and target DNA, as well as for donor and target DNA containing scrambled LD–RD and LT–RT regions, were determined using the MST tertiary binding function. Single-stranded DNA was purchased from IDT and annealed in buffer containing 10?mM Tris pH 8.0, 5?mM MgCl2 and 5?mM KCl. For MST, 20?nM RNP consisting of labelled IS621 recombinase and LE-encoded ncRNA were incubated with a dilution series of duplexed donor or target DNA oligonucleotides (10?µM to 0.076?nM). MST was performed at 37?°C using premium capillaries (NanoTemper Technologies) at medium MST power with the LED excitation power set to automatic (excitation ranged from 20% to 50%).

    In vitro recombination assay

    The in vitro activity of IS621 recombinase was evaluated by incubating 10?µM IS621 with 20?µM LE-encoded ncRNA and 0.5?µM duplexed target and donor DNA oligonucleotides (Supplementary Information) in buffer containing 20?mM Tris-HCl, 300?mM NaCl, 5?mM MgCl2, 1?mM DTT, 0.05?U?µl?1 SUPERase•In RNase Inhibitor (Invitrogen) at 37?°C for two hours. Reactions were then treated with 40?µg Monarch RNaseA (NEB) for one hour and then treated with 1.6?units of Proteinase K (NEB) for a further hour before clean-up of DNA with AMPure XP Beads (Beckman Coulter) using a 2× bead ratio. To detect recombination products, 0.5?µl of the purified reaction product was PCR-amplified with primers designed to amplify the LT–RD and LD–RT recombination products. PCR products were visualized by running PCR reactions on 8% TBE gel (Invitrogen) and staining with SYBR Safe (Thermo Fisher Scientific), and were imaged on a ChemiDoc XRS+ (Bio-Rad). PCR products were sequenced using Oxford Nanopore sequencing (Primordium Labs).

    Plasmid recombination assay in E. coli

    BL21(DE3) cells (NEB) were co-transformed with a pTarget plasmid encoding a target sequence and a T7-inducible IS621 recombinase and a pDonor plasmid encoding a bridge RNA, a donor sequence and a GFP CDS upstream such that after recombination into pRecombinant GFP, expression would be activated by the synthetic Bba_R0040 promoter adjacent to the target site. When expressing the bridge RNA in cis, pDonor encodes a full-length RE–LE sequence (298?bp), which naturally encodes the donor, the bridge RNA and a promoter to express the bridge RNA. When expressing the bridge RNA in trans, pDonor encodes a shortened donor sequence (22?bp) and a bridge RNA driven by the J23119 promoter and followed by the HDV ribozyme.

    To measure excision, a Bba_R0040 promoter is separated from the GFP CDS by the donor site, 1?kb of intervening DNA sequence including an ECK120029600 to terminate transcription, and a target site on the same strand. Co-expression of a second plasmid encoding a bridge RNA and a T7-inducible IS621 recombinase results in the excision of the intervening 1-kb sequence, yielding GFP expression.

    To measure inversion, a Bba_R0040 promoter is encoded adjacent to a top-strand donor sequence, followed by a GFP CDS and target sequence encoded on the bottom strand. Co-expression of a second plasmid encoding a bridge RNA and a T7-inducible IS621 recombinase results in the inversion of the GFP CDS (around 900?bp), yielding GFP expression.

    In all GFP reporter assays, co-transformed cells were plated on fresh LB agar containing kanamycin, chloramphenicol and 0.07?mM IPTG to induce recombinase expression. Plates were incubated at 37?°C for 16?h and subsequently incubated at room temperature for 8?h. Hundreds of colonies were subsequently scraped from the plate, resuspended in TB and diluted to an appropriate concentration for flow cytometry. Around 50,000 cells were analysed on a Novocyte Quanteon Flow Cytometer to assess the fluorescence intensity of GFP-expressing cells. The mean fluorescence intensity of the population (including both GFP+ and GFP? cells) is plotted as analysed with NovoExpress software (v.1.5.6). pRecombinant plasmids were isolated by picking GFP+ colonies under blue light, seeding in TB containing kanamycin and chloramphenicol, incubating for 16?h at 37?°C with shaking at 200?rpm, and isolating using the QIAprep Spin Miniprep kit. The isolated plasmids were sent for whole-plasmid sequencing to confirm recombination (Primordium Labs).

    Design of the oligo pool for systematic pairwise screening of bridge RNA target-binding loops and targets

    A pooled screen was designed to test target and target-binding loop mismatch tolerance and relative efficiency across diverse guide sequences. Several categories of oligos were designed to answer different questions. First, 10,656 oligos were designed to test hundreds of different target guides with single-mismatch pairs. That is, for a given target, one position in the guide and the corresponding position in the target to generate all 4?×?4?=?16 combinations of nucleotides. Target guides were selected to reduce genomic off-targets. Next, 3,600 oligos were designed to test different combinations of double mismatches between target-binding loop and target. Next, 2,000 oligos were designed as an internal set of negative controls by ensuring that none of the 9 programmable positions (excluding the CT core) matched in the target-binding loop and the target. Next, another 1,800 oligos were designed to test more single-mismatch combinations, but did not include all 4?×?4 combinations in the target and the target-binding loop. Finally, 1,610 oligos were designed to test how mismatches in the dinucleotide core of the bridge RNA sequences affected the recombination efficiency. One unique barcode per amplicon was assigned at random, ensuring that no two barcodes were within two mismatches of each other. Each oligo encoded a synthetic Bba_R0040 promoter followed by a target sequence, a unique barcode, the J23119 promoter and the first 104 bases of the bridge RNA, which includes the 5′ stem-loop and target-binding loop. The oligos were ordered as a single pooled library from Twist Bioscience.

    Cloning of the oligo pool for systematic pairwise screening of bridge RNA target-binding loops and targets

    A vector encoding the final 73?bp of the bridge RNA (the WT donor-binding loop) and a T7-inducible IS621 recombinase was digested using BsaI. The oligo library was amplified with primers encoding overhangs compatible with the digested vector for Gibson cloning. In brief, the library was cloned into the vector by Gibson cloning, and electroporated in Endura DUO electrocompetent cells (Biosearch Technologies). Hundreds of thousands of colonies were isolated for sufficient coverage of the oligo library, and plasmids containing library members were purified using the Nucleobond Xtra Midiprep kit (Macherey Nagel).

    Recombination assay with the library of bridge RNA target-binding loops and targets

    The plasmid library encoding thousands of target and bridge RNA target-binding loop pairs was co-electroporated into E. cloni EXPRESS electrocompetent cells (Biosearch Technologies) along with a donor plasmid and an inactive kanamycin resistance gene. Recombination between the two plasmids results in the expression of the kanamycin resistance gene, allowing cell survival. After co-electroporation and recovery, cells were plated on bioassay dishes with LB agar. One plating condition, serving as the control, was LB agar with chloramphenicol and ampicillin, which maintain the plasmids but do not induce or require recombination. A second condition was LB agar with chloramphenicol, ampicillin, kanamycin and 0.1?mM IPTG; IPTG induces recombinase expression, prompting recombination, and kanamycin selects for cells that have induced recombination between the donor and the target plasmid. Both conditions were performed in two replicates. Recombination indicates a compatible target–target-binding loop pair within the library.

    Hundreds of thousands of colonies were scraped from the bioassay dishes and had plasmid DNA extracted using the Nucleobond Xtra Midiprep kit (Macherey Nagel). After plasmid DNA isolation, samples were prepared for NGS. For DNA isolated from the control conditions, a PCR was used to amplify the barcodes specifying target and bridge RNA pairs to measure the distribution of barcodes without selecting conditions. For DNA isolated from selection conditions, a PCR was used to amplify the barcodes specifying target and bridge RNA pairs, with one primer priming from the donor plasmid and the other priming from the target plasmid such that only barcodes from recombinant plasmids were measured. The distribution of barcodes from recombinant plasmids was subsequently compared to the distribution of barcodes under control conditions.

    Analysis of target specificity screen

    Amplicon sequences were processed using the bbduk tool66. Amplicon sequencing data were then aligned to their respective wild types using bwa-mem, with ambiguous nucleotides at all variable positions67. Barcodes were then extracted from the amplicons using custom Python scripts. Barcodes were mapped to the designed barcode library, tolerating single mismatches when making assignments. This resulted in a table of barcode counts per biological replicate. Using custom R scripts, the counts were normalized within each replicate using counts per million (CPM), which converts raw barcode counts into barcode counts per million barcodes. CPM values were then averaged across the two biological replicates in each condition. For the recombinant barcodes, CPM values were then corrected by the control barcode CPM values using a simple correction factor for each barcode, calculated by dividing the expected barcode CPM (assuming a uniform distribution) by the observed barcode CPM. These corrected CPM values were subsequently used in many of the individual analyses. Mismatch tolerance was assessed by limiting the analysis to the top quintile of the most efficient 4?×?4 single-mismatch sets, in which each set was ranked according to the barcode with maximum efficiency, and then averaging the percentage of total CPM within each set at each position. The motif of enriched nucleotides at each position was generated by determining the nucleotide composition of the top quintile of the most efficient target-binding loop–target pairs (without mismatches), and comparing this to the nucleotide composition of the entire set.

    IS621 genomic insertion assay with long-read sequencing

    A plasmid was prepared that encoded a donor sequence adjacent to a constitutively expressed kanamycin resistance gene and a temperature-sensitive Rep101 protein. Plasmid replication of this donor plasmid was eliminated in cells upon growth at 37?°C, ensuring that cells encode a single copy of the donor plasmid. A cell line was prepared encoding this donor plasmid by transforming BL21(DE3) and making the resultant cell line chemically competent using the Mix & Go preparation kit (Zymo). The temperature-sensitive donor plasmid was then transformed with a second plasmid encoding a T7-inducible recombinase and a constitutively expressed bridge RNA. The donor-binding loop of the bridge RNA was programmed to recognize the donor sequence within the donor plasmid and the target-binding loop of the bridge RNA was programmed to recognize a target sequence in the BL21(DE3) E. coli genome. After transformation, cells were recovered and plated on 10-cm LB agar plates containing 0.02?mM IPTG, chloramphenicol and kanamycin; insertion of the donor plasmid and expression of the kanamycin resistance gene from the genome is required for cell survival. The thousands of resulting colonies, each with an insertion of the donor plasmid into the genome, were scraped from the plate. Genomic DNA was extracted from the pool of colonies using the Quick DNA Miniprep Plus kit (Zymo). Genomic DNA was then cleaned up using AMpure XP (Beckman Coulter) and sequenced using bacterial genome nanopore sequencing to at least 100× genome coverage.

    Sequencing data were downsampled to a sequencing depth of 200× in reprogrammed bridge RNA experiments, and to a depth of 1,400× in the WT bridge RNA experiments. To identify long reads containing potential insertion junctions between the plasmid donor and the E. coli genome (NZ_CP053602.1), all individual reads were programmatically scanned for the presence of the terminal 20 nucleotides of the donor sequence, excluding the core. If a 20-bp sub-sequence of a read matched the 5′ terminus or 3′ terminus (allowing for up to two mismatches), then the read was split and the flanking sequences were written to separate files. These flanking sequences were then mapped back to the plasmid sequences and the E. coli genome using minimap2 (Li 2018), and assigned as originating from the plasmid or the E. coli genome according to whichever had the higher alignment score. Reads were then assigned to specific insertion junctions in the E. coli genome to identify precise insertion sites. Insertion sites that were within 5?bp of each other were merged together using bedtools merge68 and a representative insertion site was selected. For the reprogrammed bridge RNA genome insertion experiments, additional filters were applied to remove low-quality alignments and account for a low rate (<1%) of cross-sample contamination (possibly owing to index hopping). Low-quality predicted insertion sites were excluded only if they met certain criteria: either (1) occurring at a total insertion frequency of less than 1%; occurring at a Levenshtein distance of more than 2?nt from the 11-nt target and donor; and supported by a large fraction of clipped reads (more than 25%, indicating low alignment quality); or (2) occurring at a total insertion frequency of less than 1%; occurring at a Levenshtein distance of more than 2?nt from the 11-nt target and donor; and matching a high frequency (more than 1%) and close target match (Levenshtein distance of less than 3?nt) in a different sample (suggesting that index hopping across samples is likely). The total number of reads per site was subsequently used to determine the insertion specificity for each site.

    Off-target sites were evaluated by calculating the Levenshtein distance between the 11-nt off-target and the 11-nt target and donor sequences. Sequences with a Levenshtein distance of more than 2?nt from the target and donor were further evaluated by searching for shared k-mer sequences in the 14-nt off-target, the 14-nt expected target and the 14-nt donor. To determine whether the off-target sequences were enriched for shared target or donor k-mers, the maximum-length shared k-mer distribution was generated and compared to a null distribution in which the 14-nt off-target sequences were randomly shuffled. This shuffling procedure was repeated 1,000 times to calculate the null distribution.

    A computational pipeline was developed to identify potential structural variants (50?bp or greater in size) that were independent from the donor plasmid. All long-read nanopore sequences were aligned to the BL21(DE3) E. coli genome (NZ_CP053602.1) and the pDonor and pHelper plasmid sequences. Reads that aligned to the pDonor or pHelper sequences were then excluded from the E. coli genome alignment. These filtered alignments were analysed using fgsv v.0.0.1 (ref. 69). The tool geNomad was used to annotate a structural variant involving a possible prophage element70.

    For the WT bridge RNA, REP elements were also identified and annotated to determine how frequently they were targeted. REP elements were identified by a BLAST search of three different known REP sequences collected from two different studies11,16. These query sequences were TGCCGGATGCGGCGTAAACGCCTTATCCGGCCTAC, GCCTGATGCGCTACGCTTATCAGGCCTACG and GCCTGATGCGACGCTGGCGCGTCTTATCAGGCCTACG.

    Design of the oligo pool for systematic screening of bridge RNA donor-binding loops and donors

    A pooled screen was designed to test donor-binding loop programmability, mismatch tolerance and relative efficiency across diverse guide sequences. Several categories of oligos were designed to answer different questions. Donor sequences were selected to reduce predicted genomic off-targets. First, 13,593 oligos were designed that included complete single-mismatch scans across 100 distinct donors, including all position 4?×?4?=?16 mismatches with the donor at the corresponding position. Next, 5,000 completely random donor guides were selected and paired with a perfectly matching donor for the analysis of a high number of diverse donor sequences. Finally, 2,297 oligos to test single-mismatch and double-mismatch scans of the WT donor sequence and 4 other functional donors were included. Next, 50 negative control oligos were included that ensured that none of the 9 programmable positions (excluding the CT core) matched in the donor-binding loop and donor. Each oligo encoded a partial sequence of the IS621 RE (52?bp 5′ of the CT core), the reprogrammed donor sequence and a full-length LE (191?bp) encoding a bridge RNA as found in the WT system, such that expression of the bridge RNA would be mediated by the natural promoter in cis. The donor site sequence and donor-binding loop sequence of the bridge RNA were modified in each member according to the description above, whereas the target-binding loop of the bridge RNA was constant and programmed to recognize the target sequence T5, which is orthogonal to the BL21(DE3) E. coli genome. The oligo was flanked on both ends with sequences suitable for Golden Gate cloning into a desired plasmid backbone. All oligos were ordered as a single pooled library from Twist.

    Cloning of the oligo pool for screening of bridge RNA donor-binding loops and donors

    First, a vector was constructed encoding a kanamycin resistance gene with no promoter on the bottom strand, followed by the first 61?bp of the IS621 RE sequence. This was followed by a BsaI landing pad site for Golden Gate cloning, an HDV ribozyme sequence and a unique molecular identifier (UMI) of length 12. The UMI backbone was pre-digested by BsaI and the oligo library was cloned into the backbone through Golden Gate cloning after amplification with appropriate primers, such that the full-length IS621 RE was reconstituted and the LE containing the bridge RNA was directly adjacent to the HDV ribozyme sequence. The resulting library was electroporated in Endura DUO electrocompetent cells (Biosearch Technologies). Hundreds of thousands of colonies were isolated for sufficient coverage of the oligo library, and plasmids containing library members were purified using the Nucleobond Xtra Midiprep kit (Macherey Nagel).

    Recombination assay with the library of bridge RNA donor-binding loops and donors

    The plasmid library encoding thousands of donor and bridge RNA donor-binding loop pairs was co-electroporated into E. cloni EXPRESS electrocompetent cells (Biosearch Technologies) with a target plasmid encoding the T5 target sequence and a T7-inducible IS621 recombinase. Recombination between the two plasmids results in the expression of the kanamycin resistance gene, allowing cell survival. After co-electroporation and recovery, cells were plated on bioassay dishes with LB agar. One plating condition, serving as the control, was LB agar with chloramphenicol and ampicillin, which maintain the plasmids but do not induce or require recombination. A second condition was LB agar with chloramphenicol, ampicillin, kanamycin and 0.07?mM IPTG; IPTG induces recombinase expression, prompting recombination, and kanamycin selects for cells that have induced recombination between the donor and the target plasmid. Both conditions were performed in two replicates. Recombination indicates a compatible target–target-binding loop pair within the library.

    Hundreds of thousands of colonies were scraped from the bioassay dishes and had plasmid DNA extracted using the Nucleobond Xtra Midiprep kit (Macherey Nagel). After the isolation of plasmid DNA, samples were prepared for NGS. For DNA isolated from the control conditions, a PCR was used to amplify the UMI specifying donor and bridge RNA pairs to measure the distribution of UMIs without selecting conditions. For DNA isolated from selection conditions, a PCR was used to amplify the UMIs specifying donor and bridge RNA pairs, with one primer priming from the donor plasmid and the other priming from the target plasmid such that only UMIs from recombinant plasmids were measured. The distribution of UMIs from recombinant plasmids was subsequently compared to the distribution of UMIs under control conditions. UMIs were initially mapped to donor–bridge RNA pairs by amplifying a region of the input donor library such that information about all variable sites within the full length of the RE–LE was captured in addition to the adjacent UMI.

    Analysis of the donor specificity screen

    All amplicon sequence data were preprocessed using bbduk to remove adapters. Next, UMIs were mapped to their respective oligos. This was done by aligning to the expected amplicon sequence with ambiguous N nucleotides in all of the variable positions using bwa-mem67. UMIs were then determined from the alignments, and combined with the variable LDG and RDG to guarantee the uniqueness of each UMI to each oligo. Next, control and recombinant samples were analysed in much the same way as the previously described target screen, but UMIs were counted rather than assigned barcodes. Next, UMI counts were converted to CPM, averaged across two biological replicates and normalized according to the correction factors calculated in the control condition. These CPM values were then analysed across different oligo categories to assess mismatch tolerance, how distance from the wild-type donor affects efficiency and which nucleotide sequences were favoured or disfavoured at each position in the donor.

    Additional analyses of natural IS110 sequences

    Natural IS621 target sites were extracted from the genomic sequence database by searching for exact matches to the 1,277-bp IS621, excluding the core. These target sequences were then clustered using mmseqs2 and the parameters ‘easy-cluster --cov-mode 0 -c 0.800 --min-seq-id 0.800’52. This search and clustering identified 272 distinct target sites, which were then analysed to identify a conserved target motif and compared with the experimental observed IS621 target sequences in the E. coli BL21(DE3) genome.

    A paired alignment of target sites and bridge RNA sequences was analysed to determine how the target site motif changed as the guide RNAs were varied. All aligned bridge RNA sequences that lacked gaps in the nine-base LTG and the four-base RTG were first identified. Next, only LTG and RTG sequences with CT core guides were selected. Next, only target-binding loops with more than 20 associated target sites were kept. For each of these unique remaining target-binding loops, a consensus sequence of the motif was constructed by selecting the most common nucleotide at each of the 11 target positions. If there were ties, then the position was represented by the ambiguous IUPAC character N. These consensus target sites were then compared with the expected target sites to determine how closely they matched.



日本在线不卡v二区| 人人干人人操人人爱| 奇米四色影视777久久久| 干干干天天| 老司机福利青青草| 久久无码一区二区二三区性色| 91强热人妻| 亚洲激情综合| 五月丁香成人网| 欧美精品成人一区二区在线观看| 日日夜夜天天| 高清在线偷拍自拍视频| 免费αV在线视频| www…国产操逼| 六月婷婷激情| 亚洲国产欧美另类自拍| 国产强奸乱伦第1页| 一区二区三区免费视频入口| 国产67194| 欧美精品成人在线播放| 蜜臀无码视频在线观看| 五月丁香六月综合缴清无码| 国产精品无码论坛| 国产女人91精品嗷嗷嗷嗷| 精品人体无圣光凹凸| 亚洲制服aⅴ中文字幕| 婷婷久久五月天| 天天视频黄| 手机看av网站在线看| 性videos欧美熟妇hdx| 国产兽交视频在线播放| 噜噜吧,噜噜色,噜噜| 精品伊人久久久大香线蕉小说| 日韩性爱免费观看视频| 手机在线人成免费视频| 思思热国产在线视频| 1000部熟女视频在线观看| 中日韩一区二区三区欧美| 二男一女成人A片| 亚欧视频在线| 深夜激情无码| 国产高清26uuu| 密乳无码| 成人26uuu| 无码最新| 又粗又长又大国产不卡| 26uuu国产| 日本三级韩三级99久久| 丁香六月啪啪| 国产精品熟女九色九色蜜臀| 日本在线播放不卡一区| 亚洲天堂精品日韩电影| 激情丁香婷婷| 男人天堂黄片| 亚洲Av无码成人精品国产| 综合久久久久久久久91| 天天干人人看综合| 我要去看2个日本美女.com曹逼| 中文字幕性感少妇av| 久久性视频| 人人干黄色| www国产天美久久久| 亚洲无992tv| 粉嫩av在线一区二区| 男女做爰猛烈动高潮A片免费应用 国产成人精品亚洲日本 | 色嗨嗨在线| 色综合久| 色久综合| 手机在线免费看的av| 久久五月天婷婷| 黄片com.| 国产人妻天天干精品| av日韩在线观看电影| 天天色天天干天天爱| 婷婷五月成人| 色99在线| 人人操人人摸avav| 激情综合婷婷| 91色色色| 久久久久久网址| 人人考人人摸人人干| 国产内射爽爽大片| 2019精品国产无码成人| 狠狠操综合| 欧美大香蕉同搞| 日本欧美成人片AAAA| 一类无码操逼视频| 高清国产性猛交xxxx乱大交| 巨爆乳肉感一区二区三区竹菊影视| 久久久久久国产精品免费网站| 欧美一级黄片视频在线| 日本亚洲熟女视频| 秋霞免费无码视频日韩A片| 熟女激情综合网| 欧美国产欧美在线观看| 国产精品毛片?v一区二区三区| 久久久精品中文字幕爱豆| 99热精品在线在线| 天天色天天干天天爱| 亚洲无码?第一页| 家庭乱伦国产精品| 日韩无码黄色片| 五月天社区| 免费成人自拍视频在线| 久久精品久久久久久久久| 国产Aα| 精品人妻一区二区三区四区石在线| 久久性爱视频| 青娱乐淫乱1314| 一个国产在线综合网站| 人妻天天爽天天爽三区| 国产精品人妻无码久久久互動交流| 日韩免费福利在线观看| av天堂精品久久| 亚洲AV成人精品网站在AV| 欧美精品久久96人妻无码| 亚洲成人性爱网站在线播放| 强奸xx国产| 日韩激情中文字幕有码| 99亚洲精品| 色欲三区| 日韩中文字幕视频在线观看 | 国产精品探花色| 欧美成人一级麻豆| 1769国内精品视频| 超碰在线人妻中文字幕| 无码精品久久| 91影库| 五月天大香蕉| 天天干人人干天天日97| 国产极品粉嫩馒头一线天av| 国产操伦| 国产91亚洲精品一区二区三区| 欧美一级A片在线看视频性色| 国产综合日韩伦理| 强歼乱伦资源网| 成人八戒网站| 91精品久久久久久综合五月天| 久久久久久99AV无码免费网站| 日韩少妇一区二区三区| 激情小说日韩无码| 婷婷久久五月| 欧美偷拍区| 开心五月天激情网| 亚洲精品第一| 性色av网站| 在线视频亚洲无码| 99久久久久久亚洲精品不卡| 欧美性爱伊人| 久久久精品国产亚洲AV无码| 九九亚洲| 成人AV在线电影| 婷婷av在线中文字幕| 日本一区二区中文字幕久久| 久久东京国产精品视频| 少妇无码av专区线| 99热最新| 欧美精品日韩一区二区| 免费毛片在线播放| 超碰99在线观看| 乱老女人一区二区视频| 天天干天天操天天操夜夜操天天操| 日韩免费在线观看不卡| 婷婷色婷婷| 在线日韩精品一区二区三区| 国产精品一级毛片不卡视| 亚洲国产精品无码AV久久久| 亚洲国产一级黄色视频| 久久99精品视频| 亚洲一区二区久久久久| 国产51色综合久久免费| 高清国产性猛交xxxx乱大交| ?亚洲伊人伊成久久人综合网 | 国产 日韩,欧美 自拍| 天堂69亚洲精品中文字| 99热9| 在线无码视频| yiqicaoav| 日韩熟女三十乱伦| 欧美巨大性舒爽顶到了| 亚州操逼网| 色欲天天婬色婬香WWW夜色| 亚洲 欧美 精品专区 极品| 欧美日韩精品一区二区三区高清| 亚洲 一区二区 自拍| 国产男人又猛又粗又爽| 国产 丝袜 欧美中文 另类| 国模艳艳啪啪一区| 99re国产中文字幕| 乱论91| 五月天婷婷综合网| 黄色av一区二区在线| 无码日韩网站| 最新av网站在线观看| 国产精品一二三区18| 美国精品国产精品| 性久久| 51久久夜色精品国产麻豆| 天天草天天日| 无码少妇精品一区二区60岁老人| 色五月婷婷五月天| 黄片com.| 水多多映视AV| 久久久久久精品免费看A级| 免费观看网黄| 国产精品熟女九色九色蜜臀| 色婷婷基地| 狠狠综合| 久久精品国产亚洲5555| 中文字幕av亚洲在线| 久久精品国产欧美日韩亚洲欧美日韩中文久久国产一区 | 国模少妇一区二区三区| 三级激情网站| 日本综合色图| 免费观看成人www精品视频| 国产日韩在线播放| 99热只有| 国产熟女乱论| 亚洲Av无码成人精品国产| 99啪啪视频| 国语av狠狠色丁香婷婷综合激情| 2018色综合天天操| 伊人五月天| 依人大香蕉| 强奸乱伦中文字幕AV| 黄色操人| 国产激情视频一区区三区| 色欲天天综合网| 色综合网1| 日产操逼| 奇米四色影视777久久久| 黄色小视频日本txt| 日韩欧美亚欧在线视频| 国产极品精品美女视频| 久久久久国产精品片区无码直播 | 人人人干干人人干| 国产精品盗摄 偷窥盗摄| 国产精品一级特黄aaa大片在线观看 | 亚洲一区二区三区中文字幕| 亚洲成av人片色午夜乱码| 精品少妇人妻av久久免费| 99久久婷婷国产综合| av在线播放国产一区| 五月天丁香婷婷综合网站| 精品国产一区二区三区四区在线看| 乱性AV| 99精品在线| 伊人久久综合影院精品久久久| 凹凸视频在线观看伊人| 成年人性爱日韩| 天天综合色| 精品一区二区三区国产| 免费av高清无码| www.91理论| AV污污污污| 婷婷亚洲综合| 东京热一区二区中文字幕| 干干干天天| 国产精品熟女九色九色蜜臀| av无线看| 日韩操p| av资源在线播放天堂| 立川理惠被中出无码| 一级性爱啪啪视频| 看大黄色大片原件| 日韩精品在线放| 很很操在线| 精品一区二区亚洲国产| 国产91福利小视频在线观看| 亚洲男人久久综合天堂| 超碰99在线| AA丁香综合激情| 五月婷婷丁香中文字幕| 麻豆这里只有精品| 久久精品 六十路 熟女 欧美| 99婷婷| 久久性爱城| 97人人模人人爽人人| 久久久久久久亚洲Av无码| 韩国黄片aaaa| 人妻一区视频| 91在线视频国产网站| av网站免费线看| A一区片| 99久久亚洲精品无码毛片潘甜甜| A V视频日本| 高清无码一区二区三区| 97人人模人人爽人人| 伊人久久综合影院| 色香色欲天天综合网天天来吧| 婷婷五月成人| 国产在线强奸视频| 51一区二区三区| 日本 欧美 国产一区| av网站在线看| 北约熟女超碰| 国产强奸乱伦欧美| 欧美日本成人一区二区| 天天爽夜夜操| 免费一级精品啪啪视频| 午夜一区| 樱花蜜乳av| 日韩99精品视频综合区| 天天看精品动漫视频一区| 欧美亚洲日本激情在线| 99精品在线播放| 黑人精品一区二区在线播放| 欧美精品成人一区二区在线观看| 中文日本免费高清| 在线看的av| 国产免a费看黄片在线| 五月天色色网站| 国产精品麻豆免费视频| 国产精品一级特黄aaa大片在线观看| 欧美v亚洲v综合v国产v妖精| 91操人| 岛国视频一二三区| 久思思热视频在线观看| 极品国产内射| 婷婷五月色| 九九视频黄色片| 亚洲成人综合在线| 国内毛片热久久思思热| 性爱精品一区| 99热精品在线观看| 人妻中文在线| 久久久久久AV无码免费网站| 精品中文字幕一区二区l - 百度| 91丨九色丨大屁股| 色婷婷丁香| 国产日逼视频| 亚洲精品人妻在线| 亚洲AV无码AV吞精久久久久| 激情综合网激情五月天| 日韩成人精品视频自拍| 狠狠操综合| 国产精品熟女乱伦| 能在线播放的国产三级| 福利在线观看一区二区| 国产高清精品福利| 国产日韩欧美操逼视频| 岛国片国产成人亚洲播放| 国产av激情无码久久天堂| 欧美成人性爱视频在线播放| 十八禁av无码免费网站APP| 色五月综合网| 综合网色| 中文字幕 码精品视频网站| 抽插无码高清一区| 久久肏大逼| av天堂5| 啪啪啪东京| 久久精品国产欧美日韩亚洲欧美日韩中文久久国产一区 | 国产无马av| 人妻 中文 日韩| 国产一区二区av综合| 色婷婷狠狠18禁| 91人精品妻入口| 色五月激情网| 日人妻视频91| 妺妺跟我一起洗澡没忍住| 在线观看亚洲专区| 国产一区二区在线电影| 五月天激情综合网| 欧美一级美片在线观看免费| 亚洲图片偷拍视频区| 国产亚洲精品美女久久久m| 久久婷婷电影网| 欧美性爱伊人| 五月天社区| 亚洲操逼视频网站| 亚洲人妻中文在线视频| 91粉嫩萝控精品福利网站_精品影音先锋国 | 久久久久久国产手机AV| 国产69精品久久久久99尤物| 色综合一本| 九色视频91| 久久xxxx| 国产精品分类在线观看| 欧美日韩在线视频网站| 99久久九九| 亚洲成人激情小说视频| 视频分类 国内精品| WWW.操逼.COM| 秋霞曰韩R级| 操逼操逼操| 思思热国产高清| 国产精品美女在线一区| 国产福利视频精品视频| 日本色色色色色视频| 精品国产91av一区二区三区| 黑人精品成人一区二区三区 | 久久性爱网站| 国产精品制服丝袜清纯唯美| 在线国产一区二区av| 91日本在线观看| 国产精品一级片在线看| 国内毛片婷婷六月色| aa片毛片| 1024人妻熟女一区二区三区| 久久九九99| 亚洲AV无线| 精品高清av中文字幕| 伊人五月天| 亚洲婷婷综合网| 欧美性爱在线无码| 日韩三级伦理中文字幕| 久久久久国产一区二| 婷婷大香蕉| 大色综合网| 五月天综合在线| 欧美久久人人网| Sekablack无码一区| 99热在线观看| 粉嫩av在线| 国产精品久久久久久久毛片1| 欧美大香蕉同搞| www亚洲免费| AV无码久久久精品| 一级黄色性爱裸体视频| 亚洲第2页| 国产高清成人传媒影视| 丁香六月婷婷久久综合| 色哟哟综合| A级在线视频| 亚洲97成人在线观看| 国产 三级自拍| 久久人妻视频| 91热爆在线| 欧美性爱五月天| 欧美强奸乱| 国产免a费看黄片在线| 东亚亚洲无码高清| 国产亚洲精品美女久久久m| 五月丁香综合啪啪| A片大香蕉在线| 日韩欧美操逼xxx| 欧美高清无码免费视频高清版| 日韩视频精品在线观看| 五月丁香社区婷婷日韩欧美精品影院| a'v在线资源| 欧美日综合| 欧美的性爱网站免费| 色官网色综合| 日本成人免费一区二区三区| 一区二区三区 日韩欧美| 另类一区| 一级免费啪啪片| 丰满人妻大屁一区二区| 亚洲?V高清一区二区三区尤物| 欧美不卡五十路| 日本色色色视频| 高清在线不卡一区二区 视频| JIZZJIZZ国产精品喷水| 亚洲精品一区二区三区在线播放| 996热| 顶级丝袜熟女一区二区三区| 国产强奸AV在线| 国产久久免费精品视频| 97av在线视频| 天天插夜夜操| 国产精品爽爽v| 激情六月天| 中文字幕久久婷婷丁香五月天| 婷婷久月| 国内毛片无遮挡国产| 久久怡红院| 九九RE视频在线精品| 26uuu国产| 综合一区中亚洲国产成人综合精品| 亚洲制服欧美另类内射| 97任你吞精| 欧美爆乳精品一区二区| 久久精品性| av在线播放国产一区| 香蕉人人操tv| 1024人妻| 久久精品熟女亚洲AV麻豆软件| 五月香婷婷| 强免费黄色网址| 欧美性爱五月天| 欧美日韩插逼视频| 一区二区三区免费视频入口| 91精品丝袜久久久久久| 欧美成人一级麻豆| 天天插天天射| www.色婷婷| 一区二区三区黄色片a| 激情五月天色色| 艹精品| 日韩乱插| 成人自拍三级在线观看| 亚洲精品啪视频| 啪啪性爱免费视频| 国产亲戚伦亲在线| AV一起草在线| 国产精品自拍xxxx| 久久草草亚洲蜜桃臀| 91麻豆va国产精品| 久久亚洲精品成人av| 尤物网站91| 人妻另类 专区 欧美 制服 | 大色综合网| 97在线视频观看| 免费黄色A片| 国产三级多多影院2022国产AA一级毛片无码 | 婷婷伊人五月| 亚洲国产精品无码AV久久| 午夜精品久久久久久久99蜜桃一| 五月激情天| 色婷婷激一区二区三区| 久久精品久久久久久久久| 久久久九九网站| 欧美日韩在线视频网站| 激情久久av一区av二区av| 国产精品人妻免费精品| 黑人狂躁日本妞一区二区三区 | 性开放中文AV高清无码免费看| 正在播放国产精品一区| 免费在线看黄片av| 日韩无码一级黄色av片| 亚洲福利影院一区久久| 性爱精品一区| 久久久精品中文字幕爱豆| 探花熟女,姿勢到位,體驗感也到位| 三级日本一区二区三区| 凹凸久久人人| 黄片免费看黄片免费看| 色婷婷丁香五月| 国产精品久久久久久无码红治院| 欧美亚洲日本视频久久久| 黑人免费福利视频| 久久精品店| 久久久久久99AV无码免费网站| 欧美色性爱| 欧美aa一级片| 手机在线免费看的av| 午夜美女福利视频| 99热在线观看| 欧美婷婷| 欧美综合国产精品久久丁香| 亚洲字幕一区二区| 乱伦强奸区日韩| 久久只有精品| 桃花色综合影院 | 国产91乱伦| 国产在线视频二区| 国内精品嫩模A∨私拍小视频| 99视频只有精品| 亚洲人体视频在线观看| 欧美成人一级麻豆| 亚洲日韩国产欧美综合v| 亚洲无码AV九九九| 日本午夜福利影院| 秋霞成人做爱| 国产精品亚洲高清在线| 啪啪啪精品视频| 欧美日韩精品一区二区三区高清 | 婷婷色中文字幕| 日本女优在线视频福利| 日本熟女中文字幕一区| 强奸乱伦动态污图免费| 性久久久| 国产亚洲精品自在线亚洲情侣| 激情五月综合网| 51国产午夜精品视频| 成人AV超碰免费在线| 日本久久精品| 亚洲成人福利电影免费| 亚洲性爱无码乱伦av| 国产精品一区二区 尿失禁| 1000部熟女视频在线观看| 手机在线免费看的av| 久久国内| 精品无码一区二区三区色欲| 欧美日韩另类在线播放| 五月婷婷综合激情| 人人操人人摸人| 丁香六月激情| 91free福利| 成人日韩3| 欧美午夜一区二区三区| 亚洲精品a人片在线观看视| 91社区拍啪人妻| 国产强奸91| 大香蕉强奸乱伦| 亚洲精品国产熟女| 国产精品熟女乱伦| 午夜男女爽爽爽在线视频| 尤物网站91| 人人操人人插 - 百度 - 百度| 中日韩熟女| 亚洲 欧美 日韩 国产一区二区| 日本黄大片在线观看视频| 日本人妻中文字幕| 偷拍偷窥与盗摄视频专区| 极品尤物女神在线观看| 日本亚洲熟女视频| 激情第四色| 日韩不卡在线一区二区| 久久精品国产亚洲AV片多多| 免费观看日本操逼视频| 亚洲午夜未满十八勿入网站日本又色又爽又黄 | 永久免费av无码网站国产app| 午夜120视频在线观看| 国产成人无码a| 日日日日做夜夜夜夜无码| 十八禁视频网站| 蜜乳性色无码专日粉嫩骚逼AV | 中文字幕狠狠玩| 日韩成人无码| 国产精品高清2021在线| 亚洲蜜桃V妇女| 亚洲精品亚洲人成在线麻豆| 日本一区视频在线观看| 91一起操| 成人精品在线免费视频| 九九热免费视频| 国产成人精品必看| 强奸乱伦日韩AV| 日韩中文字幕视频在线观看| 1769成人国产精品视频| 久久久久久久久成人av解说| 人妻22p| 国产无码精品成人| 亚洲第一免费视频| 久久男人精品| 日韩色欲久久一二三四区| 国产农村一一级特黄毛片| 操逼无码操逼| 欧美日韩国产高清在线一二三区 | 精品亚洲成人免费在线| 色色色欧美| 黄色一级视| 乱伦av麻豆| 99热只有| 边做饭边操逼逼| 日本操嫩b网| 无码精品蜜桃一区二区三区ww| 久热在线精品免费观看| 波多野结衣被操50分钟免费视频 | 国产无马视频| 日本www操操操| 亚欧高清v| 日韩中文字幕国产| 操逼逼无码| 大香蕉av在线| 亚洲av综合色区图片亚洲| 欧美性爱精品一区二区| 麻豆三极片| 免费无码婬片AAAA片直播色戒| 国产97色在线| 欧美日韩日产免费网站看| 中文字幕亚洲永久精品| 亚洲激情综合| 人妻aa| 一区二区三区视频国产免费| 久久99热这里只频精品6学生| 亚洲天堂无码| 91人妻尻屄视频| 五月婷婷六月天| 操操吧亚洲乱伦视频| 狠狠狠狠狠狠| 欧美经典一区二区三区| 国产精品一区二区a| 女人高潮大叫一级毛片| 色青青久久影视| 国产成年免费大片黄在线观看| 动漫爆乳3D奶水一区在线观看| 亚洲中文字幕三级在线| 一区二区三区黄色片a| 国产h片在线观看视频| 精品国产一区二区三区香蕉欧美| 日韩黄片视频试看| 大香蕉久| 日韩欧美亚洲一区二区三区影院| 91九色蝌蚪在线观看| 国内亚洲精彩视频在线| 五月婷在线| 欧美日韩操操操| 色波多| 久久无码成人| 日韩激情电影中文字幕| av一区二区三区 中文| 亚洲超碰在线| 天天草天天干天天日| 极品尤物女神在线观看| 婷婷丁香六月| 中文字幕av亚洲精品| 欧美日韩99| 尤物视频一区| 国产a级午夜毛片| yiqicaoav| 日韩无码专区| 亚乱色| 婷婷精品视频| 国产www色在线观看| 激情五月天婷婷| 日日日日做夜夜夜夜无码| 亚洲性爱电影| 国产精品一区二区校花| 成人性爱电影网| 美中日韩无码| 亚洲在线网站| 国产福利av精彩对白| 人人透人人操| 久久精品国产Aⅴ| 日韩亚洲中文字幕在线| 色婷久久| 国产一区二区精品久久99 | 国产呦精品一区二区三区下载| 丁香五月天啪啪| 亚洲欧洲无码bt精品合集| 性久久| 91在线无码精品秘 软件| 久久首页| 蜜乳中文字幕a在线| 国产精品自在自拍视频 | 日韩性爱视频在线免费观看| 国产亚洲中文不卡二区| 2020中文在线一区二区三区| 2024年最新色情网站在线观看| 狠狠激情综合狠狠操中文字幕 | 久久久一区二区三区三州| 一线黄色免费性爱片| 午夜福利激情在线视频| 操逼逼无码| 欧美日韩午夜精品一区二区三区| 婷婷色婷婷| www超碰| 五月婷在线| 综合色色婷婷| 五月天丁香婷婷综合网站| 亚洲欧洲另类| 一块操欧美性爱| 欧美日韩国产三级黄色| 九九热精品| 人妻另类 专区 欧美 制服| 国产精品69人妻无码久久久| 五月激情啪啪| 一区二区三区免费视频入口| 五月婷婷综合网| 五月婷网站| 91狼人| 久久99午夜精品一区人妻| 国产精品一区二区久久精品| 成年女人黄网站| 日本视频在线观看污污污| 色婷婷五月综合激情中文字幕| 国产超碰AV在线精品| 水多多映视AV| 欧美爆乳精品一区二区| 色婷婷激情| 插欧洲美女欧美精品| 性在久久久久久| 色婷婷五月综合| 亚瑟国产精品久久无码| 四虎国产精品永久入口| 欧美日韩操逼动图| 午夜美女福利视频| 中文AV制服乱伦| 五月丁香六月综合缴清无码| 国产精品女久久久久av爽| 日产操逼| 思思在线免费视频| 亚洲制服aⅴ中文字幕| 欧美性爱1080p| 欧美岛国精品在线观看| 夜夜草网站| 亚洲欧美中文日韩视频中国语| 婷婷综合久久| 最新亚洲人成网站在线影院| 2020久久免费视频| 久久 国产 无码| 91小视频| 免费又黄又裸乳的视频| 思思热免费视频观看| 免费观看成人www精品视频| 91站街按摩店老熟女熟女| 色婷婷av在线观看| 国产不良强奸视频免费看| 伊人天天久久动态图| 日韩亚洲精品一区二区| 国产日韩精品无码去免费专区国产 | 3028国产精品| 久久av无码| 国产亚洲日韩在线三区黑人| 免费看片黄| 亚洲无码电影久久久| 国产av美女被艹的乱叫| 国产综合网站在线播放| 97干在线| 日本成人在线不卡一区二区三区| www99热| 五月天激情综合网| 欧美性爽xyxOOOO| 丁香六月婷婷| 天天草天天日| 岛国在线国产| 黄色一区二区秘书性感| 91久久久久久| 成人一道本免费视频| 一区二区免费电影久久| 美国精品国产精品| 欧美日韩成人在线| 91狠狠色丁香婷婷综合久久精品| 日本孕妇一区二区视频操逼免费看| 中文字幕五月婷婷免费| 全球成人中文在线| 激情五月婷婷综合| 在线毛片片免费观看| 99久久婷婷国产综合精品草原| 被男人添B超爽视频| 国产精品com| 国产精品一区二区 尿失禁| 欧美顶级黄色大片免费| 中文字幕精品探花视频| 一级片在线观看高清无码| 韩日无码在线观看| 强被迫伦姧在线观看无码网站| 婷婷丁香熟妇综合网| 亚洲日韩视频二区| 成年人黄色小视频网站| 日人妻视频91| 国产精品无码久久久久2025| 国产精品嫩草久久久久| 一级黄色性爱A级片| 被男人吃奶很爽的毛片| 黄色无码高清黄色无码网站| 色情综合| 亚洲精品性爱片| 国产AV久久久蜜爱影集| 黄片免费久久久久久久| 丁香婷婷五月| 88xx成人精品视频| 国产乱伦一二三区| 亚洲熟女中文字幕在线| 婷婷色在线| 亚洲婷婷丁香在线| 欧美呦呦性爱| 天天看精品动漫视频一区| 婷婷综合视频| AV一起草在线| 一区二区三区亚洲| 欧美性爱www免费版| 午夜成人爽爽爽爽A片李冰冰 | 91在线视频观看国产| 久久永久无码人妻视频| 美国aaaaa一级黄片| 大香蕉一级黄色片久久| 日韩一级欧美一级国产一级台湾 | 啪啪啪精品视频| 日本孕妇一区二区视频操逼免费看| 久草福利在线资源站| 国产一进一出视频网站| 激情干在线| 中文字幕狠狠玩| 久久人人爽爽人人爽人人片αV| 五月天婷婷成人网| 色综合久| 二男一女成人A片| 丁香五月自拍| 日韩av不卡在线看| 午夜男人一级A片7777| 日韩日本欧美在线观看| 免费无码婬片AAAA片直播色戒| A片大香蕉在线| 26uuu国产| 久久一区二区三区入口| 亚洲乱码国产乱码精网站| 五月丁香狠狠爱| 天天搞在线综合网| 免费看黄视频亚洲网站| 天天插天天射| 东京热视频网| 色九九综合| 色综合色色| 中国少妇XXXX做受| 大香蕉久操| 91精品亚洲内射孕妇| 日韩久久激情精品| AV网站高清无码在线观看| 插欧洲美女欧美精品| 成人短视频在线观看| 精品亚洲国产成人av网站| 色爱综合网| 欧美成人午夜免费福利785| 日本伦乱九九九综合| 色一色综合网| 国内亚洲精彩视频在线| 国产人伦a片信息免费片| 激情网色| 涩涩这里只有精品视频| 97在线视频观看| 果冻国产精品麻豆成人av| 亚洲国产麻豆一区二区三区| 亚洲精品一区二区日本| 91小视频| 99精品久久久久久久婷婷蜜桃| 久99| 日韩中文字幕精品一区在线| 91精品国产高清久久久久久,亚洲成人| 91成人高清在线观看| 国产欧美一级在线观看| 中文字幕在线高清男人的天堂| 看一级特黄a大一片| 日韩欧美俄罗斯A片| 这里有精品| 欧美特大AA级黄片| 国产乱伦视频污| 精品无码久久久久久久杏吧| 国产成人www免费人成看片| 激情综合五| 狠狠色色| 99九九久久| 99啪啪视频| 超碰在线欧美性爱激情| 久久这里是精品| 天堂种子在线www网资源| 久久亚洲AV无码专区国产精品| 国产在线播放成人免费| 亚洲精品无码少妇久久| 99热欧美| wwwss在线观看| 九九色精品| 国产一区二区av综合| 色牛牛AV| 激情久久久| 99综合网| 99国产精品视频尤物| 欧美人妻一区| 凹凸视频在线观看伊人| 国产精品嫩草影院午夜两性| 黄片视频观看| 四虎午夜影院| 一区二区三区亚洲| 国产女人和拘做爰视频| 乱伦一二三| 亚洲av无码成人精品国产| 国产激情在线| 国产精品久久久久久 百度| 欧洲黄色网| 久久黄色视频一区二区三区| 亚洲天堂中文字| 大香蕉久| 久久精彩视频| 天天草天天干天天日| 色哟哟AⅤ| 翔田千里A片一区二区| 久草婷婷| www.yeyecao| 私人尤物在线精品不卡| 久久社区一区二区三区| 久久久九九网站| 亚洲国产精品99久久久| 国产人妖的免费的视频| 欧美日韩色综合网| 操操吧亚洲乱伦视频| 日逼国产| 草草影院最新网址| 99热这里是精品| 香港澳门日本三级网站| 亚洲无码国产精品久久| 天天视频黄网站| 永久免费观看的毛片的网站| 韩国午夜理伦三级好看| 日本色色色色色视频| 中日韩欧美精品无码AⅤ一区二区| a在线观看| 亚洲AV成人精品网站在AV| 成人黑料社久久| 国产性久久久| 视频一区二区三区精品| 亚欧视频在线| 色五月综合网| ..日韩av毛片精品久久久 | 亚洲一欧洲中文字幕在线| 久久婷综合| 丰满高潮18xxxx| 五月丁香久久| 九九色精品| 九月激情婷婷| 国产精品久久久久综合| 一级黄色牲爱A级片| 亚洲成人免费中文字幕| 激情无码日韩| 啪啪性爱免费视频| 亚洲成a人片在线观看中文!!! | 欧美性爱十八禁| 偷拍三区| 国产亚洲精品第一最新| 五月丁香激情综合| 亚欧色图在线激情| 五月花婷婷| 精品人成视频在线观看| 人人考人人摸人人干| 国产日韩欧美操逼视频| 天天日天天色| 久9re热视频这里只有精品| 99国产在线 精品 视频| 青青伊人这里只有精品| 婷婷九月色| 一区二区三区黄色片a| 国产第25页在线观看| 91粉嫩萝控精品福利网站_精品影音先锋国| 国产老太乱伦一区| 欧美精品三区| 国产免费一区2区3区| 国产久久久久久久久一区二区| a天堂视频| 日韩一级欧美一级国产一级台湾 | 欧美中文字幕日韩在线| 午夜福利在线合集| 99热精品在线在线| 日本人妻丰满熟妇久久久久久| 免费a级毛片av无码久久精品中文字幕 | 午夜精品久久久久久久久久蜜桃 | 日韩欧美成人午夜福利| 大肥女高潮bbwbbwhd视频| 人妻人人做人人澡人人爽欧美一区 | 精品国产丝袜一区二区三区乱码| 岛国在线一区二区三区| 性爱边摸边日免费AV| 免费精品人妻一区二区三| 精品超碰国产| 尤物黄色在线观看网站| 水多多映视AV| 99久在线精品99re8热| 欧美日韩在线国产在线| 成人无码影片视频在线| 久99久视频| 久久色情| 丰满美女一级毛片在线播放| 亚洲 国产 精品一区| 中文字幕一区二区免费在线| 五月天我淫我色av| 在线毛片片免费观看| 婷婷在线播放| 操逼国产免费| 操国产逼| www.五月天| 日韩免费av片高清无码| 一区二区免费电影久久| 日韩操人| 日韩免费簧片| 免费av在线播放二区| 欧美的性爱网站免费| 国产精品亚洲免费| 日本人妻中文字幕| av强奸乱轮| 欧美日韩第一页| 欧美日韩国产一区二区小黄片大全| 久久九九视频九九视频| 久操热线| 三级色综合| 日本一区视频在线观看| 日韩亚洲国产视频| 91免费看一区二区三区| 欧美日综合| 18啪啪手机免费性爱| 岛国毛片在线观看免费| 日韩黄片影院| 亚洲最新中文字幕免费| 亚洲中字幕日本一区二区三区| 欧美 精品国产制服第一页| 肉动漫无遮挡h在线观看| av天堂5| 日韩av在线免费网站| 婷婷视频在线免费观看| 中文字幕日韩精品久久| 人人操人人爽人人操人人| 亚洲一区二区精品福利| 免费一级a毛片久久久久久鸭绿欲| 激情接吻视频久久久久久| 超碰成人最新最好看| 免费人成?大片在线播放| 日韩在线国产字幕| 久久 国产 无码| 亚洲日韩成人性爱视频| 高清国产精品福利网站| 放黄片放3级黄片没穿衣服| 岛国在线免费视频| 丰满人妻一区二区三区大胸懂色 | 中文字幕日韩电影人妻| 精品伊人久久久大香线蕉小说| 91中文在线| 精品区国产区一区二区三区| 精品久久人妻成人网| 日本www操操操| 成人三一级一片aaa| 国产97视频免费观看| 人妻熟女av国产网站| 亚洲一区二区在线观看91| 亚洲图片欧洲图片aⅴ| 青娱乐欧美激情一区二区| 亚欧国产无码精品在线| 一级日本牲交大片好爽在线看| 69XX一中文字幕人妻91| 欧美1区二区三区公司| A片A5445444| 久久久久九九九| 懂色av色欲av蜜臀av| 天天懆天天日| 国产女人和拘做爰视频| 色五月综合| 黄片免费日韩| 国内外色色色色色成人视频| 免费无码国产精品v片在线观看| 丝袜人妻av一区二区| 国产无码三级视频在线观看| 无码国产精品午夜不卡(| 一起草欧美| 日韩性爱网址| 一块操欧美| 欧美成人四级在线播放| 免费在线黄片视频| 二色av| 99精品热| 亚川综合视频| av国产无码| 91超级碰碰| ..日韩av毛片精品久久久 | 五月丁香成人网| 国产曰批免费观看久久久| 91干熟女| α√在线| 天天肏夜夜肏| 97在线视频观看| 伊人热综合| 国产黄片在线免费观看| 91精品女厕偷拍视频| 99在线观看| 国产精品一区二区亚洲人成毛片| 色黄污美女啪啪啪免费网站| 操逼逼中文字幕| 操操操日本的逼| 四虎精品永久在线播放| 激情 欧美 亚洲 小说| 亚洲国产欧美中日韩成人综合视频| 日日插夜夜| 91丨九色丨东北熟女| 黄色AV免费| 色综合五月天| 久久国内| 成人AV在线电影| 天天干天天干天天| 午夜小电影在线插入淫高潮| 久久婷婷伊人| 日韩性爱再线视频| 国产树林里野战在线看| 国产黄色剧情影片麻豆免费播放| 国产91精品在线免费| 日韩精品一区二区三区色欲| 成人片在线播放| 国产自偷自拍一区| 美女视频尤物网在线看 | 日韩一性一交一A片俄罗斯| 欧美性爱在线无码| 2017av无码免费无线播| 丝袜大香蕉|