Quantitative identification and genomic visualization of integron-resistance gene cassettes in environmental samples

Abstract

The present invention discloses a method for identifying, quantifying and visualizing the genome of integron-resistance gene boxes in environmental samples, which is a high-throughput identification, quantification and genome visualization of mobile genetic element integrons and drug-resistant genes mediated by them based on third-generation sequencing; including: 1) identification and quantification of integrons: using the "collection-screening-verification" strategy to construct an integrase gene database and apply it to integron identification and quantification; 2) identification of drug-resistant gene boxes and genome visualization: using the preferred variable region gene box primers combined with PCR technology to obtain its DNA, using PacBio third-generation full-length amplicon sequencing to obtain the gene box diversity long sequence, identifying the integron-mediated drug-resistant gene box category based on a professional database, using a custom R script to visualize the integron-multiple drug-resistant gene box array, and quantification based on short sequence and drug-resistant gene box long fragment mapping. The present invention provides technical support for monitoring and controlling the development of bacterial drug resistance.

Description

Technical Field

[0001] The present invention belongs to the field of environmental microbial detection, and specifically relates to a method for high-throughput accurate identification, quantification and multiple arrangement characterization of mobile genetic elements in the environment - integrons and multi-drug resistance genes carried by them based on the third-generation high-throughput long-read sequencing technology. It is crucial for the study of the environmental fate, pollution proliferation and efficient control of drug-resistant genes. Background Art

[0002] The excessive use and abuse of antibiotics have accelerated the environmental proliferation and spread of antibiotic-resistant bacteria (ARB) and resistance genes (ARG). Mobile genetic elements (MGEs, such as plasmids, integrons and transposons) mediating horizontal gene transfer (HGT) of resistance genes is one of the main mechanisms for the increased risk. Among them, integrons and their gene cassettes (GC) systems, as mobile genetic elements, are one of the main internal mechanisms for the formation and spread of bacterial multidrug resistance. They can also be located on other mobile genetic elements such as plasmids or transposons, and have a higher risk of horizontal gene transfer. The basic structure of integrons consists of conserved sequences at both ends and a variable region in the middle. The expression of the gene cassettes they carry depends on the integrase gene (intI), promoter (Pc) and recombination site (attI). Under the action of integrase, they recognize the recombination site, cut, capture and integrate diverse exogenous genes in cells through site-specific recombination to load them into their own gene cassettes, and achieve gene expression with the help of their strong promoters. Integrons can carry multiple drug resistance gene cassettes, making the host multi-drug resistant, potentially exacerbating the level of drug resistance pollution in the system. For example, researchers detected dfrA1-qacEΔ1-cmlB-blaVEB1-aadB-arr2-cmlA1-blaOXA10-aadA1a from isolates of Acinetobacter baumannii, which is resistant to multiple antibiotics such as aminoglycosides, β-lactams, chloramphenicol, rifampicin and trimethoprim. Due to the wide distribution of the integron-gene cassette system and its high correlation with disinfectants, heavy metals and drug resistance genes, the integrase gene (indicating the abundance of integrons) is considered to be an indicator of artificial pollution. The traditional quantitative analysis method of detecting integrons by PCR technology based on a single primer seriously underestimates the risk characterization of integrons, and there is a lack of high-throughput characterization methods for the diverse drug resistance genes it carries. Therefore, it is crucial to construct a high-throughput, accurate identification and quantitative method for integrons and their drug resistance gene cassettes to study the environmental fate, pollution proliferation and efficient control of drug resistance genes.

[0003] In the identification and quantification of integrons, integrase genes are often used as marker genes. However, integrase gene sequences are diverse, and identification based on a single or partial sequence can easily lead to a serious underestimation of system integrase genes. Therefore, creating a complete integrase gene diversity sequence database and using it as a reference sequence set to identify integrons provides an opportunity for accurate quantification of integrons; in the characterization of integron gene cassettes, low-throughput characteristics (such as Sanger sequencing) or limited assembly efficiency (such as second-generation sequencing) limit the full analysis of integrons and the diverse drug resistance gene arrays they carry. Until recent years, with the gradual reduction of sequencing error rates of long-read third-generation sequencing technologies (PacBio / Nanopore) and the gradual improvement of big data analysis tools, accurate analysis and comprehensive understanding of integron-drug resistance gene cassettes have become possible.

