Whole genome sequence analysis of multi drug resistant community associated methicillin resistant Staphylococcus aureus from food fish: detection of clonal lineage ST 28 and its antimicrobial resistance and virulence genes

Fish samples collection

The MRSA strain7 isolated from Trichiurus lepturus (ribbon fish), a sea caught fish of Veraval, Gujarat, India. This sun dried fish is supplied to local markets and other parts of India mainly to North Eastern states and neighbouring countries.

Isolation, identification and antimicrobial susceptibility testing (AST) of S. aureus

ISO 6888-1 and ISO 6888-2:2003 (ISO, 2003) was followed for the isolation and identification of S. aureus from the fish samples. In brief, 10g of samples were homogenized in sterile Normal Saline and 0.3, 0.3 and 0.4 ml was poured on to Baird Parker Agar (BPA) plate after making the serial dilution. Characteristic colonies of S. aureus were purified and streaked onto Brain Heart Infusion (BHI) agar slant for further confirmation on the readymade MRSA plates BBL CHROMagar MRSA II (BBL, Difco).AST was carried out by disc diffusion method as per CLSI guidelines (CLSI, 2018) with 24 antibiotics (Dodeca Staphylococci-1 and 2, HiMedia, Mumbai) Mueller Hinton agar supplemented with 4% NaCl were used for the test.The inhibition zones were measured and to find the sensitivity and resistant. S. aureus ATCC 25923 and ATCC 43300 were used as reference strains in this study.

DNA isolation and identification MRSA by triplex PCR

DNA isolation was carried out using GenElute Bacterial Genomic DNA Kit (Sigma- Aldrich, Spain). The Multiplex PCR was carried out (Al-Talib et al., 2009) for simultaneous detection of S. aureus and MRSA with mec A (mecA- F: ACGAGTAGATGCTCAATATAA, mecA-R: CTTAGTTCTTTAGCGATTGC, 293 bp), fem A gene (femA-F: CGATCCATATTTACCATATCA, femA-R: ATCACGCTCTTCGTTTAGTT, 450 bp) and 16SrRNA primer specific to Staphylococcus genes specific detection (16S rRNA-F: GCAAGCGTTATCCGGATTT, 16S rRNA-R: CTTAATGATGGCAACTAAGC, 597 bp) and were Sanger sequenced for the confirmation. 50 µl of reaction mixture contained 200 µM dNTPs, 2.5mM MgCl₂, 1X PCR buffer, 0.5U Tag DNA polymerase (Sigma, Spain), 100 µg DNA/ µl with primer concentration of 0.6 pmol 16S rRNA, 0.8 pmolnuc, and 1.0 pmol mecA. The reaction was carried out using SureCycler 8000 (Agilent, Santa Clara, CA, USA) with initial denaturation at 94 °C (3 min) 34 cycles of denaturation at 94 °C (30 s), annealing at 60 °C (30 s), extension at 72 °C (30 s) and a final extension at 72 °C (5 min) and electrophoresis on 1.5% agarose gels with ethidium bromide in Gel Doc (BioRad, Hercules, CA, USA) (Fig. S1).

Whole genome sequencing and data analysis

The genomic DNA was isolated by bacterial genomic DNA isolation kit (Sigma-Aldrich, France), and its quality was checked on a NanoDrop. The whole-genome sequencing was carried out in IlluminaHiSeq 2500(paired end) (Illumina Inc., Cambridge, UK)and was assembled to determine the genetic structure and its multiple drug resistance.

The number of paired-end reads was approximately 7 billion short-read sequences in pairs of ∼300 bp, the number of bases (Mb) was 650.26, and there was 35.28% G+C content. De novo contig assembly was performed using MaSuRCA (Gutiérrez et al., 2012), and further downstream processing was performed. Coding sequences (CDSs) were predicted from the contigs using Glimmer (Zimin et al., 2013), and 2,693 predicted CDSs were found. Organism annotation, gene and protein annotation to the matched genes, gene ontology annotation, and pathway annotation were carried out with the use of the NCBI database. Overall, we observed that 2,669 (99.10%) of the predicted CDSs had at least one hit in the NCBI database. Among the total significant BLASTX hit CDSs, 1,703 genes were annotated using the UniProt database. The total number of Gene Ontology annotations identified for molecular functions was 889, with 604 annotations having to do with a biological process and 236 annotations having to do with cellular components. We predicted tRNA genes from the contigs using tRNAscan-SE (Delcher et al., 1999) and found 70 genes and the detailed analysis of the WGS are provided in the (Sivaraman et al., 2017).

