Serratia marcescens strains, designated here ano1 and ano2, were isolated from two female A. stephensi (LBT and LB1) in our laboratory colony by aseptic dissection of the midgut using sterile tuberculin syringes and needles, followed by plating midgut contents plus sterile saline on Luria-Bertani agar (BD, USA) (Chen et al., 2016). Serratia fonticola strain MSU001 was isolated similarly from a female Aedes triseriatus Say from our laboratory colony. History of mosquito strains is reported elsewhere (Chen et al., 2014, 2015, 2016). Serratia marcescens strain ano1, S. marcescens strain ano2, Serratia fonticola strain MSU001 and their derivatives (see below) were cultured in LB by shaking at 28°C. A dual luciferase expression system was used to mark bacterial strains for quantification. A NanoLuc luciferase-tagged E. anophelis (SCH814) used in this study was previously established and grown on CYE medium (Chen et al., 2015); the firefly system was developed for strains isolated here (see below). Escherichia coli DH5α and E. coli λ pir were used for cloning and conjugative transfer of DNA, respectively. E. coli cells for experiments were grown aerobically in LB broth at 37°C. Bacto agar (Difco, Detroit, MI) was added to a final concentration of 20 g/l. Kanamycin (100 μg/ml) was added for plasmid selection in E. coli, but a higher concentration of kanamycin (600 μg/ml) was required for plasmid selection in Serratia spp.

Isolation and purification of bacterial genomic DNA were performed with the Wizard Genomic DNA Purification Kit (Promega, CA, USA). Integrity and quantity of DNA were assessed using gel electrophoresis, a NanoDrop 2000 UV-Vis spectrophotometer (Thermo Scientific, MA, USA) and Qubit 2.0 fluorometer (Life Technologies, MA, USA), respectively.

The firefly reporter vector was developed based on the wide-range plasmid pBBR1 MCS2 (Chen and Hickey, 2011). The fluc gene was cloned from vector GL4.50 (Promega, WI) with forward primer Walker219 (aggatcctttaagaaggagatatacatatggaagacgccaaaaacataaag) and reverse primer Walker208 (agcatgcttacaatttggactttccgcccttc) with BamHI and SphI on its 5′ and 3′-end, respectively. The amplicon was gel purified and ligated into the T-easy vector (pSCH955). The insert was released from pSCH955 by the restriction enzymes BamHI and SphI and cloned into the same sites on pBBR1 MCS2 downstream of the Plac promoter (pSCH956). Plasmid pSCH956 was introduced into E. coli S17 lamda pir to facilitate conjugation, creating the donor strain for the fluc reporter plasmid (SCH981). SCH981 was mixed with recipient S. marcescens ano1 or S. fonticola MSU001 conjugatively to transfer plasmid pSCH956, leading to firefly reporter strains SCH983 (S. marcescens) or SCH985 (S. fonticola), respectively.

Laboratory rearing procedures for A. stephensi mosquitoes were previously described (Chen et al., 2015). To investigate persistence of infection of bacteria in host mosquitoes, we performed experiments with S. marcescens alone, S. marcescens and E. anophelis together, or S. fonticola and E. anophelis together, using strains with NanoLuc and firefly reporters as described above to quantify presence of bacteria in the mosquito gut lumen. The first experiment was designed to investigate bacterial persistence in live larvae, across molts and metamorphosis past the pupal stage into the adult stage, when bacteria were fed to larvae. When larval A. stephensi reached 3rd instar, individual reporter bacteria SCH983 (S. marcescens with firefly reporter), SCH985 (S. fonticola with firefly reporter) or SCH814 (E. anophelis with NanoLuc reporter) or combinations of them (SCH814/SCH983 or SCH814/SCH985) were added to sterile water in plastic dishes (final concentration adjusted to approximately 4 × 10⁸ CFU/ml) containing 50 larvae, and larvae allowed to feed on the suspension for 24 h, after which the regular larval food regime (Tetra fish food) was continued. Upon metamorphosis, pupae were retrieved, rinsed with sterile water, and transferred into distilled water in cages for adult emergence. After holding the adults for 4 days (during which time they were provided 10% sucrose solution prepared with sterile water), at least 24 adults were randomly sampled and proceeded as described previously (Chen et al., 2015). A second experiment was conducted to investigate persistence of bacterial infection when introduced to adult mosquitoes. The same bacterial preparations as those described in the larval experiment above were fed to adults (1 day of post-emergence) in 10% sucrose solution. After ad lib. feeding for 24 h, the bacterial solution was replaced with fresh, sterile 10% sucrose solution and mosquitoes were sampled at intervals thereafter, and surface-sterilized by immersion in 70% ethanol followed by extensive washings in sterile Milli-Q water. Mosquitoes were homogenized in sterile phosphate buffered saline (PBS), diluted if necessary, and immediately subjected to the Nano-Glo Dual-Luciferase Reporter Assay (Promega, WI).