[0004] Therefore, studying the types of mobile resistance genes mediated by integrons and their evolutionary mechanisms is crucial to controlling the spread of bacterial resistance. It not only expands the cognitive system of environmental resistance pollution behavior, but also provides an important scientific basis for the targeted control of high-risk resistant pollutants. Summary of the invention

[0005] The purpose of the present invention is to provide a method for high-throughput identification, quantification and genome visualization of mobile genetic element integrons and drug-resistant genes mediated by them in environmental samples based on third-generation sequencing.

[0006] The present invention firstly uses the "collection-screening-verification" strategy to construct a complete integrase gene diversity sequence database, and establishes an integron precise quantification method based on the database based on a custom script; then, the second-generation metagenomic sequencing, PCR amplification and third-generation long-read sequencing technology are integrated to obtain the diversity sequence information of the amplified products of the drug-resistant gene box carried by the integron, and a high-throughput identification and visualization method of the multi-drug resistance gene box carried by the integron is established in combination with the self-constructed script.

[0007] The specific technical solution for achieving the purpose of the present invention is:

[0008] A method for identifying, quantifying and visualizing the genome of an integron-resistance gene cassette in an environmental sample comprises the following steps:

[0009] Step 1: Construct the integrase gene database of three types of integrons, including:

[0010] 1.1: Collect all integrase gene sequences;

[0011] 1.2: Manually screen the fully matched sequences based on the sequence description information, and remove the sequences of the non- / missing integrase genes intI1, intI2, and intI3 in the description content;

[0012] 1.3: Remove redundant sequences and verify sequence accuracy, and then construct an integrase gene nucleotide sequence seed library;

[0013] 1.4: Use Diamond tool to compare the seed library sequence with the non-redundant protein library (NR) of NCBI database, and set the parameter to e value ≤ 10 -5 , based on the description, non-integrase sequences were screened and eliminated, and then the sequences were de-redundant and verified for accuracy, and then an integrase gene amino acid sequence database was constructed;

[0014] Step 2: Extract microbial DNA from environmental samples for second-generation metagenomic sequencing and perform quality control to obtain a metagenomic dataset;

[0015] Step 3: Identification and quantification of integrons: Using the metagenomic dataset obtained in step 2, based on the integrase gene amino acid sequence database constructed in step 1, a custom Python script was imported to perform identification and quantitative analysis of the three types of integrase genes I, II, and III in environmental samples. The relative abundance of integrase genes is equivalent to the relative abundance of integrons in environmental samples.

[0016] Step 4: Screen three types of candidate primers for integron gene cassettes based on literature, use candidate primers combined with PCR technology to amplify integron gene cassettes in environmental sample DNA, use pulsed field gel electrophoresis to visualize gene cassette sequence diversity and length, and obtain three pairs of primers and their amplification products with the highest diversity of gene cassette amplification products through comparison and selection;

[0017] Step 5: The PCR amplification product of the integron gene cassette obtained in step 4 is subjected to third-generation long fragment sequencing and processed to obtain the Q30 sequence;

[0018] Step 6: Use the Q30 sequence obtained in step 5 to identify drug-resistant genes, use Prodigal to predict open reading frames (ORFs), compare the ORF sequences to the structured drug-resistant gene database (SARG), and screen and identify drug-resistant genes according to the standard of "similarity ≥ 80% and coverage ≥ 75%", and then obtain the gene cassette sequences of class I, II, and III integrons carrying drug-resistant genes; use the Diamond tool to compare the ORF sequences of the gene cassettes of class I, II, and III integrons carrying drug-resistant genes with the comprehensive non-redundant database NCBI NR library to annotate all functional genes in the sequences, and set the parameter to e value ≤ 10 -2 ; Organize the (multiple) drug resistance genes and other functional gene information carried by the integron, summarize them into structured text data, and use a self-built script based on the R language gggenes package to visualize the diversity arrangement of the integron-antibiotic resistance gene cassette;

[0019] Step 7: Based on the metagenomic dataset obtained in step 2 and the gene cassette sequences of class I, II, and III integrons carrying drug-resistant genes obtained in step 6, use the contig mode of CoverM software to map the metagenomic sequences to the gene cassette sequences of class I, II, and III integrons carrying drug-resistant genes to characterize the relative abundance of gene cassettes of class I, II, and III integrons carrying drug-resistant genes in environmental samples.

