成人做爰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.



18禁超污无遮挡无码免费网| 欧美aa一级片| 国产超碰AV在线精品| 伊人一级免费黄片| 欧美一级在线观看成人| 久久精品高清AV | 国产高清亚洲日韩一区 | 激情综合二| 中文字幕乱在线伦视频中文字幕乱码在线 | 欧美性爱第一页久久| 免费啪啪啪网站18岁| www..com操老师| 先锋色眉乱伦资源| 亚欧色图在线激情| 九九黄色视频在线观看| 亚洲日韩精品一区视频在线| 翔田千里av一区二区三区| 多乙久久久久久| 在线观看成人性爱免费小视频| 国产 无码 一区二区| 国产精品一二三区18| 国产在线精品偷| 亚洲性爱成人| 精品亚洲黄色片 国产精品导航一区二区 | 无码精品一区二区三区潘金莲| 思思热在线视频免费| 日韩特级毛片免费观看全集| 一级婬片120分钟试看| 青青草天天亲夜夜操网| 大香蕉乱伦视频网| 亚洲福利影院一区久久| 一牛影视成人片免费| 黄色操人| 色欲久久久久综合网| 国产婷婷综合在线观看| 翔田千里A片一区二区| 欧美视频在线视频免费va| 国产精品无码论坛| 亚洲图片激情综合另类| 国产精品嫩草久久久久| 777AV电影| 亚洲欧美日韩中文播放| 秋霞无码av鲁丝片一区| 久久美女国产| 精品国产91内射久久| 亚洲精品国产无码高清| 美国人人操人人操| 福利社区午夜一区二区| 中日韩熟女| 精品一区二区亚洲国产| 岛国毛片手机在线观看| 韩国一级AAA| 91综合色噜噜| 在线观看国产黄色| 亚洲国产精品成人综合| 国产精品自产拍在线观看社区| 欧美午夜视频精品久久| 2024人人操人人摸| 人人操,操人人| 一级黄色性爱A级片| 我要去看2个日本美女.com曹逼| 欧美婷婷五月天| 美女视频尤物网在线看 | 日韩欧美国产一区二区三区四区| 超碰在线欧美性爱激情| 欧美色性爱| 久久久久亚洲熟妇熟女| 国产毛片久久久久久久| 美女国产一区二区久久| 亚洲欧美一区二区网址| 97在线视频观看网站| 内射老妇BBWX0C0CK| 伊人久久大香大香线蕉中文| 国产AV久久久蜜爱影集| 青青操网| www亚洲免费| 久久精品毛片免费不卡| 囯戸精品高潮呻吟旡码| 在线视频免费播放一区| 久久五十路熟女人妻| 国产尹人在线视频免费| www.亚洲成人一区| 思思热久久成人| 在线观看亚洲专区| A V视频日本| 久久美女福利是上海美女| 岛国1区2区3区在线观看| 一道α片欧美| 一区二区三区视频在线观看免费| 免费一级a毛片久久久久久鸭绿欲| 懂色天天爱天天日天天射天天澡| 艹精品| 久久久精品电影| 国模无码一区二区三区在线| 色呦呦国产精品免费看| 91精品国产91久久福利| 无码免费在线观看黄色片| 亚洲欧洲综合av在线| 日本精品不卡一二三区| 新亚洲无码| 欧亚乱色熟女一区二区| 国产精品无码成人精品| 男人女人18禁片免费看网站| 成人aⅴ一区二区三区| 天天激情干| 日韩去日本高清在| 人人操,人人插| 国产乱人伦AVA麻豆软件.| 国产精品4p在线观看| 久久久久久国产精品免费网站| 亚洲字幕一区二区| 18禁止看精品中文字幕| 天天插夜夜操| 亚洲另类小说卡通动漫| 密乳视频在线| 日韩在线一区二区| 久久久久久久免费A片国产成a人亚洲精∨品无码 | 国产精品女aA片爽爽视频| 国产树林里野战在线看| 日韩AV一起草| 2025年A片视频精品| 色波多| 丰满人妻一区二区三区性色| 国产自产自拍| 国产精品福利资源在线尤物| 国产在线视频午夜精华在| 亚洲av综合伊人久久| 粉嫩粉嫩一区性色AV片 | 人人操,操人人| 综合在线导航一区| 竹菊影视国产一区二区| 天天色综亚洲91污| 欧美影院一区二区三区| 亚洲精品美女久久久久久久久| 伊人久久大香大香线蕉中文| 91在线无码精品秘 软件| 欧美v亚洲v日韩v最新在线二区| 日韩欧美经典在线观看| 人、人、摸,人、人、草| a片在线播放| 国产又操| 顶级丝袜熟女一区二区三区| 18禁止看精品中文字幕| 久久精品色欧美aⅴ一区二区| 五月丁香社区婷婷日韩欧美精品影院 | 午夜理论片在线观看免费| 91精品人妻电影| 少妇高潮特黄A片| 免费成人自拍视频在线| 91深夜夜| 最近2019中文字幕国语免费版| 无码外流操逼视频| jazzjazz国产精品麻豆| 天天视频黄网站| 啪啪资源网| 国产日韩精品suv| 国产精品com| 亚洲精品尤物yw在线影院| 人人操人人摸avav| 日韩成人大片在线观看| 尤物av网站免费在线播放| 国产特级毛片AAAAAA高潮流水| 国产中文大片资源中文字幕| 粉嫩粉嫩一区性色AV片| 在线a亚洲视频播放在线| 国产欧美一级在线观看| 日韩黄片视频试看| 国产精品com| 国产精品女生av| 一线黄色免费性爱片| 秋霞曰韩R级| 国产日韩欧美中文在线播放| 美女黄频a美女大全免费皮| 丁香色狠狠色综合久久小说| 精品一区二区人妖| 亚洲中文字幕在现观看| 一区AV| 黄色小视频日本txt| www国产无码| 岛国色情视频在线观看| 九月激情婷婷| 日本国产亚洲一区在线观看| 人人性爱视频免费| 亚州熟女乱伦| 国模少妇一区二区三区| 老子午夜伦不卡影院| 人人操我人人干| 欧美日韩 强奸乱伦| 性爱乱伦一区| 日韩性爱高清免费视频| 国产成人欧美精品在线| 欧美系列在线一区二区| 人妻一二三区| 精品国产一区二区三区香蕉欧美| 日韩av在线播放不卡| h4610国产人妻| 青青操狠狠撩| 日韩av不卡在线观看| 色色色色日本| 殴美在线AⅤ| 成人日韩3| 97最新在线播放视频| 国产精品日本无码A片| 欧美日韩操操操| 欧美日韩性爱无码| 精品成人亚洲午夜电影| 八戒午夜福利理论片| 久久久久亚洲熟妇熟女| 男女性无套 免费九一| 韩国三级理论在线| 免费无码国产精品v片在线观看 | 成人福利视频网| 在线观看亚洲专区| 波多野结衣一级视频| 太久视频| 免费少妇一区二区| 欧美顶级黄色大片免费| 午夜免费福利视频一区| 欧美国产日韩高清在线| 久久精品国产亚洲AV先锋| 人人操人人狠狠操| 亚洲av性爱电影| 肉动漫无遮挡h在线观看| 国产精品麻豆视频网站| 黑人狂躁日本妞一区二区三区| 欧美一区二区福利在线| 欧美熟妇成人一区二区| 日韩性爱视频免费在线| 香伊人在线| 人人摸人人入| 秋霞怕怕片| 日本一级特级毛片视频| 97在线精品观看视频| 国产精品露脸在线观看| 国产一级高清免费观看| 久久双插| 操逼免费视频无码国产| 亚洲一区二区三区中文字幕| 91在线视频免费播放| 欧美婷婷五月天| 国产精品又黄又猛又粗| 操逼无码操逼| 久久精品性| 人人人干干人人干| 成人免费看吃奶视频网站| 五十路六十路素人熟女| 欧美精品1区2区3区| 国产99精品一区二区三区免费| 亚洲免费在线探花| www.99中文字幕| 色婷亚洲五月在线观看| 超碰精品国产无码| 久久久免费的精品| www亚洲免费| 免费精品中文字幕| 操日韩第| 強姦亂倫a| 色欲无码人妻日韩欧美精品| 可能人人看人人摸| 日本一级二级三级网站| 久久无码一区二区二三区性色| 午夜福利1区2区3区| 一个国产在线综合网站| 在线观看高清AV| 91网九色蝌蚪操熟女| 欧美第一页| 无码最新| 91精品亚洲内射孕妇| 97精品熟女少妇一区| 亚洲va有码在线天堂| 国产无码精品成人| 精品欧美А∨无码黑人大荫蒂| 免费一级毛片在线视频观看| 美女午夜福利免费视频| 亚洲成人在线高清| 国产精品无码AV网站| 新亚洲无码| 天天色综合天天操| 一类av片在线看| 日日日日做夜夜夜夜无码| 91在线精品| 最新日产中文在线麻豆 | 欧美顶级黄色大片免费| 强奸抽插av| 国产51色综合久久免费| 欧美性爱1080p| 亚洲av国产av综合av卡| 国产一区二区欧美日本| 国产熟女完整版中字| 国产午夜在线观看| 曰韩成人免费视频| 国产精品亚洲无码| 天天操人人操狠狠插| 精品国产91内射久久| 沈阳熟女高潮对白视频| 亚洲欧美日韩制服另类| 日韩性爱毛片操骚逼| 国产1769在线| 欧美色视频在线| 国产欧美日韩在线不卡第一页| 日本十八禁免费看污网站| 无码人妻一区二区三区色欲aⅴ| 欧美精品23| 99亚洲国产精品色一区二区三区| 国产午夜无码片在线观看影视| av最新免费中文字幕| 95自拍视频在线观看| 午夜福利视频在线一区| 亚洲无码?第一页| 国产精品国产拍高清AV| 人人操,人人插| 久久国产成人精品国产成人亚洲 | 亚洲阿v天堂无码z2018| 亚洲aV性爱| 成人麻豆av电影网站| 日va操| 国产精品ww久久| 屁屁影院一区二区三区国产| 绯色一区二区三区不卡少妇| 另类TS人妖一区二区三区| 被窝影院午夜看片无码| 亚洲伊人久久综合97| 无码日韩人妻av一| 国产女性无套 免费观看| 日本高清一区二区在线| 久久毛卡| 国产av美女被艹的乱叫| 亚洲欧洲自拍图片专区满春格| 亚欧日韩成人| 国产女人与拘做受视频免费| 免费福利视频中文字幕| 中文一区二区三区影院| 操逼片中文| 自慰白浆在线观看| 中文字幕精品资源在线| 在线a亚洲视频播放在线| 日韩在线地址一| 麻豆视频国产一区二区| 在线观看亚洲专区| 亚洲av影音先锋| 国产理论视频在线播放| 亚洲国产成人精品久久久国产成人一区二区三. | 韩国一级婬片A片无码天美| 国产99久久99热这里只有精品15| 国产亚州高清国产拍精| 95自拍视频在线观看| 99久国产精品午夜性色福利| 狠狠久久手机视频精品| 精品美女久久久久| 国产高潮AA片免费看| 国产女人和拘做爰视频| 午夜国产乱伦视频| 久久精品国内Av熟女高清| 内射老妇BBWX0C0CK| 777AV电影| 伊人激情五月天一区二区| 国产成人无码a| 欧美精品三级黄片| 国模吧 一区二区三区| 免费啪啪一级视频| 夜夜草我| 免费岛国一级片| 蜜乳AV网址| 日本黄色精品专区网站| 国产精品久久99日日| 久久精品国产亚洲AV嘿嘿| 欧美性巨大╳╳╳╳╳高跟鞋| 中文字幕五月婷婷免费| 天天视频综合在线观看视频| 日韩三A大片在线观看| 日韩欧洲操屄视频| 久久一级无码精品毛片6| 高清国产性猛交xxxx乱大交| 精品无码久久久久久久久果冻糖心| 乱伦熟妇一区二区| 亚洲好看强奸乱伦| WWW操逼| 久久东京伊人一本到鬼色| 99亚洲国产精品色一区二区三区| 日韩福利综合一区| 人妻天天爽天天爽三区| 国产综合色精品在线观看| 亚洲成人激情小说视频 | 99re免费视频精品全部| 亚洲精品一二牛牛| 在线色导航| 人妻啪| 不卡av免费在线网址| 亚洲国产精品无码AV在线| 天天躁日日躁AAA片李宗瑞| 蜜乳视频网站| 日韩av一级黄片| 欧美精品三级黄片| 在线强奷到舒服的无码视频| 1000部熟女视频在线观看| 精品久久无码午夜福利| 日韩人妻一二三区视频| 96久久精品一二三区色欲| 免费观看国产小粉嫩喷水精品午| 2018天天干在线视频| 久久草草亚洲蜜桃臀| WWW.操逼.COM| 国产激情av女片自拍| 亚洲av无线观看| 欧美性爱中文字幕无线码| 欧美一区二区三区蜜桃| 天天操天天舔| 三上悠亚在线毛片91| 国产JDAV无码视频在线观看| 91成人精品在线播放| 国产Aα| 国产成人 综合亚洲 天堂| 欧美丰满少妇xx高潮| 台湾一区国产高清在线| 秋霞网—男女啪啪亚洲免费体验区| 青娱乐国产盛宴视频| 少妇无码av专区线| 免费福利视频中文字幕 | a片偷拍视频| 93人人操人人| 欧美国产伊人久久久久| 久久91视频| 天天做天天爱夜夜爽毛片试看| 一区二区免费电影久久| 人人人人插| 无码人妻精品一区二区中文| 美国黄片aaa| 亚洲日韩视频二区| 亚洲欧美一区二区网址| 唐山老熟妇露脸啪啪叫| 两性色网| 免费操逼视频下载| 入口操逼网站| 狠狠色狠狠色狠狠五月| 国产精品人人爽人人做可爱福利| 国产 日韩 欧美 中文 另类,国产 欧美 另类 制服 变态,高清 日韩 欧美 中文,高 | 脫衣舞一区二区三区| 免费伦费视频在线观看| 欧美性爱精品七区| 欧美黄色片在线播放| 91欧洲国产成人久久精品网站| 中国操逼无码| 国产精品福利视频播放| 香蕉人人操tv| 欧美视频在线视频免费va| 国产激情在线| 亚洲精品成人动漫在线| 国产精品爽爽v| 精品国产无码中文| 欧美丰满少妇xx高潮| 午夜国产乱伦视频| 人人操人人舒服| jazzjazz国产精品麻豆| 偷窥自拍A片| 欧美精品日韩久久久九| 搡老人老9丨女老熟人| 婷婷去俺也去六月色| 亚洲黄色视频在线观看视频| 在线综合 亚洲 欧美中文字幕| 伊欧美综合视频| 午夜美女福利视频| 精品伊人久久久大香线蕉小说| A级在线视频| 岛国成人av在线播放网址| 无码137片内射在线影院| PMv在线观看| 女人天堂av在线播放| 欧美成人一级免费电影| 色欲av一区二区三区蜜芽| 国内亚洲高清无码| 国产福利av精彩对白| 亚洲福利影院一区久久| 日本熟女中文字幕一区| 呦呦影院| 久久鲁干| 偷拍 欧美 日韩| 性夜影院爽黄A爽免费动漫| 人妻色偷色噜| 中文操嬖片。| 国产黄片在线免费观看| 国产中文字幕曰本毛片| 性爱视频无打码在线观看| 无码137片内射在线影院| 亚洲精品美女久久久久久久久| 一区二区三区机械有限公司 | AV和黑人在线播放| 国产精品99久久久www| 一线黄色免费性爱片| 最新中文字幕精品在线| 乱伦图一区| 成人影院永久免费观看网址| 久久久精品电影| 亚乱色| 欧亚性爱视频免费看| 国产色产精品在线观看| 国产在线综合网| 欧美成人午夜免费福利785| AV乱伦国产| 免费黄色A片| 全球成人中文在线| JIZZJIZZ国产精品喷水| 在线啊v一区| 亚洲欧洲精品视频发布| 亚洲AV无码成人精品久久| 亚欧性爱无码| 亚洲av国产av综合av卡| 人人干人人搞人人摸| 欧美一区二区三区互相| 美国aaaaa一级黄片| 国产精品无码av在线| 久久直播国产| 欧美性爱精品一区二区| 伊人精品久久网站| 日本性爰一道本| 亚洲精品性爱片| 丁香五月综合| 日韩亚洲Av人人夜夜澡人人爽| 黄色片A级一区二区三区| 综合在线导航一区| 亚洲乱码尤物193YW| 唐山老熟妇露脸啪啪叫| 性爱乱伦一区| 亚洲精品性爱片| av在线免费一区二区| 国模精品娜娜一二三区| 国产高清成人传媒影视| 精品毛片av一区二区| 人人看欧美性爱| 日韩兔费看黄片| 亚洲一区日韩精品| 日本人妻最新在线中| 超碰国产精品无码| 少妇无码av专区线| 日产操逼| 18禁止看精品中文字幕| 国产精品无码AV网站| 日韩精品 视频一区二区| 日本理论在线| 99老司机精品视频在线观看| 男女一进一出视频久久| 国产精品视频精品一二| 偷窥自拍A片| 亚洲欧美自拍偷拍| 亚洲欧洲精品视频发布| 一块操欧美| 国产无码久久高清| 人妻啪| 乳欲人妻办公室奶水| 欧美aaaaaaa| 蜜乳av一区二区| 欧美日韩在线国产在线| 香蕉久久AⅤ...| 国产日韩欧美中文在线播放| 麻豆综合一区av| 26uuu国产亚洲综合| 永久免费观看的毛片的网站| 欧美中出1| 一区二区三区精品视频| 久久久免费的精品| 绯色一区二区三区不卡少妇| 国产福利影视| 亚洲A曰本VA欧美VA视频| 久久一级无码精品毛片6| 免费伦费视频在线观看| a一区二区三区乱码在线| 乱伦图av| 午夜福利av电影在线| 可以在线观看AV的网站| 日本操逼视频在线| 欧美第一页| 国产一区二区在线电影| AV一起草在线| 久久成年片色大黄全免费网站| 91人妻尻屄视频| 精品欧美老熟女一二区| 以及麻豆国产入口在线观看免费| 免费人成毛片乱码| 永久免费发布性爱网| 操迟操逼在巾线Fre看| 亚洲精品精品一区二区| 麻豆国产96在线| 欧洲特黄毛片免费看欧洲毛片 | 亚洲国产麻豆一区二区三区| 五月婷网站| 国产精品一区av在线| 天天躁日日躁成人字幕aⅴ| 国产精品嫩草久久久久| 不卡啪啪视频| 乱老女人一区二区视频 | 欧美强奸一区二区诱惑| 青娱乐二区免费| 秋霞无码av鲁丝片一区| 另类小说综合网| 天天射天天操天天干天天吃2018| 人人妻人人爽 97人人看碰人免费公开视频 | 天堂69亚洲精品中文字| 中文字幕久久婷婷丁香五月天| 婷婷干黄色| 日本护士高潮| 欧美岛国精品在线观看| 国产亚洲美日韩Aⅴ中文字幕无码成人| 亚洲精品日日夜夜52| 国产视频一区二区免费| 干我久操| 色就色综合| 久久久久国产无av| 久久毛卡| 国产精品久久久三级无码| 老熟女乱子伦中文字幕一区二区| 无码最新| 亚洲最新Av| 国产精品一区午夜福利| 亚洲第一黄色av网站| 午夜AV人气不卡| 亚洲婷婷综合网| 中文字幕国产在线天堂| 一区二区视频在线播放| 成人无码影片视频在线| 欧美三级中文字幕hd| 老熟女搡BBBB搡BBBB视频| 大色综合网| 台湾一区国产高清在线| 精品性爱一二三区| 秋霞免费AV| 免费观看网黄| 久操网视频| 女人 A一级| AⅤ片水多多| 91人妻素女| 操逼网站视频漫画国产| 91丨豆花丨熟女| 狠狠操夜夜| 国产绿奴视频在线观看| 狠狠婷婷亚洲中文综合久久 | 2025年A片视频精品| 26uuu久久| 强奸乱伦资源| 欧美强奸一区二区诱惑| 国产高清在线观看欧美| 国产传媒操逼视频| 性videos欧美熟妇hdx| 九九久久一区二区伦理| 18禁免费视频| 欧美性爱www免费版| 欧美日韩精品一区二区三区高清| 日韩无码操逼片| 人妻中文字幕日韩电影| 8050无码八戒| 亚洲国产奇米影视久久| 免费在线观看国内色片网站网址| 色欲蜜臀AV| 一区二区娱乐网站| 性在久久久久久| 18禁免费视频| 精品欧美老熟女一二区| 大香蕉强奸乱伦| se吧提供国产乱老熟视频胖女人| 国产 三级自拍| 亚洲AV性爱电影| 一区操逼日比视频| 免费A片三p视频| 高清国产成人无码| 精品超碰国产| 欧美精品日韩久久久九| 久久成人午夜精品影院| 激情无码日韩| 东北丰满熟女国产一区| 欧美操逼熟女| 亚洲欧美精品国产一区二区| 五月丁香六月综合缴清无码| 欧美很很操视频| A级毛片在线看免费| 国产精品久久久久久久久AV大片| 人干人人人操人人摸| 26uuu偷拍亚洲欧洲综合| 国产CHASE男男GAYGA 毛多色婷婷 久久无码一区二区二三区性色 | 亚洲精品毛片在线观看| 3d成人精品一区二区| 欧美视频一区二区在线| 被操高清无码视频| 久久亚州精品成人Av无| 精品国产肉丝袜在线拍国语| 国产精品久久久久中文字幕| 青娱乐久久艹| 国产又黄又粗又猛大片| 国产性久久久| 一区二区三区视频在线观看免费| 熟妇综合一区二区三区| 天天视频网站黄| 欧美日韩精品一区二区三区高清 | Sekablack无码一区| 国产午夜在线观看| 亚洲va有码在线天堂| 国产精品无码成人精品| 性做久久久久久免费观看软件| 中日韩久久久| 夜夜天天噜狠狠爱2021| 亚洲国产精品无码AV在线| 国产在线播放成人免费| 人妻-91porn| 久久精品国产亚洲AV成人直播| 俺去俺来也在线www| 色哟哟AⅤ| 精品人体无圣光凹凸| 亚洲精品影视老司机| 欧美91视频| 亚洲视频中文一区| 2021久久国产综合精品青草| 成人线上超碰| 国产三级在线现体验区| 日韩人成网站在线播放| 亚洲国产精品久久AV| 欧美一级特黄淫片在线观看| 亚洲无码电影久久久| 爽 好舒服 无码刺激久久| 国产极品美女高潮无套在线观看 | 亚洲古典另类欧美在线| 丁香六月激情| 无码不卡亚洲成?人片 | 亚洲高清视频在线免费观看| 亚洲第一免费视频| 午夜福利成人免费视频 | 艹精品| 久久久久人妻二区精品叶可怜| 情趣丝袜无码操逼视频| 久精品无码av一区二免费国产在线观看 | 国产精品一区二区手机看片| 天天日天天干天天色| 高潮的A片激情扒开一区| 日本熟女免费視颖| 特级丰满少妇一级AAAA爱毛片| 精品国产人成在线| xxx0国产在线播放| 啪啪啪东京| 把腿张开老子CAO烂你| 大屁股xxxxx| 在线毛片片免费观看| 在线岛国新天堂8| 亚洲中文字幕精品一区| 亚洲午夜福利在线影院| 404操逼福利视频| av黄图片在线观看| 人人玩人人添人人澡免费| 操B视频日韩无码| 欧美精品一区二区少妇免费A片| 国产成人精品日本视频| 无遮挡猛进视频免费无限观看| 夜夜影视四色| 377p欧洲日本亚洲大胆| 国产一区二区av综合| 久久产精品一区二区三区电影| 插入逼91| 性爱av在线免费观看| 一起草三级AV电影在线观看| 自拍偷拍2025在线观看| 欧美日韩国产一区二区小黄片大全| 一本久道在线综合视频| 日韩无码一级黄色av片| 91一区二区三区蜜桃| 久久亚洲AV无码白度| 亚洲五码一区二区三区| 秋霞视频一区二区| 欧美十八禁在线看| 免费福利视频中文字幕 | 五月天黄色激情视频| 亚洲伊人a线观看视频| 浪人综合网| 亚洲国产精品无码AV在线| 综合激情二| 人人 操人人 操人人| 大香蕉强奸乱伦| 久久直播国产| 欧美第一页| 亚洲va有码在线天堂| 久久黄色性爱视频| 亚洲国产欧美一区二区潘金莲| 中文字幕91页| 欧美精品三区| 久久国产免费激情视频| 欧美的性爱网站免费| 亚洲欧洲日本精品中文a∨| 噜噜噜无码AV一级一级久久影院| 黄色一区二区秘书性感| 蜜乳性色无码专日粉嫩骚逼AV | 在线国产福利网址导航| 日韩免费看黄片| 久热影视| 日韩字幕一区| 亚洲 欧美日韩 另类| 精品九九国产无码| 99久久久无码国产精品性男| 欧美国产精品久久九九| 久久久久久久强迫| 岛国片在线观看视频亚洲| 干干干天天| 综合 欧美 亚洲 日本| 无遮挡猛进视频免费无限观看| 日韩一级成人毛片免费观看| 成人精品一区二区91毛片不卡 | 女人被添高潮免费视频| 97国产高清视频在线观看| 强奸乱伦大香蕉网| 日韩成人精品| 精品国产91av一区二区三区| 深夜激情无码| 91精品丝袜在线观看| 人、人、摸,人、人、草| 一级黄色视频网| 国产免费一区2区3区| 精品区国产区一区二区三区| 99老司机精品视频在线观看| 91操人| 人人弄人人摸| 日本性爰一道本| 国产一级高清免费观看| 国内毛片无遮挡国产| 欧美日韩国产色图在线| 国产日韩欧美中文在线播放| 91精品亚洲内射孕妇| 3p国产欧美99热| 操逼不卡中文字幕| 欧美性爱第一页久久| 国产99精品一区二区三区免费| 中国一级操逼视频| 亚洲无码精品AV久久久| 日韩 国产 欧美自拍| 淫荡熟女乱伦网 | 国产激情在线观看| 我要去看2个日本美女.com曹逼| 热久久国产| 在线有码中文字幕| 高清无码久操视频| 青青免费在线视频一区| 一区操逼| 麻豆国产视频精品观看| 国产精品福利视频播放| 天天操天天干一区二区| 蜜乳中文字幕a在线| 国产成人精品亚洲日本| 国产一区二区三三视频| 在线观看一级α片刺激高潮视频| 永久免费观看的毛片的网站| 无码久久亚洲高清,| 日本操大逼| 国产1769在线| 综合欧美日本三级| 中文乱码字字幕在线第5页| 国产精品青草综合久久| 亚洲AV永久无码一区仙野| 日韩免费高清大片在线| 国产精品一区二区久久精品| 婷婷中文网| 国产乱伦亚洲| 成年无码动漫av片无尽在线 | AV污污污污| 久久美女福利是上海美女| 视频一区二区免费在线| 嫖老熟女A片一二三区| 日韩欧美亚欧在线视频| 久久riav中文精品| 激情综合网激情综合| 欧美真人抽搐一进一出gif| 美女被艹尤物视频| 色网在线视频观看免费| 精品在线观看视频在线| 久久久免费的精品| 中国东北熟女老太婆内谢| 国产高清精品一区二区三区毛片| 亚洲精品一二牛牛| 国产精品亚洲免费| 精品国产精品一区二区| 国产亚洲精品自在线亚洲情侣| 91国产操逼视频| 超碰日韩人妻| AV一区观看| 亚洲中文字幕精品久久久久久直播| 久久性爱免费送| 国产精品人妻无码久久久老鸭窝| 免费毛片在线播放| 翔田千里爆乳巨臀无码| www.黄色在线| 国产精品干干干| 国产成人精品无码久久| 国产成人 综合亚洲 天堂| 男人的天堂午夜av| 搡老熟女老女人老熟妇免费视频| 日韩亚洲中文有码视频| 久久国产成人精品国产成人亚洲| 少妇3P性爱自拍| 国产家庭乱伦网址| 日韩熟女三十乱伦| 性做久久久久久免费观看软件| 无码 黑人一区二区三区| 美国日韩黄色片| 91丨九色丨国产丨人妻在线| 精品国产肉丝袜在线拍国语| 男女猛烈无遮掩视频免费软件| a片偷拍视频| 激情视屏国产乱伦强奸| 国模不卡| 国产成人在线观看网址 | 国产白嫩漂亮KTV在线| 美女一区二区国产精品| 亚洲成人一区二区精品| 欧美激情久操网| 操逼逼无码| 啊视频在线| 91成人在线免费视频| 无码操逼网| 亚洲无码?第一页| 亚洲欧洲日韩国产自在线| 在线观看无码三级少妇| 999岛国大片| 国产无码精品成人| 国产亚洲日韩欧| 国产女人成人精品视频| 精品国产一区二区三区四区在线看| www.99色| 三级日本一区二区三区| www.91理论| av亚欧| 性videos欧美熟妇hdx| 国产久久成人| 亚洲欧美国产va在线播放频| 睡产熟女乱伦| 国产九九九九九九九九| 爱我干综合| 美女操逼A A| 国产又色又粗又黄又爽| 国产黄片在线免费观看| 乱伦图av| 26uuu国产日韩综合在线观看| 天堂中文资源在线bt| 国语精品内射在线观看| 日韩一级欧美一级在线观看| 乱伦日本色图AⅤ| 岛国大片国产| 岛国片国产成人亚洲播放| 久96热在线观看视频| 收看日本人日bb| 国产精品岛国片在线观看| 日韩精品亚洲专区在线影视| 岛国片国产成人亚洲播放| 91粉嫩萝控精品福利网站_精品影音先锋国 | 亚洲第一精品在线视频| 精品久久久久久中文字幕视频免费 | 高清无码在线播放网站| 超碰在线欧美性爱激情| 张柏芝国产一区在线观看| 中文字幕中文字幕一区二区| 日韩欧美大片免费高清啪啪| 欧美性爱系列| 国产精品成人无码av| 亚洲毛片一级带毛片基地| 热久久国产| 国产午夜在线观看| 18禁在线视频| 岛国黄| 国产高清午夜成人在线观看| 九九综合久久中文字幕| 国产精品久久久久久久久AV大片| 伊人精品久久网站| 日韩肏逼视频| 亚洲欧洲中文日韩女优乱码| 黄色电影观看久久9| 91午夜无码| а√天堂资源官网在线资源| 国产 亚洲 一二三四| 天天夜躁日日躁狠狠2002 | 国人欧美精品一区二区| 秋霞Av理论一级在线| 在线无码网站| 午夜一区二区三区国产| 日熟女| 国产亚洲精品一区二区三区| 久久久久久亚洲中文| 日韩免费大片一级播放| 久久久久久久久久久免费精品| 天天谢天天干| 国产自产自拍| 操逼操逼视频操逼| 在线性黄高清免费视频| 私人尤物在线精品不卡| 亚洲图片日本AⅤ欧美在线| 在线中文字幕| 亚洲国产成人福利在线观看| 国产对白刺激视频| 免费看一级a性色生活片久久无| 免费亚洲黄色视频在线观看| 收看日本人日bb| 久久欧美性爱视频| 亚洲色图日韩精品| 在线人人人人人人精品超| www.久久最新地址| 五十路六十路七十路熟婆| 特级丰满少妇一级AAAA爱毛片| 国产9 9在线 | 亚洲| 中文字幕在线免费观看2| 思思性爱| AV乱伦国产| 69少妇一区二区| 国产亚州精品美女久久久免费| 久操网视频| 全球成人中文在线| 人人看人人摸人人色| 正在播放国产精品一区| 去干网最新版| A V视频日本| 久久99亚洲精品久久99果| 俄罗斯一区二区视频在线观看| 欧美亚洲| 久久久久久久伊人精品| 一本大道不卡一二三区| 在线观看色视频| 亚洲综合另类小说色区亚洲成av人片在www| 国产精品不卡一区二区电影| 欧美一级黄片免费播放| 日本不卡高清视频| 国产精品熟女九色九色蜜臀 | 视频国产精品未满十八禁止在线观看| 美国人人操人人操| 国产亚州日韩欧美看片| 欧美精品久久久久久久久88| 色五天伊人| 性暴力欧美猛交在线直播| 国产成人无码网站在线视频 | 国产亚洲精品美女久久久久久2021| 岛国片在线视频网站| 国产精品无码久久久久2028| 亚洲av综合伊人久久| 亚洲精品毛片在线观看| 日人妻视频91| www.91人妻.com| 国产小黄片在线免费观看| 91艹B视频| 天天色综合天天操| 日韩精品一区二区日韩| 亚洲第一在线视频| 日韩美女操b| 亚洲中文国际强奸字幕| 欧美性爱中文字幕无线码| 男男H黄动漫啪啪无遮挡网站| CCYY草草影院地址入口| 国产午夜精品理论片一二三区区| 国产亚洲精品激情| 一区二区三区在线美女| 日语五十路和六十路亚洲国产精品| 国产亚洲精品激情| AVE乱伦| 日韩一级片在线看| 91精品久久久久五月天精品| 日韩A优精品在线观看| 亚洲欧洲网站免费观看| 中日韩免费看男女操逼大全| 国产最新小视频在线播放下载| 欧美亚洲一级在线观看| 亚洲欧美激情在线视频| 中文字幕欧美日韩三级| 色综合久久av| 欧亚无码视频| 亚洲人妻在线一区| 中文字幕一区二区三区四区在线视频 | 国产精品久久久久久亚洲色欲| 国产操伦| 手机在线视频国内精品| 久久欧美性爱视频| 户外裸露刺激视频第一区| 东京热视频网| 静品嫩模一区二区| 久久精品国产精品一区| 九月丁香综合网| 日韩欧美一级特黄大片| 双插性欧美一二三区| 国产精品久久久久综合| 人人摸人人摸人人干| 大香焦A片| 操逼网站视频漫画国产| 思思热久久成人| 无码137片内射在线影院| 国产欧美日韩精品中文| 国产三级日产三级韩国三级| 黄色香蕉视频网站一区| 在线可观看的黄色网址| 岛国黄色大片网站| 色官网色综合| 性色AV蜜色av色欲av| 日产欧美电影一区二区三区| 亚洲免费人妻在| 国产无码三级视频在线观看| 91足交| 亚洲无码一二三区| 国产精品 视频| 无码高清少妇久久| 十八禁啪啪视频| 性暴力欧美猛交在线直播| 欧美性爱18观看| 亚洲无码com| 一级黄色性爱裸体视频| 日本污ww视频网站| 亚洲欧洲无码bt精品合集| 亚洲中文字幕在现观看| www.亚洲成人一区| 国产精品三级视频网站| 激情视屏国产乱伦强奸| 在线免费观看高清无码视频| 一线黄色免费性爱片| 毛片久久| 亚洲天堂7777| 日本高清有码网址视频| 国产高清MV操逼视频| A 在线网址| 亚洲囯产精品女人久久久| 牛黄色久午久| 人妻另类 专区 欧美 制服| 天天躁日日躁xxxxx| 国产第12页| 亚洲 暴爽 AV人人爽日日碰| 91狼人| 亚洲激情在线一区二区| 国产成人免费观看在线视频| 亚洲激情综合另类男同| 自拍偷拍国产欧美日韩韩| 18禁在线视频| 亚洲最新中文字幕免费| www.伪伪| 天天看,天天做| 超碰成人公开| 日本在线不卡v二区| 人妻99p| 亚洲熟女av日韩熟女| 99久久久久久久久| 日韩欧美成人综合在线| 伊欧美综合视频| 男人的天堂 在线一区| 性爱视频啪啪啪啪| 午夜福利在线合集| 国产成人自拍视频在线| 最新av网站在线观看| 国产在线观看91精品一区| 国产CHASE男男GAYGA 毛多色婷婷| 操逼www.| 在免费jIzzjIzz在线视频| 男人的天堂一区三区| 91av熟女人妻| 韩国一级做a久久久久| 无码抄逼网| 成人性爱高清视频免费看| 亚洲最大网站av| 色色网91| 日日插夜夜| 国产精品 视频| 久草福利在线资源站| 亚洲成人黄色在线观看| 国内亚洲精彩视频在线| 激情黄色片在线观看| 能直接看AV的网站| 国产内射爽爽大片| 国产一级黄色片在线观看| 免费观看成人www精品视频| 在线观看AV片|