key: cord-0868723-fhjmygpu
authors: Lin, Yufei; Dong, Xiaohong; Sun, Rui; Wu, Jiao; Tian, Lejin; Rao, Dawei; Zhang, Lihua; Yang, Kun
title: Migratory birds-one major source of environmental antibiotic resistance around Qinghai Lake, China
date: 2020-05-30
journal: Sci Total Environ
DOI: 10.1016/j.scitotenv.2020.139758
sha: 88ada397ec9be5f9454b8f8307d00282e8230f15
doc_id: 868723
cord_uid: fhjmygpu

Abstract Migratory birds are potential transmitters of bacterial antibiotic resistance. However, their role in the environmental dissemination of bacterial antibiotic resistance and the extent of their impact on the environment are not yet clear. Qinghai Lake is one of the most important breeding and stopover ground for the migratory birds along the Central Asian Flyway. Here, we investigated the bacterial antibiotic resistance in the environment and among the migratory birds around the lake. The results of culture-based analysis of bacterial antibiotic resistance, quantitative PCR and metagenomic sequencing indicate that migratory birds are one major source of bacterial antibiotic resistance in the environment around Qinghai Lake. Network analysis reveals the co-occurrence patterns of antibiotic resistance genes (ARGs) and bacterial genera. Genetic co-localization analysis suggests high co-selection potential (with incidence of 35.8%) among different ARGs, but limited linkage (with incidence of only 3.7%) between ARGs and biocide/metal resistance genes (BMRGs). The high genetic linkage between ARGs and mobile genetic elements (MGEs) is still largely confined to the bacterial community in migratory birds (accounting for 96.0% of sequencing reads of MGE-linked ARGs), which indicates limited horizontal transfer of ARGs to the environment. Nevertheless, the antibiotic resistance determinants carried by migratory birds and their specific genetic properties (high co-selection and mobility potential of the ARGs) remind us that the role of migratory birds in the environmental dissemination of bacterial antibiotic resistance deserves more attention.

J o u r n a l P r e -p r o o f 3 The number of migratory birds worldwide is as high as 5 billion a year (Guenther et al., 2012) . Their migration flyways cover all continents including the Antarctica. The role of migratory birds in facilitating the transfer of antibiotic-resistant bacteria and resistance genes has been gradually recognized (Guenther et al., 2012; Atterby et al., 2016; Ramey et al., 2018; Marcelino et al., 2019) . The potential of intercontinental transportation of antibiotic-resistant bacteria and resistance genes via avian migration was also proposed (Liakopoulos et al., 2016; Mohsin et al., 2016; Lin et al., 2020) .

Qinghai Lake is the largest salt-water lake in China and one of the most important breeding and stopover grounds for many migratory birds along the Central Asian Flyway. Qinghai Lake is in the northeastern part of Qinghai province. The lake is less affected by human activities due to the official protection and National Nature Reserve policies. However, more than 150,000 migratory birds use this lake as breeding or stopover ground annually (Zheng et al., 2018) . In our previous work, we reported the transportation of bacterial antibiotic resistance from a polluted urban river (Jin River, Chengdu, China) to surrounding environment (Wangjianglou Park) mediated by the wild birds (the egrets) (Wu et al., 2018; Lin et al., 2020) . In the similar manner, the migratory birds may transport antibiotic-resistant bacteria or resistance genes from other places along the migration route to the environment around Qinghai Lake. The antibiotic-resistant bacteria and resistance genes, as well as their potential transportation between the environment and the migratory birds under such a unique geographical background were seldom reported.

J o u r n a l P r e -p r o o f 4 ( Supplementary Fig. S1 ). The aim of this study is to explore the possibility and extent of bacterial antibiotic resistance transmission through migratory birds in the wild environment such as Qinghai Lake watershed. The genetic linkages among different antibiotic resistance genes (ARGs), between ARGs and biocide/metal resistance genes (BMRGs), and between resistance genes (ARGs or BMRGs) and mobile genetic elements (MGEs) were evaluated.