Spa typing, Multilocus sequence typing (MLST), and prediction of virulence and antimicrobial genes

The spa typing, MLST analysis and VirulenceFinder2.0 of the MRSA was done by comparing the WGS of the MRSA using the Center for Genomic Epidemiology (http://www.genomicepidemiology.org/) (Larsen et al., 2012) where in raw sequencing reads are uploaded KMA is used for mapping Illumina sequencing platforms. The database includes following genes: hlb, hlgABC, tst, lukED, lukFS-PV, etAB, edinABC, aur, splABE, scn, sak, ACME and enterotoxins A-E, G-O, R, U, Q. The ARGs detection by The Comprehensive Antibiotic Resistance database (CARD) by Protein homolog criteria (Alcock et al., 2020).

Whole genome of CA- MRSA strain

The draft genome of a CA-MRSA strain7 contains 7,731,542 number of reads and 731,267,948 bases assembled into 120 contigs, 35.11% G + C, 80 tRNAs with the NCBI accession number NBZX00000000.1 (BioProject: PRJNA352109; BioSample: SAMN05967212), the N50 estimate was 66022 bp. It contained an average scaffold of 23759.70 bp, largest scaffold of 170965 bp and smallest of 301 bp. The predicted and annotation of the genome (NCBI Prokaryotic Genome Annotation Pipeline) shows 2797 coding sequences (CDSs) and the number of predicted CDSwith significant BLASTX match was 2755, number of predicted CDs with UniProt annotation was 1741. We observed that 100% of the predicted CDSs have a similarity of more than 60% at the 237 protein level with the existing proteins at the NCBI database. The majority of the top BLASTX hits 238 belonged to Staphylococcus species (top 15 organisms).The average gene size of 902, longest gene of 30258 and shortest gene of 72 bp in length (Tables S1 and S2). The gene ontology sequence distribution of MRSA7 strain contributed 51%, 35% and 14% towards molecular functions, biological process and cellular components, respectively (Fig. 1 and Supplementary Excel data file and Fig. S3). The molecular and biological process has almost equally made up the majority of the GO annotations, followed by cellular components. The subcategories like ‘metabolic process’ and ‘cellular process’ are in the cluster of biological process; two sub categories ’binding function’ and ‘catalytic function’ were clustered in Molecular Function; two sub categories of ‘cell’ and ‘cell part’ were clustered in Cellular Component and the details are provided in the Table S2. There is only 65 number of CA-MRSA Bioproject is available in the NCBI data (http://www.ncbi.nlm.nih.gov) and those were mainly from Healthcare Associated MRSA (HA-MRSA), whereas no reports from the fishery environment are available and it is the first kind of CA- MRSA in food fish of India.

Antimicrobial resistance pattern and ARGs by CARD Analysis of WGS

Higher level of resistance to rifampicin (76–100%), trimethoprim/sulfamethoxazole (69–81%), linezolid (69–75%), gentamicin (53–72%), piperacillin-tazobactam (53–72%), ampicillin/sulbactam (23–37%), cefoxitin (16%), and ofloxacin (12–23%) among these isolates from these fish samples were noticed (CLSI, 2018) and recognized as multidurg resistant strain (MDR) as they showed resistance to more than three classes of antibiotics. Complete susceptible to the following antibiotics such as azithromycin (AZM) 15 µg, clarithromycin (CLR) 15 µg, ciprofloxacin (CIP) 5 µg, gatifloxacin (GAT) 5 µg, clindamycin (CD) 2 µg, tigecycline (TGC) 15 µg, moxifloxacin (MO) 5 µg, lomefloxacin (LOM) 10 µg, norfloxacin (NX) 10 µg, novobiocin (NV) 30 µg, teicoplanin (TEI) 30 µg, nitrofurantoin (NIT) 300 µg, pristinomycin (RP) 15 µg were observed.

The CARD is a rigorously curated collection of known resistance determinants and associated antibiotics, organized by the Antibiotic Resistance Ontology (ARO) and Antimicrobial Resistance (AMR) detection models at McMaster University (McArthur et al., 2013; Jia et al., 2016). Results revealed that eight major ARO genes encoded for antimicrobial resistance such as mepR, mgrA, ar1R with perfect Resistance Gene Identifier (RGI) hit and S. aureus lmrS, fosB, norA and G1pT mutation with strict RGI hit. This strain possessed these resistance genes to 12 drug classes viz., aminoglycoside antibiotic, cephalosporin, diaminopyrimidine, fluoroquinolone, fosfomycin, glycylcycline, macrolide, oxazolidinone, peptide, phenicol and tetracycline antibiotic and also for an acridine dye resistance through mechanisms such as antibiotic efflux, antibiotic inactivation and antibiotic target alteration. ARGs showed 100% sequence similarity except S. aureus LmrS (96.67%) and S. aureus norA (95.62%) as compared with the % length of reference sequence (Table 1).

The isolate possessed the ARGs such as mepR for glycylcycline and tetracycline antibiotic, mgrA for fluoroquinolone, cephalosporin, penam, tetracycline and peptide antibiotic, arlR for fluoroquinolone and acridine dyes, lmrS for macrolide, aminoglycoside, oxazolidinone, diaminopyrimidine and phenicol, fosB for fosfomycin, norA for fluoroquinolone, glpT and murA transferase for S. aureus with mutation conferring resistance to fosfomycin resistance (Fu et al., 2015; Lee et al., 2020). Raygada & Levine (2009) stated that most CA-MRSA were susceptible to doxycycline, minocycline, clindamycin, trimethoprim-sulfamethoxazole (TMP-SMX), chloramphenicol, rifampin, and linezolid. Whereas, the susceptibility to different antimicrobial agents were reported in varied ranges in different parts of the world (Dennis et al., 2002; Ruhe et al., 2005; Diep et al., 2006; Kallen et al., 2007; Raygada & Levine, 2009).

The three ARO genes mepR, mgrA, and arlR were 100% identical to reference and aligned with 100% coverage to the reference sequence. In addition, analysis of S. aureus glpT determining fosfomycin resistance showed a point mutation which led to a higher level of resistance to fosfomycin. The AMR genes analysed through CARD from the draft genome of MRSA7 strain revealed that they arein correlation with the antimicrobial resistance profiles of disk diffusion tests. Moreover, analysis of the draft genome sequences of this MRSA7 strain showed methicillin resistance belonged to ST28/t021 and is highly resistant to oxazolidinone, sulfonamide, cephems, fluroquinolones and aminoglycoside. In silico analysis revealed high level of conservation in the AMR genes. A few reports from China documented the low incidence of CA-MRSA with high proportionof multidrug resistance (Yao et al., 2010; Wang et al., 2012). It is further suggested that the MDR ST28/t021 type of MRSA strain might have acquired these AMR genes from the workers and possibly circulating in thisdried fish environment. The phylogenetic trees are constructed based on the virulence genes and ARGs (Figs. S2 and S3).The Phylogenetic tree was constructed in NCBI with its data base based on the relations between isolates to virulence genes and ARGs. Data are visualized as minimum spanning tree (MST) and showed that methicillin resistant S. aureus (MRSA 7) strain of this study are closely related to other S. aureus strains mainly MRSA strains in the NCBI data.

Putative virulence factors of CA- MRSA detection

The VirulenceFinder v 2.0 analysis of the core-genome virulence factors MRSA7 strain revealed different virulent factors for cytotoxins, leukotoxins, lekocidins and superantigenic toxins- enterotoxins and virulence factor for exozyme genes such as aureolysin (aur) and serine protease (spiE) (Table 2).

The GOs virulence pathways mainly to the cationic antimicrobial peptide (CAMP) resistance, S. aureus infection (aureolysin) hemolysis by symbiont of host erythrocytes, extracellular component, pathogenesis, toxin activity (hemolysin), peptidase degP and htrA serine protease and GOs for ARGs are two- component system (arlR) for adhesion, autolysis, MDR and virulence genes, metallothioltransferase (fosB), chemical binding (LmrS), MATE family multidrug efflux pump (mepA), major facilitator superfamily (MFS) antibiotic efflux pump (mgrA, norA), N-acetylmuramic acid 6-phosphate etherase (murA) and MFS transporter, DHA2 family, multidrug resistance protein (norB and glpT). The virulence genes of MRSA included gamma-hemolysin chain II precursor (hlgA), gamma-hemolysin component B precursor (hlgB), gamma- hemolysin component C (hlgC), Panton Valentine leukocidin F component (lukF-PV) and Panton Valentine leukocidin S component (lukS-PV) with 100% identity on its genes. The presence of leukocidins- Panton-Valentine leukocidin (PVL) (LukS and LukF proteins) in this MRSA7 strain may cause cytotoxicity to different leukocytes and macrophages and is reported to be associated with CA-MRSAinfections. The present study revealed the presence of leukocidins- Panton-Valentine leukocidin (PVL) (LukS and LukF proteins), reported to be associated with CA-MRSA infections. The incidence of CA- MRSA in seafood was also reported (Sivaraman et al., 2016). The hemolysins viz., hlgA, hlgB and hlgC are well known toxins to blood cells causing cell lysis and cell death. Several reports documented the presence of increased virulence factors in CA-MRSA (Novick & Muir, 1999; Wang et al., 2007; Wardenburg et al., 2007; Li et al., 2010; DeLeo et al., 2011; Hetem et al., 2012). CA-MRSA strains possess a varying repertoire of virulence factors/ toxin genes mainly encoded on mobile genetic elements and could be transferred between strains by horizontal gene transfer and recognized as a potential pathogen.

Staphylococcus aureus are the most frequently occurring food-borne pathogens worldwide with the presence of heat stable preformed staphylococcal enterotoxins (Seo & Bohach, 2007; FDA, 2012) and could cause foodborne intoxication with less than 1.0 microgram toxin from the contaminated food (CDC, 2010). The WGS of this MRSA7 strain (S. aureus 55/2053) shows the presence of Staphylococcal enterotoxins (SEs) such as enterotoxin G (seg), enterotoxin I (sei), enterotoxin M (sem), enterotoxin N (sen), enterotoxin O (seo) and enterotoxin U (seu) with 99.87, 100, 99.86, 97.43, 100, and 100% identity, respectively with the Accession number CP002388.1. S. aureus is considered to be one of the potential food-borne pathogens due to the presence of SEs: SEA to SEE, SEG to SEI, SER to SET (CDC, 2010; Argudín, Mendoza & Rodicio, 2010). SEs may cause vomiting and diarrhoea and the toxins are one of the most common causes of food-borne diseases. The toxins are secreted by entero-toxigenic S. aureus strains in food; they are heat-stable and are not degraded by cooking processes. The SEs are super antigens which trigger T-cell activation and proliferation; their mode of action probably includes activation of cytokine release and cell death via apoptosis and potentially lethal toxic shock syndrome (Balaban & Rasooly, 2000; Lin et al., 2010). This strain also harbours the virulence factor for secreted exoenzyme genes such as aureolysin (aur) and serine protease (spIE) with 100% identity and accession numbers (CP009554.1 and BX571856.1), respectively and they mainly interfere with host metabolic or signalling cascades. The protease aureolysin (S. aureus neutral proteinase) cleaves many proteins including insulin B, with a preference of cleaving after hydrophobic residue and also inactivates PSMs, which cause pathogenesis of osteomyelitis (Larkin et al., 1982). Aureolysin, glutamylendopeptidase, and the cysteine proteases staphopain A and B all interfere with complement factors, leading to evasion of complement-mediated bacterial killing (Baliga et al., 2008). The exfoliative toxins-serine proteases specifically cleave desmosomalcadherins of the superficial skin layers (Recsei et al., 1986; Amagai et al., 2000), leading to staphylococcal scalded skin syndrome (SSSS), a severe skin disease presenting with rash, blisters, and severe lesional damage of the skin. Hence, the co-existence of these different virulence factors for cytotoxins, leukotoxins, lekocidins and superantigenic toxins- enterotoxins and for exozyme genes such as aureolysin (aur) and serine protease (spiE), suggesting that MRSA7 strain may be of significant pathogenic of importance to the public health.

Spa typing and Multi Locus Sequence Typing (MLST)

The availability of the spa types on the central spa server (http://spaserver.ridom.de) and its unified nomenclature provide to understand clonal diversity and transmission of MRSA in the hospital and community settings. The most common spa type in the strain was t021 and its repeats are 15-12-16-02-16-02-25-17-24 in the position between 973–1226 bp in the contig position of MRSA7_gen_1564. Ridom spa server (https://spa.ridom.de/spa302t021.shtml) shows a total of 4568 strains with the frequencies of 1.06% with the comment on CA-MRSA (lukS-luk F) with the CC30 and ST-30 and 3 spa types are reported in India.

MLST analysis of the total genome sequenced (NBZX00000000.1) revealed that MRSA7 strain typed to a ST 28 in fish and fishery products in India with arcC, aroE, glpF, gmk, pta, tpi, ygiL genes with 100% identity and coverage with the alleles of 46, 402, 9, 37, 67, 57, and 557, respectively (Table 3). In the present study, MRSA7 strain Type t021 represents so called ST 28 MRSA, which was found in the fish samples in Gujarat State, India.

The earliest reportedCA-MRSA infections acquired from the communityin 1980 when outbreaks of invasive infections occurred among the intravenous drug users in Detriot (Saravolatz et al., 1982; Levine, Crane & Zervos, 1986) and a unique CA-MRSA was isolated in Australia among aborigines with the history of exposure to antibiotics in 1993. Then the incidence of CA-MRSA was reported in different parts of the world, New Zealand, UK, France, Finland, Canada, USA, and Asia. Major sequence types disseminated across the globe includes ST1, ST5, ST8, ST15, ST30, ST59,ST72, ST80, ST88, ST93, ST121, ST152, ST188, ST398, ST508,ST728, ST834, USA300, etc. (Sowash & Uhlemann, 2014). There is a paucity of information on CA-MRSA in Asian countries as compared to the Europe and USA due to lack of information on the molecular epidemiology of S. aureus and detailed genotyping based on WGS. It is revealed that a low burden of CA-MRSA except ST59 (Taiwan clone) in Taiwan, China and ST239 in South Korea, Malaysia, China, India and Pakistan (Huang & Chen, 2011). A number of studies from India, mainly conducted at hospital centres, suggestedrelatively high prevalence of MRSA, but the contribution of CA-MRSA to S. aureus infections is less clear. To date, strains ST22 and ST772 have been identified as major CA-MRSA clones among infectious isolates and the incidences reported were very low (D’Souza, Rodrigues & Mehta, 2010). However, no report of CA- MRSA in fish is available in India and this is to be a first of its kind based on WGS analysis.

A human pandemic community-associated CC97 lineage MRSA harbouring the antimicrobial resistance genes mecA and mecC has been shown to have originated from animals (Chon, Sung & Khan, 2017). However, in most of the cases reported, a single risk factor could not be ruled out. The presence of MRSA ST28/t021 with luk genes, virulence genes and ARGs in fish is the first kind of report from India. Further comparative genomic analysis provide insights to understand the genetic background information namely genetic structure, dissemination, emergence, genomic diversity of this CA-MRSA strain on virulence, antibiotic resistance. This would further help in the control of such serious public health concern pathogen through the implementation of strict antimicrobial stewardship programme at the local and national level.

Nucleotide sequence accession number

The draft genome sequence of the MRSA7 strain was submitted in the NCBI GenBank under the accession number NBZX00000000.