Susceptibility to different antibiotics was tested by minimal inhibition concentration methods in LB broth (European Committee for Antimicrobial Susceptibility Testing of the European Society of Clinical Microbiology and Infectious Diseases, 2003). Cultures were incubated at 28°C for 24 h under aerobic conditions and the OD_600nm interpreted according to the manufacturer's instructions. Assays were performed in triplicates.

Urease production was tested by using BBL urase test broth (BD, MD, USA). Serratia marcescens strain ano1 or E. coli DH5α was inoculated to the broth and incubated for 48 h.

Hemolysin production in S. marcescens strain ano1 was tested by inoculating bacterial cells on Remel Blood Agar plate (Thermo Scientific, KS). Hemolytic activity was evaluated following incubation at 28°C for 48 h. Staphylococcus aureus and Elizabethkingia meningoseptica were used as the positive controls of beta- and alpha-hemolysin, respectively. Bovine whole blood cells (Hemostat Laboratories, CA) were washed with phosphate buffered saline (PBS) and re-suspended in 1 mL of PBS. Serratia marcescens strain ano1 (final concentration, 4.5 × 10⁸ cells/mL) was incubated with above washed blood cells for 48 h. The erythrocytes were counted using a hemocytometer under microscopy. The negative control was the same as the treatment without introduction of bacteria.

Next generation sequencing (NGS) libraries were prepared using the Illumina TruSeq Nano DNA Library Preparation Kit following standard procedures recommended by the manufacturer. Completed libraries were evaluated using a combination of Qubit dsDNA HS, Caliper LabChipGX HS DNA and Kapa Illumina Library Quantification qPCR assays. Libraries were combined in a single pool for multiplexed sequencing and this pool was loaded on one standard MiSeq flow cell (v2) and sequencing was performed in a 2 × 250 bp paired end format using a v2, 500 cycle reagent cartridge. Base calling was done by Illumina Real Time Analysis (RTA) v1.18.54 and output of RTA was demultiplexed and converted to FastQ format with Illumina Bcl2fastq v1.8.4.

De novo assembly was performed using SPAdes (version 3.9.0) (Bankevich et al., 2012). The reads were assembled using SPAdes (version 3.9.0) and further edited by using DNASTAR (v1.12) (Nurk et al., 2013). Initial prediction and annotation of open reading frames (ORF) and tRNA/rRNA gene prediction were carried out with the Rapid Annotation using Subsystem Technology server (RAST) (Overbeek et al., 2014). Gene annotation was carried out by NCBI Prokaryotic Genome Automatic Annotation Pipeline (PGAAP 3.3) (http://www.ncbi.nlm.nih.gov/genome/annotation_prok/).

The functional categorization and classification for predicted ORFs were performed by RAST server-based SEED viewer (Overbeek et al., 2014). Multi-drug resistance genes were predicted in the Comprehensive Antibiotic Resistance Database (McArthur et al., 2013). Prophage prediction was done with PHAST (Zhou et al., 2011) and Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR) were predicted by CRISPRfinder (Grissa et al., 2007). For genome similarity assessment, average nucleotide identity (ANI) was computed using EDGAR 2.0 (Blom et al., 2016).