[0020] Further, the step 1.1 specifically includes: selecting the "gene" sub-library of NCBI and searching and obtaining the corresponding gene sequences with the keywords of class I, II, and III integrase genes "intI1, intI2, intI3, class 1integrase, class 2integrase, and class 3integrase"; and extracting the integrase sequence from the SEARCH / BLAST module of the integrin database INTEGRALL by keyword search.

[0021] Further, the step 1.3 specifically includes: clustering the sequences obtained in step 1.2 using CD-HIT software according to 100% sequence similarity (parameter: -c 1) to obtain non-redundant integrase nucleotide sequences, and then using the Pfam protein family domain database and CDsearch tool for comparison based on the principle of conserved domain similarity to verify sequence accuracy, thereby forming an integrase gene nucleotide sequence seed library.

[0022] Further, step 3 specifically includes: firstly, using the makedb command to format the three types of integrase gene databases constructed in step 1, using the BLASTP tool to compare the environmental sample metagenomic data set obtained in step 2 with the database, and adding the annotation results that meet the requirements of "similarity ≥ 90%, coverage ≥ 50% and e-value ≤ 10 -5 " were identified as class I, II, and III integrase genes; the environmental sample metagenomic dataset was compared with the 16S database Greengenes using BWA and BLASTN tools to calculate the number of 16S genes in each sample; Diamond was used to compare the bacterial and archaeal single-copy marker gene databases to calculate the number of cells in each sample; the following formulas were used to quantify the abundance of the three types of integrase genes in the samples:

[0023]

[0024] Where: A is the abundance of class I, II or III integrase genes, which is equivalent to the abundance of integrons in environmental samples, and the unit is copy number / cell (copy / cell); N i(intI-like sequence) is the number of sequences annotated to the intI1, intI2 or intI3 reference sequences; L readsis the read length of the metagenomic sequence (bp); L i(intI reference sequence) is the length of the intI1, intI2 or intI3 reference sequence; N 16S sequence L is the number of 16SrRNA genes identified by comparing metagenomic sequencing data with the Greengenes database; 16S sequence is the average length of 16S rRNA gene in Greengenes database, where L 16S sequence The value is 1432bp; n is the number of intI1, intI2 and intI3 reference sequences mapped; N 16S copy number is the average 16S rRNA gene copy number in cells, calculated using CopyRighter software.

[0025] Further, the step 4 specifically includes: searching for literature in the literature database Web of Science using the keywords "genecassette* or intI1 or intI2 or intI3 or integron*" and "antibiotic resistance gene*", and summarizing candidate primers based on the base sequences and amplification conditions of the three types of integron gene cassette primers used or recorded in the literature.

[0026] Further, the step 7 specifically includes: outputting the relative abundance of the gene cassettes of class I, II, and III integrons carrying drug resistance genes in environmental samples in the contig mode of CoverM software, normalizing them to RPKM (reads per kilobase per million mapped reads) values, i.e., reads per kilobase per million mapped reads, and the calculation formula is as follows:

[0027]

[0028] Where: Exon Mapped Reads is the number of reads mapped to the target drug resistance gene cassette; Total Mapped Reads is the total number of mapped reads; Exon Length is the length of the target drug resistance gene cassette.

[0029] In the technical solution of the present invention, the PCR reaction is performed by a gradient PCR instrument (Bio-Rad T100, Bio-Rad).

[0030] In the technical solution of the present invention, fastp is used for quality control and filtering of raw data.

[0031] In the technical solution of the present invention, Prodigal software is used to predict ORF of the sequence after quality control.

[0032] In the technical solution of the present invention, the functional genes are annotated in the NR database by BLAST comparison, and the target drug-resistant gene sequences are compared by BLAST to identify the drug-resistant genes.

[0033] In the technical solution of the present invention, the R package gggenes is used to visualize the arrangement of the integron gene cassette resistance genes; in addition, an automated R script is optimized and configured for abundance visualization.