Sampling efforts were carried out at 8 sites around Qinghai Lake in June 2017 (Supplementary Fig. S2 and Fig. S3) . It was in the early middle stage of the breeding season by when most of the migratory birds using Qinghai Lake as breeding ground had been back in the place. The sampling specimen include water, soil/sands and bird feces. Water samples were collected in sterile 1-L bottles. Beach soil/sands were collected from the foreshore using a trowel. Multiple soil/sand samples were collected within a half-meter-radius circle, pooled, and stored in the brand-new self-sealing bags. Fresh droppings of migratory birds were collected in autoclaved 50-mL Falcon tubes. Detailed information for the sampling efforts is summarized in Table S1 . In total, we collected 26 bird feces, 18 water samples and 24 sand/soil samples. Among the 26 bird feces samples, 18 were from bar-headed geese (Anser indicus, AI), 2 were from great cormorants (Phalacrocorax carbo, PC), 2 were from J o u r n a l P r e -p r o o f 6 inoculate the testing 96-well plates containing antibiotics of different concentrations.

A flame-sterilized 48-pin replicator was used to accelerate the inoculating operation.

The MIC endpoints were determined as the lowest concentration with no visible growth after 20 h of incubation at 37 o C. Duplicate tests were conducted for each antibiotic concentration. Positive and negative controls were conducted in antibiotic-free media to assess the growth of E. coli isolates and sterility of the assay, respectively. Quality control of the procedure was conducted by using the susceptive E. coli ATCC 25922.

The multiple-antibiotic resistance index (MARI) was calculated for each type of sample (Krumperman, 1983; Wu et al., 2018) . The antibiotic resistance patterns of the E. coli isolates from different samples were compared via non-metric multidimensional scaling (NMDS). The MICs of E. coli isolates were log-transformed, z-scaled, and used as input data for NMDS. The NMDS was performed in Matlab_2016Ra (MathWorks, Natick, MA, United States) with the mdscale function (Wu et al., 2018) .