The pan genome, core genome, and specific genes of S. marcescens ano1 and ano2 were analyzed by comparing with those in 11 representative Serratia genomes using EDGAR 2.0 (Blom et al., 2016). The development of pan genome and core genome sizes was approximated using the core/pan development feature. The core genome calculated by EDGAR 2.0 was used to infer a phylogeny for the 13 Serratia genomes in this study. The 27,144 amino acid sequences (2,088 per genome) of the core genome were aligned set-wise using MUSCLE v3.8.31 (Edgar, 2004), resulting in a large multiple alignment with 9,015,188 amino acid residues in total (693,476 per genome). This large alignment was used to construct a phylogenetic tree using the neighbor-joining method as implemented in the PHYLIP package (Felsenstein, 1989).

The data from these Whole Genome Shotgun projects have been deposited at DDBJ/ENA/GenBank under the accession. The version described in this paper is version MJVB00000000 and MJVC00000000 for S. marcescens ano1 and ano2, respectively. The BioProject designations for this project are PPRJNA340333 and PRJNA340334, and BioSample accession numbers are SAMN05712591 and SAMN05712592 for S. marcescens ano1 and ano2, respectively.

Statistical analyses were performed using SAS (version 9.2; SAS Institute, Cary, NC).

When S. marcescens ano1 was fed to 3rd instar A. stephensi larvae, the infection persisted to the adult stage for 62.5% (15/24) of adults at 4 days after adult emergence (Table 1), as evidenced by positive luciferase assay relative to negative controls. When S. marcescens ano1 was fed in a sugar meal to adult mosquitoes for 24 h, infection rates of the guts were 98.5% (39/40) when dissected 4 days later (Table 1). By contrast, 17.1% (6/35) and 32.5% (13/40) of the guts of A. stephensi adults were positive for S. fonticola MSU001 when fed to larval and adult stages, respectively, demonstrating that S. fonticola isolated from A. triseriatus had relatively poorer colonization and persistence than did S. marcescens in A. stephensi. When E. anophelis strain was fed to larvae and adults respectively, 87.5% (35/40) and 100% (40/40) infection rates were observed (Table 1). When S. marcescens was introduced together with E. anophelis to the adult host by sugar meal at a cell concentration ratio of 1:1, both of the bacteria had similar infection rates (92.5 and 95%, respectively), using the Dual Luciferase assay. These infection rates were comparable to those in adults introduced with a single bacterial species (Table 1). However, when they were co-introduced in the larval stage, the S. marcescens infection in A. stephensi (22.9%) was lower than that of the single S. marcescens species (62.5%), while NanoLuc-labeled E. anophelis after feeding to larvae reached an infection rate of 100% in adults whether fed with S. marcescens or fed alone. A similar trend was also observed in the co-infection experiment between S. fonticola and E. anopheles (Table 1). These results further demonstrated that infection rate of S. fonticola in A. stephensi was lower than that of S. marcescens when they were co-introduced at either larval or adult stages, suggesting that S. marcescens was co-adapted for infection in A. stephensi. Overall, the rank of infection rate to A. stephensi for the three selected bacteria was E. anophelis > S. marcescens > S. fonticola.

The assembly of strain ano1 contained 44 contigs with a size of 5.45 Mbp (Table 2). The assembly of strain ano2 contained 71 contigs with a size of 5.47 Mbp (Table 2). The ano1 genome included 5,082 coding sequences (CDS) and 119 RNA genes. The ano2 genome included 5,058 CDS and 118 RNA genes (Table 2). Among the 13 Serratia genomes selected for comparative analysis here, ano1 or ano2 had the most CDSs (Table 2). The average GC content for ano1 and ano2 was 59.6 and 59.5%, respectively, consistent with other S. marcescens; however, the average GC content in S. marcescens was much higher than that in other Serratia spp. (Table 2). No plasmid sequence was found in either ano1 or ano2, congruent with our inability to isolate plasmids from the two S. marcescens strains. RAST analysis showed that ano1 and ano2 have at least 589 and 588 subsystems, respectively (Figure S1).