[0034] Beneficial effects: Based on the second-generation metagenomic sequencing combined with PCR amplification and third-generation long sequence sequencing technology, the present invention 1) constructs an integrase gene database and its high-throughput precise identification and quantitative method for integrons in various environmental samples, and 2) creates a high-throughput identification, quantification and genome visualization method for integron-resistance gene cassettes. This set of methods provides important methodological support for quantifying the risk of integron presence in environmental systems and revealing the risk of horizontal gene transfer of diverse resistance genes mediated by them. It also provides a new perspective for studying integron-mediated mobile high-risk resistance genes and their evolutionary mechanisms. It not only broadens the cognitive system of bacterial resistance pollution behavior in the environment, but also provides an important scientific basis for the targeted control of high-risk resistance pollutants. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 It is a gel electrophoresis diagram of the PCR products of the class I, II and III integron variable region gene cassettes of the total microorganisms and lactose-fermenting Enterobacteriaceae DNA in the sewage microbiome after chlorination disinfection in the embodiment of the present invention;

[0036] Figure 2 This is a schematic diagram showing the average length distribution of long read sequences of class I, II, and III integron variable region gene cassettes respectively carried by total microorganisms and lactose-fermenting Enterobacteriaceae in chlorination disinfection in an embodiment of the present invention;

[0037] Figure 3 This is a schematic diagram of the arrangement and abundance of a super integron carrying 6 ARGs in the chlorination disinfection system of an embodiment of the present invention. DETAILED DESCRIPTION

[0038] The embodiments of the present invention are described by referring to the specific embodiments and the accompanying drawings, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments.

[0039] The present invention comprises:

[0040] 1) Based on the "collection-screening-verification" strategy, a database of integrase gene diversity sequences of three types of integrons, I, II, and III, was constructed;

[0041] In step 1), the database construction method includes the following steps:

[0042] 1-1) Collect class I, II, and III integrase gene sequences in the NCBI and INTEGRALL databases, and retrieve a total of 7294 nucleotide sequences of integrase genes intI1, intI2, and intI3;

[0043] 1-2), manually screened the description information of each sequence to obtain the completely matched sequence, and eliminated the sequences containing the non- / missing integrase genes intI1, intI2, and intI3 in the description content, and finally obtained a total of 1027 nucleotide sequences;

[0044] 1-3) CD-HIT was used to remove 100% of the duplicates according to sequence similarity, and then the Pfam protein family domain database and CDsearch tool were used for comparison to verify the sequence accuracy. After sequence verification, a total of 241 class I, class II and class III integrase gene sequences were finally obtained to form the integrase gene nucleotide sequence seed library;

[0045] 1-4), the seed library sequence was compared with the non-redundant protein library (NR) of the NCBI database using the Diamond tool, and non-integrase sequences were screened out according to the description. Then, CD-HIT was used to remove 100% of the sequence similarity and the Pfam protein family domain database and CDsearch tool were used to compare and verify the final database, totaling 455 non-redundant integrase gene amino acid sequences: class I (424), class II (26) and class III (5).

[0046] 2) Extraction of microbial DNA from environmental samples and metagenomic sequencing;

[0047] In step 2), the analysis method comprises the following steps:

[0048] 2-1) For environmental samples, sample collection is carried out according to the sampling method requirements of different environmental samples (such as water, sludge, atmospheric particulate matter, soil, biofilm, etc.);

[0049] 2-2), extracting total DNA from environmental samples and extracting DNA from cultured target microorganisms;

[0050] 2-3) Take about 1.5 μg of each DNA sample for high-throughput second-generation metagenomic sequencing to ensure that the metagenomic data of each sample is above 10Gb after quality control;

[0051] 2-4) Use fastp to perform quality control on the metagenomic data and remove gene sequences with a length of less than 100 bp and a quality value lower than 20 as clean reads.

[0052] 3) Construct a high-throughput integron identification and quantification method based on the created integrase gene database;

[0053] In step 3), the construction method comprises the following steps:

[0054] 3-1) Use BLASTP to compare the sample metagenomic dataset to the database and select annotation results that meet the criteria of “similarity ≥ 90%, coverage ≥ 50% and e-value ≤ 10 -5 ” sequences were identified as class I, II, and III integrase genes;

[0055] 3-2) Use BWA and BLASTN tools to compare the samples with the 16S database Greengenes and calculate the number of 16S genes in each sample; use Diamond tool to compare the single copy marker gene database of bacteria and archaea and calculate the number of cells in each sample;

[0056] 3-3) Use the following formula to quantify the abundance of three types of integrase genes in the sample;

[0057]