Genomic DNA was extracted from sand/soil (1g each), bird feces (around 0.5g each) and filtration retentate of water (from about 300 mL water sample each filtered through GN-6 0.45 μm membranes) with the FastDNA ® Spin Kit for Soil (MP J o u r n a l P r e -p r o o f 7 of ARGs and MGEs via qPCR. The 16S rRNA gene was used as reference standard.

The abundances of ARGs and MGEs were expressed as the amount ratios of target genes to the reference gene (Wu et al., 2018) . Twenty-one genes were chosen for the quantification, which included typical ARGs conferring resistance to the 11 antibiotics we tested in this work, as well as some important MGE sequences (Supplementary Table S3 ). Each qPCR was operated in triplicate. Calculated abundances of target genes in different samples were visualized in a heatmap. The relative abundances of target genes were log-transformed, z-scaled, and used as input data for heatmap construction in Matlab_2016Ra (MathWorks, Natick, MA, United States).

DNA extracts of the same sample type were pooled together to make 9 composite DNA samples representing fresh water (FW), salt water (SW), sand/soil (SS), sheep feces (SH), yak feces (BG), bar-headed goose feces A (AIA), bar-headed goose feces B (AIB), great cormorant feces (PC) and great black-headed gull feces (LI), respectively. The DNA sample extracted from the ruddy shelduck feces (TF) did not pass the quality test before sequencing. Therefore, we abandoned to sequence this sample. Metagenomic libraries were generated from 1μg DNA per sample using We obtained an average of around 10.5 Gb (giga bases) of metagenomic data for each DNA sample, resulting in a total of 94.3 Gb sequencing data. The raw reads were quality filtered via Sickle (https://github.com/najoshi/sickle). We used SeqPrep (https://github.com/jstjohn/SeqPrep) to merge paired-end Illumina reads that are overlapping into a single longer read. Clean reads from each sample were assembled with SOAPdenovo software (http://soap.genomics.org.cn/, Version 1.06) (Li et al., 2008) . Assembled contigs longer than 500 bp were used for the prediction of open reading frames (ORFs) utilizing MetaGene (http://metagene.cb.k.u-tokyo.ac.jp/) (Noguchi et al., 2006) , and a non-redundant ORF set was obtained by using CD-HIT software (http://www.bioinformatics.org/cd-hit/) with thresholds of identity and coverage at 95% and 90%, respectively (Li and Godzik, 2006) . Finally, clean reads were mapped back (with an identity greater than 95%) to the non-redundant ORF set using SOAPaligner (http://soap.genomics.org.cn/) to calculate the gene abundance for each sample (Li et al., 2008) .

J o u r n a l P r e -p r o o f 9 community.

In order to analyze the abundance of ARGs, BMRGs and MGEs in the sequencing samples, we downloaded several related databases from relevant websites and constructed local databases for BLAST searching. The non-redundant ORF set was searched for ARGs against the Comprehensive Antibiotic Resistance Database (CARD, https://card.mcmaster.ca) using BLASTP with E-value ≤ 1 × 10 −5 (McArthur et al., 2013) . A query sequence was annotated as an ARG-like fragment if it had ≥80% amino acid identity with the subject sequence in the CARD and the alignment length ≥25 amino acids. To investigate the BMRGs, the ORF set was compared against the BacMet database (http://bacmet.biomedicine.gu.se) using BLASTP with E-value ≤1 × 10 −5 (Pal et al., 2014) . Similarly, the amino acid identity and alignment length thresholds were set as 80% and 25 amino acids, respectively. To determine the MGEs, the ORF set was matched against the databases of ACLAME (plasmids, pro-phages and viruses, http://aclame.ulb.ac.be) (Leplae et al., 2004) , INTEGRALL (integrons, http://integrall.bio.ua.pt) (Moura et al., 2009) and ISfinder (inserting sequences, https://www-is.biotoul.fr) with the similar method. When encountering multiple hits in the database, we just kept the subject sequence with the highest identity (score).

The abundance of an annotated ARG, BMRG or gene related to a MGE in a DNA sample was calculated according to the following equation (Li et al., 2015) . N 16S sequence is the reads number of the 16S sequence identified from the same metagenome sample. L 16S sequence is the average length of 16S sequences, which was calculated to be 1445 bp from sequences of Greengenes database.

In order to evaluate the horizontal transfer potential of resistance genes (ARGs or BMRGs), we introduced the concept of mobility incidence (Ju et al., 2019) , which calculated the probability of co-localization of a resistance gene with MGEs (mobility incidence of the resistance gene). However, in our case, the definition of mobility incidence is slightly different. For individual metagenome sample, the mobility incidence (M%) of an ARG/BMRG is defined as the ratio of detection rate (reads number) of the ARG/BMRG co-localized with at least one MGE on the same contig/scaffold against the total detection rate (reads number) of all ARGs/BMRGs from the same metagenome sample. For an individual ARG/BMRG, the mobility incidence (M%) is the ratio of detection rate of the ARG/BMRG co-localized with at least one MGE on the same contig/scaffold against the total detection rate (reads number) of the same ARG/BMRG among all metagenome samples. To evaluate the co-selection potential between different ARGs or between ARGs and BMRGs, we J o u r n a l P r e -p r o o f 11 community. This would explore the preferred gene types of co-selection events in the environment. A correlation matrix was constructed by calculating all possible pairwise Spearman's rank correlations (ρ) among ARGs and bacterial genera (Li et al., 2015) .

Only connections with a strong (Spearman's ρ > 0.8) and significant (P value < 0.01) correlation were retained for network analysis. Correlation coefficient and P value were calculated and filtered with RStudio using package Psych 

Totally, 723 E. coli strains were isolated from the collected samples (Supplementary   Table S2 ). 75 isolates were recovered from soil/sand (SS), 93 from fresh water (FW), 60 from the ruddy shelduck feces (TF), 60 from great cormorant feces (PC), 60 from great black-headed gull feces (LI), 80 from bar-headed goose feces of type A (AIA), 115 from bar-headed goose feces of type B (AIB), 90 from sheep feces (SH) and 90 from yak feces (BG). The bacterial antibiotic resistance patterns and levels of different samples are exhibited in Fig. 1 water also gives high bacterial antibiotic resistance (MARI=0.140). The bacterial antibiotic resistance in local soil/sand (FW, MARI=0.052) and ruminants (sheep and yaks, MARI=0.004 and 0.002, respectively) are low. We did not recover any E. coli from salt-water samples except for three isolates from W8-1. Most of the soil/sand (SS) E. coli (96%) were isolated from the soil/sand samples collected at sampling site 3. Most of the fresh water (FW) E. coli isolates (75%) were recovered from the freshwater samples collected at sampling sites 3 (n=36) and 4 (n=34).

The relative abundance of ARGs and MGE sequences in environmental samples determined via qPCR is exhibited as a clustered heatmap (Fig. 2 ). Among the 21 tested genes, 12 are detected the highest relative abundance in bird feces samples, 5 are detected the second highest relative abundance in bird feces samples. Some ARGs conferring resistance against the same type of antibiotics are clustered together, such as the resistance genes against tetracycline (tetW, tetO and tetL), the resistance genes against fluoroquinolone (qnrA and qnrS) and the resistance genes against beta-lactamase genes (bla NMD-1 and bla CTX-M-14 ) etc.

Taxonomic annotation of the metagenomic sequencing data against the NR database J o u r n a l P r e -p r o o f 13 indicates that the bacteria are the dominant domain, accounting for more than 96% of the DNA sequencing reads. Among them ， Proteobacteria, Bacteroidetes and Firmicutes are the major constituents, averagely accounting for 57.9%, 15.9% and 7.6% of the bacterial community ( Supplementary Fig. S4 ) respectively, followed by Actinobacteria (3.6%), Verrucomicrobia (3.3%) and Planctomycetes (2.6%). Fig. S5 ).

The samples AIA and AIB shared high similarity in both bacterial community J o u r n a l P r e -p r o o f 14 structure and metabolic module composition. Accordingly, we believe that these two types of samples can be regarded as the same one, although the appearance of the two samples is somehow different.

According to the criteria of the BLAST searching (identity >80%, E-value ≤1 × 10 −5 ), 115 major ARGs (with relative abundance >10 -4 in at least one sample) were annotated among 9 metagenome samples. Most of these ARGs are detected in the feces of migratory birds (n=105) and local ruminants (n=61) (Fig. 3A) , while ruminants and migratory birds share 53 ARGs. The ARGs detected in the feces of migratory birds are mainly from bar-headed geese (AIA and AIB) and great black-headed gulls (LI) (Fig. 3B) . Based on the annotation results of metagenomic data, we calculated the abundances of ARGs, BMRGs and MGEs in all metagenome samples. Feces samples of migratory birds show relatively higher ARG abundance ( Fig. 3C and Supplementary Dataset S2). Bar-headed geese have the highest ARG abundance, followed by great black-headed gulls (LI) and great cormorants (PC).

Although ruminants share a large portion of ARGs with migratory birds (n=53), the abundance of these ARGs in ruminants is much lower than that in migratory birds (Fig. 3D) . The lowest abundance of ARGs occurs in soil/sand ( We noticed that the freshwater E. coli isolates were mostly recovered from the water samples collected at sites 3 and 4. Site 3 is the active location of migratory birds, while site 4 is the estuary of Heima River. The estuary of Heima River is near the town of Heima River where the river may suffer from anthropogenic impacts. When we performed the metagenomic sequencing, we pooled the DNA samples from W1-1, W1-2, W3-3, W4-1, W4-2, and W5-1 together as a composite sample of FW. We also recovered E. coli isolates from freshwater samples W1-1 (n=13) and W5-1 (n=7).

Only one of the 20 isolates showed multi-drug resistance. Therefore, the high rate of antibiotic resistance among the fresh-water isolates is believed to be the result of sampling bias.

BMRGs are a class of stress resistance genes that are most likely to form co-selection effects with ARGs. Due to the situations called "cross-resistance" (Seiler and J o u r n a l P r e -p r o o f 16 the abundance of ARGs (Fig. S6) . The highest abundance is detected in AIB, followed by AIA and LI (Supplementary Dataset S4). The lowest abundance of BMRGs also occurred in soil/sand (SS). Most metal resistance genes are of the type Multimetal (n=33). Thirty-seven kinds of biocide resistance genes are detected among all metagenome samples. The resistance genes against specific metals such as copper (Cu, n=9), arsenic (As, n=7), zinc (Zn, n=5), mercury (Hg, n=3), nickel (Ni, n=3), aluminum (Al, n=2), chromium (Cr, n=2) and tellurium (Te, n=1) are also detected (Supplementary Fig. S7 and Dataset S5). Only two cases are confirmed for the co-localization of BMRGs and ARGs on the same contig/scaffold (Supplementary Table S4 ).

The abundance of MGEs among the metagenome samples is out of our expectation.

The highest abundances of plasmid, IS and integron all occur in the sample of SW ( Supplementary Fig. S8 ). Some cases for co-localization of ARGs and MGE sequences on the same contig/scaffold are listed in Supplementary Table S5 . Different from the abundance distribution of MGEs, the co-localization of ARGs and MGEs was mainly detected in the migratory birds and ruminants. According to the calculated mobility incidences (M%), the horizontal transfer potential of various resistance genes in migratory birds is the highest, followed by that in the local ruminants ( Supplementary Fig. S9A and B) . The lowest horizontal transfer potential of resistance genes occur in those local environment samples (salt water, fresh water and J o u r n a l P r e -p r o o f 17 soil/sand samples). The horizontal transfer potential of each resistance gene type is shown in Supplementary Fig. S9C and D (C for ARGs and D for BMRGs).

The co-occurrence network derived from the strong (ρ>0.8) and significant (P<0.01) correlations among bacterial genera and ARGs are exhibited in Fig. 4 . The network consists of 106 nodes and 612 edges. ARGs of the same type exhibit high co-occurrence rate ( Figure 4A ). Some bacterial genera (such as Clostridium) are closely related to ARGs (Fig. 4A) . The correlations between bacterial genera and ARGs are summarized in Table S6 . Based on the algorithm of fast unfolding of communities in large network, the entire network can be clustered into 9 major modules with modularity index of 0.478 (Blondel et al., 2008) . Three main modules occupied 76 out of 105 nodes and 548 out of 612 edges (Fig. 4B) . The most densely connected node is defined as the hub of a module, which can act as the indicator of the module (Li et al., 2015) . Its abundance reflects the abundances of other co-occurring ARGs in the same module ( Supplementary Fig. S10 ). The hub and co-occurring ARGs of the three main modules are listed in Supplementary Table S7. J o u r n a l P r e -p r o o f 18 Notably, there are some inconsistencies among the drug-resistance phenotype of E.

coli isolates, qPCR quantitative results of resistance genes and metagenomic sequencing data. Three possible reasons may explain these inconsistencies. First, we only recovered E. coli isolates from a few samples and analyzed their resistance phenotype. The antibiotic resistance of specific bacteria may not fully reflect the overall resistance of the whole bacterial community. In addition, the bacterial antibiotic resistance in individual samples (such as the freshwater samples of site 3) is not universal, which may not reflect the common characteristics of the same kind of samples. Therefore, great care should be taken when interpreting the phenotypic data of bacterial antibiotic resistance, especially for the fecal samples of migratory birds with small number, such as the fecal sample of ruddy shelduck, from which the antibiotic resistance of Escherichia coli isolated is weak. It cannot rule out the possibility of individuals with high bacterial antibiotic resistance. Second, the variety of resistance genes we selected for qPCR quantification is very limited, which may not fully reflect the level of bacterial antibiotic resistance in each sample. Third, although the method of metagenomic sequencing can obtain the most comprehensive information about the bacterial antibiotic resistance in the sample, we might ignore the possible differences between individual samples of the same type due to the fact that we mixed the DNA extracts of the same sample type for metagenomic sequencing.

Furthermore, the analysis of the metagenomic data depends on the relevant database, in our case the CARD. As we know, the database includes a lot of nontransferable nonspecific drug-efflux genes that may be just involved in some detoxification process (Martinez et al., 2009 ). However, based on the results from this study, we can still draw the conclusion that the migratory birds are one of the main sources for the bacterial antibiotic resistance in the environment of Qinghai Lake. The abundance of J o u r n a l P r e -p r o o f 19 ARGs or movable ARGs is still the highest in migratory birds without considering those drug-efflux genes ( Supplementary Fig. S11 ).

Most of ARGs (58 out of 61) detected in ruminants are also detected in migratory birds (Fig. 3A) , and the relative abundance of these shared ARGs is mostly higher in migratory birds (Fig. 3D) . The shared ARGs between migratory birds and local ruminants (Fig. 4A ) may be due to their sharing the same fresh water resources ( Supplementary Fig. S3G ), while the disparity in the abundance of ARGs among the species may reflect the transmission direction of ARGs: mostly from migratory birds to ruminants. The types and abundance of ARGs in different species of migratory birds should be related to their feeding habits and migration routes. Bar-headed geese and ruddy shelducks are herbivorous. Great cormorants prey on fish. While greater black-headed gulls utilize multiple food resources. As regarding the migration route, bar-headed geese are considered one of the highest-flying birds. Their migration includes crossing a broad front of the Himalayas and covers the wide area between India and Mongolia along the Central Asian Flyway (Takekawa et al., 2009) . Great black-headed gulls and ruddy shelducks breeding in Qinghai Lake migrate to the coast of the Bay of Bengal to overwinter (Chu et al., 2008; Muzaffar et al., 2008; Takekawa et al., 2013) . Both India and Bangladesh have suffered heavily from the contamination of bacterial antibiotic resistance (Kakkar et al., 2017; Ahmed et al., 2019) . High levels of antibiotic resistance have been detected in wild birds (gulls, Chroicocephalus brunnicephalus) in these areas (Hasan et al., 2014) . This may explain why these migratory birds carry significantly higher level of bacterial antibiotic resistance than the local environment around Qinghai Lake. Great cormorants give low degree of bacterial antibiotic resistance and accordingly low abundance of ARGs, which may be because of their special migratory habit and J o u r n a l P r e -p r o o f 20 migration route. Unfortunately, no detailed information is available for the migration of the great cormorants around Qinghai Lake. We also observed various migratory birds gathering to share fresh water sources (Supplementary Fig. S3E and G) .

Therefore, antibiotic-resistant bacteria or resistance genes may be exchanged among migratory birds of various species. Our results confirmed that migratory birds of different species shared 10 common ARGs (Fig. 3B) . Usually, ARGs shared among different environmental compartments are believed to exhibit more mobility (with higher mobility incidences). However, according to the mobility analysis of ARGs in this work, it is not the case. Procrustes analysis demonstrates close correlations between the composition of shared bacterial community and the contents of shared ARGs among different environmental compartments ( Supplementary Fig. S12 ), which indicates that the transfer of ARGs in the environment of low selective pressure is mainly through the exchange of resistant bacteria.

In the co-occurring network analysis, the abundance of the hub gene of each module reflects the abundance of co-occurring ARGs of the same module, and indicates the origin of ARGs in the module ( Supplementary Fig. S10 ). We can derive that the ARGs of module IV mainly originate from bar-headed geese B (AIB, Supplementary   Fig. S10A ). The ARGs of module V mainly originate from local ruminants (sheep and yaks, Supplementary Fig. S10B ). The ARGs of module VI most likely originate from multiple migratory birds ( Supplementary Fig. S10C ): bar-headed geese (AIA, AIB), great black-headed gulls (LI) and great cormorants (PC). As we have incorporated the bacterial genera into the network analysis, it is interesting that the bacterial genera of the ruminant-originated module V include typical gut bacteria: Prevotella, J o u r n a l P r e -p r o o f especially fiber (Wu et al., 2011) . Fibrobacter is the major rumen bacteria (Montgomery et al., 1988) . The bacterial genera in the migratory-birds-originated module VI mostly are halophilic (Halomonas), osmotolerant and psychrophilic (Psychrobacter), which is related to the characteristics of local environment and climate. Some genera are typical marine or aquatic bacteria such as Vibrio, Marinomonas, Rheinheimera (Zhong et al., 2014) , Oceanisphaera and Shewanella (Satomi et al., 2003) . It is reasonable to detect marine bacteria in the feces-droppings of aquatic birds inhabiting around the saltwater lake. There are two possible reasons for the formation of modular structure in the co-occurrence network of ARGs. First, some ARGs are related to specific bacterial genera, which present in the bacterial communities in relation to specific environmental conditions (Li et al., 2015) .

Modules IV and V, which consist of many ARGs of the type of Drug efflux, are more likely to be in this situation. Some of the bacterial genera of these two modules are possible hosts of ARGs. More attention should be paid to the bacterial genus J o u r n a l P r e -p r o o f 22 against two or more antibiotics are adjacent to each other on one MGE (Seiler and Berendonk, 2012) , present in the bacterial communities that have experienced special multi-antibiotic selective pressure, may result in the enrichment of these ARGs.

Module VI containing ARGs against specific antibiotics such as beta-lactam, glycopeptide and phenicol etc. should be of this situation. The clustering of ARGs was often discovered in the environment under heavy anthropogenic impacts, via the method of quantitative PCR (Johnson et al., 2016) and metagenomic sequencing (Li et al., 2015; Chen et al., 2019) . We detected the co-localization of some ARGs with high co-occurrence rate (with incidence of 35.8%) on the same contigs/scaffolds (Table S8 ). Most of the co-localized ARGs are of the same module in the co-occurrence network. Quite some of these co-localized ARGs encode subunits of the same efflux pump. The co-localization of ARGs conferring resistance to different antibiotics is also observed, such as the colocalized ARGs conferring resistance to tetracycline (tet(G)) and aminoglycoside (aph(6)-Id). Therefore, the strong statistical correlation between resistant genes is determined by their actual physical linkage at the DNA molecular level. Co-selection will occur between genetically linked ARGs (Johnson et al., 2016) . There are only several cases (with low incidence of 3.7%) that metal resistance genes and ARGs are located on the same contigs/scaffolds 

of antibiotic resistomes in river environment: Application to an urban river in Beijing

Satellite tracking of the migration routes of the waterbirds breeding in Qinghai Lake

Clinical and Laboratory Standards Institute. Performance Standards for Antimicrobial Susceptibility Testing-Twenty-Ninth Edition: M100. CLSI

Resistance to antibiotics of clinical relevance in the fecal microbiota of Mexican wildlife

Spatial and temporal variation in enterococcal abundance and its relationship to the microbial community in Hawaii beach sand and water

Method 1103.1: Escherichia coli (E. coli) in water by membrane filtration using membrane-thermotolerant Escherichia coli agar (mTEC)

Comparable high rates of extended-spectrum-beta-lactamase-producing Escherichia coli in birds of prey from Germany and Mongolia

The gull (Chroicocephalus brunnicephalus) as an environmental bioindicator and reservoir for antibiotic resistance on the coastlines of the Bay of Bengal

Defining and combating antibiotic resistance from One Health and Global Health perspectives

Clusters of antibiotic resistance genes enriched together stay together in swine agriculture

Wastewater treatment plant resistomes are shaped by bacterial composition, genetic exchange, and upregulated expression in the effluent microbiomes

Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences

The colistin resistance mcr-1 gene is going wild

Metadata Analysis of mcr-1-Bearing Plasmids Inspired by the Sequencing Evidence for Horizontal Transfer of Antibiotic Resistance Genes Between Polluted River and Wild Birds

Highly pathogenic H5N1 influenza virus infection in migratory birds

Meta-transcriptomics reveals a diverse antibiotic resistance gene pool in avian microbiomes

Functional role of bacterial multidrug efflux pumps in microbial natural ecosystems

The comprehensive antibiotic resistance database

First description of plasmid-mediated colistin-resistant extended-spectrum beta-lactamase-producing Escherichia coli in a wild migratory bird from Asia

Transfer of Bacteroides succinogenes (Hungate) to Fibrobacter gen. nov. as Fibrobacter succinogenes comb. nov. and description of Fibrobacter intestinalis sp. nov

INTEGRALL: a database and search engine for integrons, integrases and gene cassettes

Seasonal movements and migration of Pallas's Gulls Larus ichthyaetus from Qinghai Lake

MetaGene: prokaryotic gene finding from environmental genome shotgun sequences

Avian influenza. Evidence points to migratory birds in H5N1 spread

Global patterns of influenza a virus in wild birds

BacMet: antibacterial biocide and metal resistance genes database

Antibiotic-resistant Escherichia coli in migratory birds inhabiting remote Alaska

Shewanella marinintestina sp. nov., Shewanella schlegeliana sp. nov. and Shewanella sairae sp. nov., novel eicosapentaenoic-acid-producing marine bacteria isolated from sea-animal intestines

Heavy metal driven co-selection of antibiotic resistance in soil and water bodies impacted by agriculture and aquaculture

Clostridium perfringens as a potential indicator for the presence of sewage solids in marine sediments

Geographic variation in Bar-headed Geese Anser indicus: Connectivity of wintering areas and breeding grounds across a broad front

Movements of wild ruddy shelducks in the Central Asian Flyway and their spatial relationship to outbreaks of highly pathogenic avian influenza H5N1

Antibiotic-induced shifts in the mouse gut microbiome and metabolome increase susceptibility to Clostridium difficile infection

Global trends in antimicrobial use in food animals

Metagenomic profiling of gut microbial communities in both wild and artificially reared Bar-headed goose

Linking long-term dietary patterns with gut microbial enterotypes

Evidence for environmental dissemination of antibiotic resistance mediated by wild birds

We gratefully thank the administration of Qinghai Lake National Natural Reserve for authorizing and assisting the collection of bird droppings samples. We thank the financial supports from the Natural Science Foundation of China (No. 21677104) andthe Initiating Research Fund for Talent Introduction of Sichuan University (No.