The difference in genome features between ano1 and ano2 was not remarkable (Table 2, Figure 1, and Figure S2). The calculated ANI between ano1 and ano2 shows 100% identity, indicating that S. marcescens ano1 and ano2 were the same strain (Figure 1). ANI values indicated that these two isolates belong to S. marcescens species as they were more than 95% identical to those in S. marcescens MCB (Serepa and Gray, 2014) and S. marcescens Db1 (Flyg et al., 1980), previously isolated from nematodes and fruit flies, respectively (Figure 1). Serratia sp. TEL could be assigned to S. marcescens because it had a high ANI value (>95%) with those of S. marcescens MCB and S. marcescens Db1. It is interesting that, compared to that in S. marcescens FGI94, low ANI values (<89%) in most of the selected Serratia were found, highlighting that different Serratia sp. exist in various insects. Phylogenetic trees (based on genome analysis) showed that ano1 and ano2 were most closely related to Serratia sp. TEL and S. marcescens WW4 (Figure S3). It is interesting that Serratia sp. Ag1 and Ag2 have more CRISPR elements than do other Serratia (Table 2). We did not detect any CRISPR elements in ano1 and ano2. Remarkably, Serratia sp. Ag1 and Ag2 isolated from mosquito A. gambiae have different genome features compared to ano1 and ano2, such as GC content, ANI values and predicted CDS (Table 2, Figure 1, and Figure S3 and see below).

Serratia marcescens strains ano1 and ano2 had 5 predicted prophages (Table S1). Two of the prophages (prophage 3 and 4) were possibly complete because they consisted of tails, heads, portals, integrases, lysins and other component proteins involving phage structure and assembly. Mosquito isolates Ag1 and Ag2 have two incomplete prophages. The number of predicted prophages in other Serratia ranged from 1 to 8. Collectively, these prophages varied in size and gene organization, indicating their diversity in Serratia (Table S1).

The gene repertoire of the selected Serratia species was analyzed using their ubiquitous genes (core genome) and different homologous genes families (pan-genome) amongst the selected Serratia genomes. The pan-genome plot shows that the power trend line had not reached a plateau (Figure 2A), demonstrating that Serratia displays an open pan-genome. Serratia core genome analysis showed that the number of shared genes decreased with the addition of the input genomes and was predicted to converge against 2,127 (see Figure 2B). Singleton development plot data indicated that up to 240 new genes could be expected with every newly sequenced genome (Figure 2C). The core genome for the 13 selected Serratia was calculated to be 2,088 CDS per genome; given the assumption that the approximation slightly over-predicts the real core genome size, the current core genome likely represents the Serratia genus quite well.

Serratia marcescens ano1 and ano2 shared at least 4985 genes (data not shown). Only 13 and 12 genes were uniquely present in strains ano1 and ano2, respectively, a result within the error rate of incomplete sequencing and likely artifact. Also, among its relatives, the ano1 genome shared in common genes ranging from 2936 to 3749 in number, with the lowest number shared with Serratia sp. Ag2, indicating diverse physiological functions in the selected Serratia (Figure 3 and Figure S4). For example, S. marcescens ano1 shared at least 4198, 4178, 4172, and 3788 common genes with strains S. liquefaciens, S. marcescens Db11, S. marcescens MCB and Serratia sp. TEL, which accounts for approximately 84.0, 83.6, 83.5, and 75.8% of its encoding genes, respectively; the above 5 selected S. marcescens shared 3188 common genes (Figure 3).

Serratia marcescens ano1 was resistant to nearly all of the selected antibiotics including class β-lactams, aminoglycosides, tetracycline, aminoglycosides, macrolides, glycopeptides, and ansamycins, showing that ano1 is a multi-drug resistant strain (Table 3). At least 18 predicted enzymes/proteins conferring antibiotic resistance were predicted by CARD and RAST SEED subsystem (Figure S1 and Table S2). The predicted proteins include those conferring antibiotic resistance to β-lactams, fosfomycin, mupirocin, polymyxin, aminocoumarin, aminoglycoside, isoniazid and fluoroquinolone. Furthermore, up to 37 efflux pumps were identified, possibly contributing to tetracycline, rifampin, aminocoumarin, aminoglycoside, fluoroquinolone, macrolide as well as β-lactam resistance (Table S2). Comparative study of Serratia genomes showed that all of them have potential resistance against four classes of antibiotics: fluoroquinolone, polymyxin, mupirocin, and fosfomycin.