[0058] Where: A is the abundance of class I, II or III integrase genes, equivalent to the abundance of integrons in environmental samples, and the unit is copy number / cell (copy / cell); N i(intI-like sequence) is the number of sequences annotated to a specific intI1, intI2, or intI3 reference sequence; L reads is the read length of the metagenomic sequence (bp); L i(intI reference sequence) is the length of the intI1, intI2 or intI3 reference sequence; N 16S sequence is the number of 16S rRNA genes identified by comparing metagenomic sequencing data with the Greengenes database; L 16S sequence is the average length of the 16S rRNA gene in the Greengenes database. 16S sequence The value is 1432bp; n is the number of intI1, intI2 and intI3 reference sequences mapped; N 16S copy number The average 16S rRNA gene copy number in the cells was calculated using CopyRighter software.

[0059] 4) By comparing the metagenomic sequencing data with the constructed three-type integron database, the abundance information of the three types of integrase genes intI1, intI2, and intI3 in the samples was obtained;

[0060] In step 4), the analysis method comprises the following steps:

[0061] 4-1) Use a custom script to import the *.fa sequence library of the three types of integrons constructed;

[0062] 4-2) Combine the clean reads of metagenomic data and the imported sequence library, and use a custom gene calculation script to quantify the three types of integrase genes to obtain abundance information.

[0063] 5) Screening and determination of primers for variable regions of gene cassettes of three types of integrons:

[0064] In step 5), the analysis method comprises the following steps:

[0065] 5-1) Search the literature database Web of Science with the keywords "gene cassette* or intI1 or intI2 or intI3 or integron*" and "antibiotic resistance gene*" to obtain relevant literature;

[0066] 5-2) Summarize candidate primers based on the base sequences and amplification conditions of three types of integron variable region primers used or recorded in the literature;

[0067] 6) Perform PCR amplification of class I, II, and III integron gene cassettes (GC1, GC2, and GC3) on the DNA, and perform gel electrophoresis and imaging to determine the presence and diversity of integron gene cassette fragments in the sample;

[0068] In step 6), the analysis method comprises the following steps:

[0069] 6-1), using the three types of primers screened to amplify the target integron gene cassette fragment for amplification;

[0070] 6-2) Perform pulsed field gel electrophoresis analysis on the amplified product;

[0071] 7) Purify and quality control the amplified product, and then perform sequencing on the PacBio machine; the data is split and quality controlled to obtain the Q30 sequence;

[0072] 8) Analyze the long read sequence to identify the drug resistance gene box carried by the integron in the sample and its arrangement array information.

[0073] In step 8), the analysis method comprises the following steps:

[0074] 8-1) The long-read sequencing raw data was subjected to barcode splitting, quality control and filtering based on RS_Readsofinsert.1protocol in SMRT Portal version 2.3.0 to obtain Q30 sequence;

[0075] 8-2) Use Prodigal (v2.6.3) to predict the coding gene ORF of the sequence after quality control;

[0076] 8-3) Using BLASTP (e-value ≤ 10 -10 ) were compared with the ARGs database SARG to annotate drug resistance genes (the screening criteria were similarity ≥ 80% and coverage ≥ 75%), thereby identifying sequences carrying ARGs (ACCs).

[0077] 8-4) Based on the ACCs sequence ID, the complete sequence was extracted using seqkit (v2.1.0) software and compared with the localized NR database (e-value ≤ 10) using Diamond software (v2.0.15). -2 ) to annotate all genes on ACCs.

[0078] 8-5) Use Python self-built scripts to extract sequences carrying ARGs, obtain the annotation results of each gene on the sequence, predicted sites and gene directions, etc., draw the above information into structured text, use the self-built script based on the R language gggenes package to visualize the integron-resistance gene box, and build an automated analysis process.

[0079] 9) By mapping and aligning the metagenomic sequencing data with the long-read sequencing results, the abundance information of GC1, GC2, and GC3 in the sample was obtained, and the variable region GC with higher abundance carrying multidrug resistance genes was visualized.

[0080] In step 9), the analysis method comprises the following steps:

[0081] 9-1) Use the contig mode of the CoverM software, input the integron gene box GC1, GC2, GC3 sequences and the sample metagenomic quality control short sequence, and select the "rpkm" in the parameter "-m" to calculate the abundance of the integron variable region fragment.

[0082] 9-2) Combined with the gene arrangement results in step 8), corresponding to the variable region sequence ID, the gene arrangement and abundance visualization were performed based on the self-built script of the R language gggenes package.

[0083] Example

[0084] Analysis of relative abundance of integrons and ARG-carrying gene cassettes in the chlorination disinfection system of a sewage treatment plant.

[0085] 1. The steps for obtaining the GC amplification products of the DNA of microbial samples before and after chlorination disinfection and their third-generation sequencing Q30 data are as follows:

[0086] 1) Enrichment of total microbial DNA

[0087] The data of a certain city sewage treatment plant (138,000 m 3 / day) influent samples were used as composite microbial flora systems for chlorination disinfection process. A total of 20L of water samples were collected, transported to the laboratory within 2 hours, and stored in a 4℃ refrigerator. The chlorination disinfection treatment group system set up six groups of treatment groups with chlorination disinfection concentration ranges commonly used in water treatment processes (1–32mg / L) and a control group without chlorination disinfection (Table 1). The reaction time was 30min, and three parallel treatments were set up in each group. The chlorination disinfection system was 500mL. Immediately after the chlorination and chlorine termination reactions were completed, 450mL of water samples were taken from 21 systems before and after chlorination treatment, and filtered through a 0.22μm filter membrane using a circulating water vacuum pump (SHZ-S(Ⅲ)ABS, Lichen Technology) for DNA extraction; the remaining water samples were added to 50mL centrifuge tubes as materials for microbial culture.

[0088] Table 1 Simulated chlorination disinfection reaction system settings

[0089]

[0090] 2) DNA extraction and quality control

[0091] Total DNA of 21 samples was extracted using a DNA kit (DNeasy Power Soil Pro Kit (QIAGEN, Germany)), and the extracted DNA concentration was detected using the Qubit 4.0 Nucleic Acid / Protein Quantitative Fluorometer and 1dsDNA HS (High Sensitivity) Analysis Kit (Thermo Fisher Scientific, USA). The DNA purity was determined according to the protocol using a NanoDrop Spectrophotometer (Thermo Fisher Scientific, USA). In addition, 1% agarose gel electrophoresis (voltage 5 V / cm, 20 min) was used to evaluate the quality of the DNA, and the markers used were λ-Hind III digest (TaKaRa) and DL5,000 DNA Marker (TaKaRa). The extracted DNA was stored at -80°C.

[0092] 3) Amplification of gene cassette and third-generation sequencing

[0093] In the technical solution of the present invention, the amplification primers and conditions for the three types of gene cassette PCR required in step 6-1 are shown in Table 2:

[0094] Table 2 Amplification primers and conditions for PCR of three types of gene cassettes

[0095]

[0096] In the technical solution of the present invention, the PCR amplification system in step 6-1 is set as follows:

[0097] The GC amplification system of PCR was set to 50 μL, where:

[0098] a) Ex Taq TM Version 2.0 Taq enzyme premix (TaKaRa): 25 μL;

[0099] b) Primers for class I, class II, and class III integron gene cassettes (concentration 10 μM): forward primer 1 μL; reverse primer 1 μL;

[0100] c) Sterile deionized water: 22 μL;

[0101] d) Template DNA: 1 μL.

[0102] Leveraging PacBio The IIe sequencing system sequenced the three types of GC amplification products, and the data were split and quality controlled to obtain the Q30 sequence. Figure 1 The gel electrophoresis images of the PCR products of the class I (a), II (b), and III (c) integron variable region gene cassettes of the total microbes (Total microbes, A) and lactose-fermenting Enterobacteriaceae (LFE strains, B) in the sewage microbiome after chlorination disinfection. Figure 2 a is the average length distribution of the sequence of the variable region gene cassettes of class I, II, and III integrons of total microorganisms disinfected by chlorination, and b is the average length distribution of the sequence of the variable region gene cassettes of class I, II, and III integrons of lactose-fermenting Enterobacteriaceae.

[0103] 2. Calculation of abundance of three types of integrons and annotation, abundance calculation and visualization of drug resistance genes of Q30 sequences in variable regions. The steps are as follows:

[0104] 1) Calculation of abundance of three types of integrons:

[0105] Use a self-built Python script to align the sample metagenome clean reads with the database constructed in steps 1-4 to obtain the samples that meet the criteria of “similarity ≥ 90%, coverage ≥ 50%, and e-value ≤ 10 -5" was identified as class I, II, and III integrase genes, and their corresponding abundances were output. The abundances of the three types of integrase genes are equivalent to the abundances of the three types of integrons in environmental samples, and the unit is copy number / cell (copy / cell).