Most of the S. marcescens strains had genes encoding antimicrobial compound synthesis proteins such as hydrogen cyanide (OHT38223), bacteriocin (colicin V, OHT33745) and pyoverdine (OHT34953) (also see next), indicating that they have potential to inhibit the growth of similar or closely related bacterial strain(s) (Table S3). Ano1 and ano2 carried the gene cluster of biosynthesis of bacteriocin, consisting of 19 opening reading frames (ORFs). Most of them were conserved among the Serratia genomes examined here, including those isolated from mosquitoes and nematodes (Table S3). It is noteworthy that prodigiosin gene synthesis clusters were absent in S. marcescens ano1, also lacking in several insect symbionts (Table S4).

Strains ano1 and ano2 had a probability of acting as a pathogen of 72.9% according to the probability values assigned by PathogenFinder (Cosentino et al., 2013). Ano1 genome matched 41 pathogenic families and 12 non-pathogenic families in the database (data not shown). Furthermore, a total of 37 putative virulence factors were predicted by VFDB, which mainly account for flagella formation, lipooligosaccharide (LOS), attacking enzymes, iron uptake and transportation, different secretion systems and hemolysis of red blood cells (Table 4). It is interesting that several predicted virulence factors in ano1 and ano2 such as gluP, manC, urase, shu, and fleQ were either missed or showed low identity (<40%) in mosquito symbionts Serratia sp. Ag1 and Ag2 (Table 4). For example, gluP (encoding L-fucose:H+ symporter permease, BGV45_11240), participates in chronic infection and acts an important virulence factor in Brucella abortus (Xavier Mariana et al., 2013). manC, a mannose-1-phosphate guanylyltransferase gene (forming O-antigen) (Lee et al., 1992), was absent in Ag2 though there is a manC-like gene in Ag1 showing a significantly low identity (37%). This manC gene is absent in environmental isolate S. marcescens FGI94. An iron transporter gene, shu (Mourer et al., 2015), existed in most of the selected Serratia strains (>87%) while it was absent in those in Ag1 and Ag2 (Table 4).

At least 4 different chitinase genes encode chitinase A (OHT33328 and OHT32204), chitinase B (OHT36384) and chitinase C (OHT38806) in S. marcescens ano1 (Table 5). ChiA digests chitins from the reducing end while ChiB works on the chitin chain from the non-reducing end, indicating that they are processive enzymes (Suzuki et al., 2002; Orikoshi et al., 2005; Vaaje-Kolstad et al., 2013). Instead, ChiC acts as an endo-acting enzyme which cuts chitin chain in the middle (Vaaje-Kolstad et al., 2013). ChiC exists in most of the selected Serratia. However, an efficient chitin degradation requires a synergistic action of ChiA, ChiB, and ChiC in S. marcescens (Orikoshi et al., 2005; Vaaje-Kolstad et al., 2013). Besides the different catalytic domains (all belonging to GH18 superfamily) (Tian et al., 2014), ChiA1, ChiB and ChiC chitinases have various accessory functional domains (Table 5) such as fibronectin type III (OHT33328), chitin-binding (OHT38806), or cellulose-binding (OHT36384) domains while ChiA2 has only catalytic one (OHT32204) (Tran et al., 2011). Previous studies showed that CBD or chitin-binding domains greatly affected the catalytic activity and substrate affinity in chitinases (Dahiya et al., 2006; Tian et al., 2014). Therefore, the chitinases with various functional domain(s) may provide different physiological roles in S. marcescens (Dahiya et al., 2006). The four chitinase genes dispersedly spread in S. marcescens ano1 genomes, rather than cluster together to form an operon (data not shown). Such arrangement indicates that they may be induced under the different conditions and regulated by different mechanisms. No genes encoding chitinase (cutoff, 60% identity) were found in Serratia sp. Ag1 and Ag2, indicating that they have a poor ability to use chitin as carbon or nitrogen sources from environment or in insects. As expected, in S. marcescens FGI94, a lack of ChiA1 and ChiB as production of large amount of chitinase will disrupt its symbiotic relationship with fungi (by degrading chitins in fungal cell walls; Li P. et al., 2015). Instead, similar to ano1 or ano2, Db1, MCB, and TEL have ChiA1, ChiB, and ChiC. WW4 has the complete chitin degradation enzyme system, which is consistent with its living niche (paper machine) where the plant materials are abundant (Chung et al., 2013).