[0106] 2) Drug resistance gene annotation:

[0107] Prodigal (v2.6.3) software was used to predict the ORF of the gene cassette, and then BLASTP (screening criteria: e-value ≤ 10 -10 ) were compared with the ARGs database SARG to annotate drug-resistant genes (screening criteria were similarity ≥ 80% and coverage ≥ 75%), thereby identifying sequences carrying ARGs (ACCs).

[0108] 3) Other functional gene annotations:

[0109] Based on the ACCs sequence ID, the complete sequence was extracted using seqkit (v2.1.0) software, the ORF of the gene cassette was predicted using Prodigal (v2.6.3) software, and the localization NR database (e-value ≤ 10) was aligned using Diamond software (v2.0.15). -2 ) to annotate other functional genes on ACCs.

[0110] 4) Calculation of abundance of three types of integron variable regions:

[0111] The abundance of all sequences was calculated using the contig mode of CoverM (v0.6.1) software. Specifically, the metagenomic clean reads of the microbial samples were mapped to the integron gene cassettes GC1, GC2, and GC3 sequences, and the abundance calculation method parameter "-m" selected "rpkm".

[0112] 5) Gene arrangement representation and visualization:

[0113] The arrangement of the integron resistance gene cassettes was visualized using a visualization script based on the R package “gggenes”. Figure 3 The arrangement and abundance characteristics of the super integrons carrying 6 ARGs with high risk in the chlorination disinfection system are shown; in the sample ID, "T" represents total microorganisms, "C" represents the sample before chlorination disinfection, "1" represents 1 mg / L chlorination disinfection exposure treatment, and "2" represents 2 mg / L chlorination disinfection treatment. "T1" means chlorination disinfection treatment with 1 mg / L total microorganisms in the system.

[0114] In summary, the present invention is based on second-generation metagenomic sequencing combined with PCR amplification and third-generation long sequence sequencing technology, creating a quantitative analysis of the abundance of three types of integrons in complex microbial systems in diverse environmental samples and characterization of the (multiple) drug-resistant gene cassettes they carry. It provides a new perspective for studying integron-mediated mobile high-risk drug-resistant genes and their evolutionary mechanisms, which not only broadens the cognitive system of bacterial resistance pollution behavior in the environment, but also provides an important scientific basis for the targeted control of high-risk drug-resistant pollutants.

Claims