Serratia marcescens ano1 and ano2 were hemolytic strains with alpha-hemolysin activity (Figure 4A). Red blood cell lysis experiment was further conducted by inoculating S. marcescens ano1 in bovine whole cells (Figure 4B). At least 16% of erythrocytes were disrupted within 48 h (Figure 4B). The commensal S. marcescens living in mosquito midgut may not utilize the hemolysin(s) for pathogenesis unless they accidently invade the host through liaison and enter the insect hemolymph (Chen et al. unpublished data).

Serralysin genes are present in most of the insect associated Serratia (S. marcescens MCB, S. marcescens Db1 and Serratia sp. TEL) with high identity in amino acid sequences (>60%) (Table 6). S. marcescens ano1 possesses at least 5 genes encoding serralysins or serralysin-like proteins (OHT33392, OHT35629, OHT35861, OHT36525, and OHT38015) (Table 6). It is interesting that Serratia sp. Ag1 and Ag2 have only two of them (Table 6). At least three serralysin genes are found in plant or fungi symbionts S. marcescens, S. plymuthica S13, and S. marcescens FGI94, respectively. In this study, we identified an additional one based on the genome mining method (Table 6).

ureA, encoding an α-subunit ureases, was predicted as the virulence factor (Table 3). Among the selected Serratia, only ano1, ano2 and TEL carry ureA (Table 3). Further examination of the three selected Serratia genomes showed that there is an operon consisting of 9 genes possibly involved in urea metabolism (Figure S5). Serratia marcescens strain ano1 was confirmed to urease production positive (data not shown), indicating that the “ure” operon is functional in ano1. In this “ure” operon, three genes (ure A, B, and C) were predicted to encode α, β, and γ subunits, respectively (Carlini and Ligabue-Braun, 2016). Genes ure D, E, F, and G encode urease accessory proteins that participate in assembly and activation of ureases (Carlini and Ligabue-Braun, 2016). Additionally, urea transporter (utp) and nickel transporter (nixA) genes immediately locate downstream of urease structure protein or accessory protein encoding genes (Carlini and Ligabue-Braun, 2016). However, the ure operon is only found in 14 out of 101 Serratia genomes publically available in the IMG/MER (cutoff date, Dec.1, 2016). We found mosquito or nematode-associated Serratia (ano1, ano2, TEL) with this operon.

Enterobactin synthesis gene clusters and their organization were conserved in all S. marcescens genomes (Table 7). Several operons encoding ferrichrome ABC transporters (fhu operon), ferric siderophore ABC transporters (fep operon) and ferric citrate outer membrane transporters (fec operon) were found (Table 7). Serratia marcescens ano1 and ano2 had a delicate set of heme uptake/storage systems (Table 7). The has operon encoding the heme uptake system includes an RNA polymerase sigma factor (OHT42022.1), putative iron sensor protein (OHT42023.1), TonB-dependent heme receptor (OHT42024.1), hemophore HasA (OHT42025.1), peptidase (OHT42071.1), hemolysin transporter (OHT42026.1), and TonB-like protein (OHT42027.1). This system ensures that, after heme/hemin is released from hemoglobins or other iron containing proteins, they can be efficiently processed, stored and transported. The ferric uptake regulator (Müller et al., 2013) has a key role in modulating iron uptake, and the genome of S. marcescens strain ano1 and ano2 is predicted to encode a Fur protein (OHT38796.1). However, the fur gene is not close to any iron metabolism genes in S. marcescens (data not shown).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.