1. A method for identifying, quantifying and visualizing the genome of an integron-resistance gene cassette in an environmental sample, characterized in that: The following steps are involved: Step 1: Construct the integrase gene database of three types of integrons, including: 1.1: Collect all integrase gene sequences; 1.2: Manually screen the fully matched sequences based on the sequence description information, and remove the sequences that are not integrase or lack the integrase genes intI1, intI2, and intI3 in the description content; 1.3: Remove redundant sequences and verify sequence accuracy, and then construct an integrase gene nucleotide sequence seed library; 1.4: Use Diamond tool to compare the seed library sequence with the non-redundant protein library of NCBI database, and set the parameter to e value ≤ 10 -5 , based on the description, non-integrase sequences were screened and eliminated, and then the sequences were de-redundant and verified for accuracy, and then an integrase gene amino acid sequence database was constructed; Step 2: Extract microbial DNA from environmental samples for second-generation metagenomic sequencing and perform quality control to obtain a metagenomic dataset; Step 3: Identification and quantification of integrons: Using the metagenomic dataset obtained in step 2, based on the integrase gene amino acid sequence database constructed in step 1, a custom Python script was imported to perform identification and quantitative analysis of the three types of integrase genes I, II, and III in environmental samples. The relative abundance of integrase genes is equivalent to the relative abundance of integrons in environmental samples. Step 4: Screen three types of candidate primers for integron gene cassettes based on literature, use candidate primers combined with PCR technology to amplify integron gene cassettes in environmental sample DNA, use pulsed field gel electrophoresis to visualize gene cassette sequence diversity and length, and obtain three pairs of primers and their amplification products with the highest diversity of gene cassette amplification products through comparison and selection; Step 5: The amplified product obtained in step 4 is subjected to third-generation long-fragment sequencing and processed to obtain a Q30 sequence; Step 6: Use the Q30 sequence obtained in step 5 to identify drug-resistant genes, use Prodigal to predict open reading frames, compare the ORF sequence with the structured drug-resistant gene database, and screen and identify drug-resistant genes according to the standard of "similarity ≥ 80% and coverage ≥ 75%", and then obtain the gene cassette sequences of class I, II, and III integrons carrying drug-resistant genes; use the Diamond tool to compare the ORF sequences of the gene cassettes of class I, II, and III integrons carrying drug-resistant genes with the comprehensive non-redundant database NCBI NR library to annotate all functional genes in the sequence, and set the parameter to e value ≤ 10 -2 ; Organize the information of drug-resistant genes and other functional genes carried by integrons and summarize them into structured text data. Use a self-built script based on the R language gggenes package to visualize the diversity arrangement of integrons-antibiotic resistance gene boxes; Step 7: Based on the metagenomic dataset obtained in step 2 and the gene cassette sequences of class I, II, and III integrons carrying drug-resistant genes obtained in step 6, the metagenomic sequences were mapped to the gene cassette sequences of class I, II, and III integrons carrying drug-resistant genes using the contig mode of CoverM software to characterize the relative abundance of gene cassettes of class I, II, and III integrons carrying drug-resistant genes in environmental samples; wherein: The step 3 specifically includes: firstly, using the makedb command to format the three types of integrase gene databases constructed in step 1, using the BLASTP tool to compare the environmental sample metagenomic data set obtained in step 2 with the database, and annotating the results that meet the "similarity ≥ 90%, coverage ≥ 50% and e-value ≤ 10 -5 " were identified as class I, II, and III integrase genes; the environmental sample metagenomic dataset was compared with the 16S database Greengenes using BWA and BLASTN tools to calculate the number of 16S genes in each sample; Diamond was used to compare the bacterial and archaeal single-copy marker gene databases to calculate the number of cells in each sample; the following formulas were used to quantify the abundance of the three types of integrase genes in the samples: Where: A is the abundance of class I, II or III integrase genes, which is equivalent to the abundance of integrons in environmental samples, and the unit is copy number / cell; N i(intI-likesequence) is the number of sequences annotated to the intI1, intI2 or intI3 reference sequences; L reads is the read length of the metagenomic sequence; L i(intI referencesequence) is the length of the intI1, intI2 or intI3 reference sequence; N 16Ssequence is the number of 16S rRNA genes identified by comparing metagenomic sequencing data with the Greengenes database; L 16Ssequence is the average length of 16S rRNA gene in Greengenes database, where L 16Ssequence The value is 1432bp; n is the number of intI1, intI2 and intI3 reference sequences mapped; N 16Scopy number is the average 16S rRNA gene copy number in cells, calculated using CopyRighter software.

2. The method according to claim 1, characterized in that: The step 1.1 specifically includes: selecting the "gene" sub-library of the NCBI database and searching and obtaining the corresponding gene sequences with the keywords of class I, II, and III integrase genes "intI1, intI2, intI3, class 1 integrase, class 2 integrase, and class 3 integrase"; and extracting the integrase sequence from the SEARCH / BLAST module of the integrin database INTEGRALL by keyword search.

3. The method according to claim 1, characterized in that: The step 1.3 specifically includes: clustering the sequences obtained in step 1.2 according to 100% sequence similarity using CD-HIT software to obtain non-redundant integrase nucleotide sequences, and then using the Pfam protein family domain database and CDsearch tool for comparison based on the principle of conserved domain similarity to verify sequence accuracy, thereby forming an integrase gene nucleotide sequence seed library.

4. The method according to claim 1, characterized in that: The step 4 specifically comprises: searching for literature in the literature database Web of Science using the keywords "gene cassette* or intI1 or intI2 or intI3 or integron*" and "antibiotic resistance gene*", and summarizing candidate primers according to the base sequences and amplification conditions of the three types of integron gene cassette primers used or recorded in the literature.

5. The method according to claim 1, characterized in that: The step 7 specifically includes: outputting the relative abundance of gene cassettes of class I, II, and III integrons carrying drug resistance genes in environmental samples under the contig mode of CoverM software, and normalizing them to RPKM values, i.e., reads per kilobase per million mappings, and the calculation formula is as follows: Where: Exon Mapped Reads is the number of reads mapped to the target drug resistance gene cassette; Total Mapped Reads is the total number of mapped reads; Exon Length is the length of the target drug resistance gene cassette.

Images