Genomic identification and characterization of co-occurring Harveyi clade species following a vibriosis outbreak in Pacific white shrimp, Penaeus (litopenaeus) vannamei

Journal Pre-proof Genomic identification and characterization of co-occurring harveyi clade species following a vibriosis outbreak in Pacific white shrimp, Penaeus (litopenaeus) vannamei Paxton T. Bachand, James J. Tallman, Nicole C. Powers, Megan Woods, Danial Nasr Azadani, Paul V. Zimba, Jeffrey W. Turner PII:

S0044-8486(19)31887-3

DOI:

https://doi.org/10.1016/j.aquaculture.2019.734628

Reference:

AQUA 734628

To appear in:

Aquaculture

Received Date: 25 July 2019 Revised Date:

18 September 2019

Accepted Date: 21 October 2019

Please cite this article as: Bachand, P.T., Tallman, J.J., Powers, N.C., Woods, M., Azadani, D.N., Zimba, P.V., Turner, J.W., Genomic identification and characterization of co-occurring harveyi clade species following a vibriosis outbreak in Pacific white shrimp, Penaeus (litopenaeus) vannamei, Aquaculture (2019), doi: https://doi.org/10.1016/j.aquaculture.2019.734628. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

1

Genomic identification and characterization of co-occurring Harveyi clade species

2

following a vibriosis outbreak in Pacific white shrimp, Penaeus (litopenaeus) vannamei

3 4

Paxton T. Bachanda, James J. Tallmana, Nicole C. Powersa, Megan Woodsa, Danial Nasr

5

Azadania, Paul V. Zimbab and Jeffrey W. Turnera,*

6 7

a

8

USA

9

b

Department of Life Sciences, Texas A&M University-Corpus Christi, Corpus Christi, Texas,

Center for Coastal Studies, Texas A&M University-Corpus Christi, Corpus Christi, Texas, USA

10 11

*Corresponding author: Jeffrey W. Turner, Department of Life Sciences, 6300 Ocean Drive,

12

Unit 5858, Texas A&M University-Corpus Christi, Corpus Christi, Texas, 78412 USA. Tel: +1-

13

361-825-6206; Fax: +1-361-825-2742; E-mail: [email protected]

14 15

Abstract

16

Shrimp farming is one of the fastest-growing food production sectors but disease related

17

mortality remains a limitation to industry expansion. Bacterial species belonging to the Harveyi

18

clade, such as Vibrio harveyi and V. parahaemolyticus, have long been implicated as agents of

19

disease in shrimp aquaculture but traditional methods cannot accurately differentiate these

20

species. Following a vibriosis outbreak in Pacific white shrimp, Penaeus (litopenaeus) vannamei,

21

we hypothesized that infection promoted the co-occurrence of multiple Harveyi clade species.

22

The shrimp were reared in a zero-exchange raceway using biofloc technology and presumptive

23

Vibrio species were isolated from the hepatopancreases of moribund shrimp by culture on a

1

24

Vibrio-selective agar. Isolates were identified and characterized by PCR-based genotyping and

25

whole-genome sequencing. A phylogenetic tree based on a 43-genome pangenome analysis

26

revealed the presence of five known shrimp pathogens: Photobacterium damselae, V.

27

alginolyticus, V. harveyi, V. parahaemolyticus, V. rotiferianus, and a potentially novel Vibrio

28

species (Vibrio sp. Hep-1b-8). Multilocus sequence typing revealed that P. damselae and V.

29

parahaemolyticus isolates exhibited novel allelic profiles, although the V. parahaemolyticus

30

isolates were closely related to strains isolated from Southeastern Asia. Further, a novel

31

comparative analysis of 23 type VI secretion systems revealed three phylogenetically distinct

32

systems, which suggests that co-occurring bacteria possess diverse mechanisms for interspecies

33

competition. These results highlight whole-genome sequencing as an invaluable tool for

34

identifying and characterizing co-occurring, infection-associated bacteria at the species, strain

35

and molecular level.

36 37

Keywords: Pacific white shrimp; Penaeus (litopenaeus) vannamei; vibriosis; infection; whole-

38

genome sequencing; Vibrio; Photobacterium

39 40

1. Introduction

41

Aquaculture is the world’s fastest-growing food production sector and Pacific white shrimp,

42

Penaeus (litopenaeus) vannamei, are among the most successfully farmed aquatic species

43

(Kumar and Engle, 2015). Human population growth and a rising demand for seafood-derived

44

protein are predicted to increase farmed shrimp demand (Msangi et al., 2013). A primary

45

challenge to meeting this demand will be the prevention and control of disease (Cock et al.,

2

46

2015; Xiong et al., 2016). To this end, the accurate identification and characterization of

47

infection-associated bacteria is a critical need.

48 49

Bacteria belonging to the Vibrio genus are known agents of recurring and emerging disease in

50

shrimp aquaculture (Lightner and Redman, 1998). Vibrio carchariae, an early taxonomic

51

synonym for Vibrio harveyi, was first described as a fish pathogen nearly four decades ago

52

(Colwell and Grimes, 1984). It remains an important shrimp pathogen (Austin and Zhang, 2006)

53

causing luminescent vibriosis (Lavilla-Pitogo et al., 1990) and bright-red syndrome (Soto-

54

Rodriguez, Gomez-Gil and Lozano, 2010). More recently, Vibrio parahaemolyticus strains

55

harboring a novel virulence plasmid were shown to cause acute hepatopancreatic necrosis

56

disease (AHPND) (Soto-Rodriguez et al., 2015).

57 58

Recent studies utilizing small subunit ribosomal RNA (SSU rRNA) gene sequencing have shown

59

that infection disrupts bacterial community structure, making the gut bacteria of diseased shrimp

60

less diverse and more susceptible to opportunistic pathogens (Dai et al., 2018; Yao et al., 2018).

61

Recent studies have also reported that AHPND infection lowered intestine and hepatopancreas

62

bacterial diversity and promoted the appearance of infection-associated bacterial species (Chen

63

et al., 2017; Cornejo-Granados et al., 2017). While these studies have provided valuable insight

64

into the complexities of shrimp disease, SSU rRNA sequencing lacks resolution at the species,

65

strain and molecular level (Rossi-Tamisier et al., 2015). In particular, closely related Vibrio

66

species belonging to the Harveyi clade (e.g. Vibrio alginolyticus, Vibrio campbellii, Vibrio

67

diabolicus, V. harveyi, Vibrio jasicida, Vibrio natriegens, Vibrio owensii, V. parahaemolyticus

3

68

and Vibrio rotiferianus) cannot be accurately identified by SSU rRNA gene sequencing (Cano-

69

Gomez et al., 2011; Urbanczyk et al., 2013; Ke et al., 2017; Turner et al., 2018).

70 71

Whole-genome sequencing (WGS) is a high-resolution technology at the species, strain and

72

molecular level that is predicted to become an indispensable tool in aquaculture (Bayliss et al.,

73

2017). In this study, we used WGS to identify and characterize co-occurring bacteria associated

74

with a vibriosis outbreak in Pacific white shrimp (P. vannamei). This study did not however seek

75

to inform etiology as this outbreak was reported previously and the presence of V.

76

parahaemolyticus was implicated in reduced shrimp survival (Prangnell et al., 2016). Rather, this

77

study hypothesized that 1) infection promoted the co-occurrence of multiple Harveyi clade

78

species and 2) co-occurring species would possess distinct molecular mechanisms that may

79

facilitate interspecies competition.

80 81

2. Materials and methods

82

2.1. Sample collection

83

Bacteria cultures were provided by Drs. David Prangnell and Tzachi Samocha (Texas A&M

84

University AgriLife Research and Extension Center at Corpus Christi, TX, USA) following a

85

vibriosis outbreak in Pacific white shrimp (P. vannamei). The outbreak was reported previously

86

and the presence of Vibrio species, specifically V. parahaemolyicus, was implicated in reduced

87

survival (Prangnell et al., 2016). The shrimp were reared for 49 days with non-medicated feed in

88

a recirculating raceway (40 m3) using biofloc technology at the Texas A&M University AgriLife

89

Research and Extension Center. Briefly, hepatopancreas samples were collected aseptically from

90

four moribund shrimp (~20 g in weight) at the end of the grow out phase. Samples were

4

91

homogenized, serial-diluted with PBS, and plated on CHROMagar Vibrio media (CHROMagar,

92

Paris, France). Twenty-two bacteria colonies were randomly selected and subcultured for

93

isolation on thiosulfate-citrate-bile salts-sucrose (TCBS) agar (Oxoid, Hampshire, England) at 30

94

˚C overnight (18 hours). Isolated bacteria were archived as frozen cultures (20% glycerol final

95

concentration, -80 ˚C) for further analysis.

96 97

2.2. Genotypic identification

98

Nineteen of the 22 frozen cultures (above) were successfully resuscitated. The identities of the

99

19 isolates were assessed by the polymerase chain reaction (PCR) detection of the V.

100

parahaemolyticus thermolabile hemolysin gene (tl) (Bej et al., 1999) and the ß-lactamase gene

101

(blaCARB-17) (Li et al., 2016) using a QIAGEN Fast Cycling PCR kit (Qiagen, Valencia, CA,

102

USA). Total DNA was isolated via heat lysis (Englen and Kelley 2000) and assayed for quantity

103

(ng µL-1) and quality (A260/A280) using a BioPhotometer D30 (Eppendorf, Hamburg, Germany).

104

The cycling conditions were as follows: initial denaturation of 95 ˚C for 5 minutes followed by

105

40 cycles of 96 ˚C for 5 seconds, 55 ˚C (tl) or 50 ˚C (blaCARB-17) for 5 seconds and 68 ˚C for 10

106

seconds, followed by a 1-minute extension at 72 ˚C. The pandemic type strain V.

107

parahaemolyticus RIMD2210633 (Nasu et al., 2000) and a no template control (nuclease free,

108

PCR grade water) were included as positive and negative controls, respectively. Amplicons were

109

separated by electrophoresis on a 1% agarose gel (80V for 35 minutes) and visualized using

110

SYBR Safe stain (Invitrogen, Carlsbad, CA, USA) under UV light.

111 112

2.3. Whole-genome sequencing

5

113

The cost of sequencing all 19 of the bacteria isolates (above) was prohibitive, so 13 were

114

randomly selected for WGS. These bacteria were grown overnight (18 hours) at 30 ˚C in tryptic

115

soy broth (TSB) (Becton, Dickinson and Company, Franklin Lakes, NJ, USA) while shaking

116

(100 rpm). Cells were pelleted by centrifugation (8,000 g for 10 minutes) and washed twice with

117

1 mL PBS. Genomic DNA was isolated from the pelleted cells using a ChargeSwitch gDNA

118

Mini Bacterial Kit (Invitrogen, Carlsbad, CA, USA). The DNA was quantified (ng µL-1) using a

119

Qubit 4 Fluorometer (Invitrogen, Carlsbad, CA, USA) and stored at -20 ˚C. Two isolates (Hep-

120

2a-10 and Hep-2a-11) were sequenced at the New York University (NYU) Genome Technology

121

Center (New York, NY, USA) with an Illumina MiSeq instrument using 2 x 300 bp paired-end

122

chemistry. The remaining 11 isolates (Hep-1a-1, Hep-1a-2, Hep-1a-3, Hep-1b-8, Hep-2a-12,

123

Hep-2a-14, Hep-2a-16, Hep-2b-18, Hep-2b-19, Hep-2b-20 and Hep-2b-22) were sequenced at

124

Molecular Research LP (Shallowater, TX, USA) with an Illumina HiSeq 2500 instrument using

125

2 x 150 bp paired-end chemistry. Sequence reads from Hep-2a-10 and Hep-2a-11 were processed

126

and assembled de novo as described previously (Moreno et al., 2017). Sequence reads from the

127

remaining isolates were assembled de novo with the Maryland Super-Read Celera Assembler

128

(MaSuRCA) (Zimin et al., 2013). The draft genomes were annotated with the NCBI Prokaryotic

129

Genome Annotation Pipeline (Tatusova et al., 2016).

130 131

2.4. Genomic identification

132

The genomic identification of the 13 sequenced isolates was based on 1) a genome-scale

133

phylogenetic analysis and 2) a pair-wise comparison of average nucleotide identity (ANI) values.

134

The 43 genomes included in these analyses, including the 13 genomes from this study, were

135

listed in Table S1 in the online supplementary material. To construct the genome-scale

6

136

phylogeny, homologous genes were clustered with get_homologues (options -M -t 0 -e 1) and a

137

reduced binary matrix, representing the pangenome matrix, was compiled with the

138

accompanying script compare_clusters.pl (Contreras-Moreira and Vinuesa, 2013). The binary

139

matrix was used to generate a maximum-likelihood (ML) tree with IQ-TREE version 1.5.5

140

(Nguyen et al., 2014) with 1000 ultrafast (UF) bootstraps (Minh et al., 2013) using the best-fit

141

model as determined by ModelFinder (Kalyaanamoorthy et al., 2017). The topology of the tree

142

was tested by constructing a second ML tree with RAxML version 8.1.22 with 100 fast

143

bootstraps (-f a -m BINGAMMA -p 12345 -x 12345 -# 100) (Stamatakis, 2014). The ML trees

144

were annotated with FigTree version 1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). The ANI

145

values were calculated using JSpeciesWS (Richter et al., 2016) using 96% similarity as the

146

threshold for identifying different prokaryotic species.

147 148

2.5. Multilocus sequence typing

149

The multilocus sequence type (MLST) of the P. damselae and V. parahaemolyticus strains was

150

inferred in silico using established MLST schemes (Gonzalez-Escalona et al., 2008; Alba et al.,

151

2016). Housekeeping loci were retrieved with blastn version 2.9.0 (Camacho et al., 2009),

152

aligned with MUSCLE version 3.6 (Edgar, 2004) and trimmed with trimAl version 1.4.15

153

(Capella-Gutierrez et al., 2009). All alleles and metadata were deposited in the respective P.

154

damselae and V. parahaemolyticus PubMLST databases at https://pubmlst.org/databases/.

155 156

2.6. Detection and relatedness of protein secretion systems

157

Type VI secretion systems (T6SSs) were detected with T346Hunter (Martínez-García et al.,

158

2015) using default search parameters (blastp and HMMER3 E-value 0.0005). The complete

7

159

genome of the pandemic type strain V. parahaemolyticus RIMD2210633 (Makino et al., 2003),

160

which harbors two distinct T6SSs (T6SS1 and T6SS2, Park et al., 2004), was included as a

161

reference. The protein sequences corresponding to core system components (Shrivastava and

162

Mande, 2008), were saved for multilocus sequence analysis (MLSA). A subset of conserved core

163

protein sequences was identified with get_homologues (options -M -t 23 -e) (Contreras-Moreira

164

and Vinuesa, 2013). The protein sequences were individually aligned with MUSCLE version 3.6

165

(Edgar, 2004), trimmed with trimAl version 1.4.15 (Capella-Gutierrez et al., 2009), concatenated

166

with Galaxy (Giardine et al., 2005) and a ML tree was constructed as described above. The

167

synteny of the phylogroups was evaluated using the open reading frame (ORF) maps produced

168

by T346Hunter. Genes were named according to their common vernacular (e.g. vgrG, hcp, vasH,

169

vasI and vasL).

170 171

3. Results

172

3.1. Sample collection

173

To isolate bacteria associated with the infection, homogenized hepatopancreas samples were

174

plated in duplicate on CHROMagar Vibrio, resulting in more than 200 colony forming units

175

(CFU) per sample. The appearance of those colonies ranged from mauve (indicative of V.

176

parahaemolyticus), turquoise (indicative of V. vulnificus and V. cholerae) and colorless

177

(indicative of V. alginolyticus). Nineteen colonies with diverse morphologies were randomly

178

selected for isolation on TCBS agar and cryopreserved for further analysis.

179 180

3.2. Genotypic identification

8

181

The PCR-based detection of the tl gene identified 7/19 isolates as V. parahaemolyticus: Hep-1a-

182

2, Hep-1a-4, Hep-1b-6, Hep-1b-7, Hep-2a-12, Hep-2b-20 and Hep-2b-21. The detection of the

183

blaCARB-17 gene identified three isolates as V. parahaemolyticus: Hep-2a-12, Hep-2b-20 and Hep-

184

2b-21.

185 186

3.3. Whole-genome sequencing

187

The 13 draft genome assemblies were deposited at GenBank under the accession numbers given

188

in Table S1 in the online supplemental material. The general attributes of those genomes were

189

summarized in Table S2 in the online supplemental material.

190 191

3.4. Genomic identification

192

The identity of the 13 strains was determined by a 43-genome pangenome analysis and pairwise

193

ANI comparisons. The reduced pangenome matrix contained 21 159 binary sites that was

194

representative of 21 169 ORF clusters. The binary alignment contained 1261 constant sites, 9781

195

parsimony informative sites and 3682 distinct site patterns. Importantly, no two alignments were

196

identical and thus all strains, even strains of the same species, were considered distinct. The best-

197

fit model according to the Bayesian information criterion (BIC) scores and weights was

198

calculated as GTR2+F0+G4. The ML phylogeny resolved the 13 bacterial strains into two genera

199

and six species: P. damselae subsp. damselae (N = 4), V. alginolyticus (N = 1), V. harveyi (N =

200

1), V. parahaemolyticus (N = 2), V. rotiferianus (N = 4) and Vibrio sp. Hep-1b-8 (N = 1) (Figure

201

1). The ML phylogeny confirmed that the blaCARB-17 assay correctly identified V.

202

parahaemolyticus Hep-2a-12 and Hep-2b-20. The blaCARB-17 assay also identified Hep-2b-21 as

203

V. parahaemolyticus but that strain was not sequenced. Phylogenetic assignments were

9

204

confirmed by ANI values greater than 96% with representative strains of each species (Table 1).

205

Pairwise comparisons of P. damselae strains with P. damselae subsp. damselae KC-Na-1

206

(>97%) and P. damselae subsp. piscicida 91-197 (>96 but <97%) supported assignment to the

207

damselae subspecies. Vibrio sp. Hep-1b-8 appeared to be a novel species as the most closely

208

related strain in the GenBank database – V. brasiliensis LMG 20456 – shared only 82.38% ANI

209

but the systematics of this strain was beyond the scope of this study.

210 211

3.5. Multilocus sequence typing

212

A MLST analysis of P. damselae and V. parahaemolyticus strains tested relatedness against

213

global MLST databases. The four P. damselae strains exhibited identical but novel allelic

214

profiles (glpF 17, gyrB 25, metG 22, pntA 18, pyrC 26 and toxR 30) and sequence type (ST56).

215

This allelic profile included five new alleles and was therefore distinct in comparison to other

216

sequence types in the database. The two V. parahaemolyticus strains also exhibited identical but

217

novel allelic profiles (dtdS 126, dnaE 171, gyrB 222, pntA 4, pyrC 420, recA 113 and tnaA 23)

218

and sequence type (ST1931). The V. parahaemolyticus ST1931 sequence type differed from the

219

previously identified ST428 (isolated from a non-described environmental source in Vietnam)

220

and ST1226 (isolated from a non-described environmental source in Philippines) by single

221

mutations in the pyrC and recA alleles, respectively.

222 223

3.6. Detection and relatedness of protein secretion systems

224

To investigate mechanisms of interspecies competition, type VI secretion systems (T6SS) were

225

detected and compared. All 13 strains possessed at least one system and two distinct systems

226

(i.e., T6SS1 and T6SS2) were detected in the closed reference genome. In total, 23 T6SSs

10

227

containing between 14 and 18 core components each were detected. The get_homologues

228

analysis showed that six of those core components were present and conserved in all 23 systems:

229

ClpV, VasA, VasE, VasK, VipA and VipB (Table 2). The concatenated multilocus protein

230

alignment was comprised of 23 sequences, 3791 amino acids, 3181 parsimony informative sites

231

and 3245 distinct site patterns. The best-fit model according to the Bayesian information

232

criterion (BIC) scores and weights was calculated as LG+F+I+G4. The ML phylogeny resolved

233

the 23 T6SSs into three major clades that were each further resolved into two minor clades:

234

T6SS1_A and T6SS1_B, T6SS2_A and T6SS2_B, and T6SS3_A and T6SS3_B (Figure 2).

235

Clades 1 and 2 were homologous to T6SS1 and T6SS2, respectfully. Clade 3 represented a

236

distinct system. The major clades included multiple species: T6SS1: P. damselae, V.

237

alginolyticus, V. harveyi, V. parahaemolyticus and V. rotiferianus, T6SS2: V. alginolyticus, V.

238

harveyi, V. parahaemolyticus and V. rotiferianus and T6SS3: V. harveyi and V. rotiferianus.

239

Three of the minor clades were restricted to a single species: T6SS1_B containing P. damselae,

240

T6SS2_B containing V. rotiferianus and T6SS3_A containing V. rotiferianus. V. alginolyticus

241

Hep-1a-2 was the only strain that possessed systems homologous to both T6SS1 and T6SS2. The

242

gene organization of the T6SSs showed strong agreement with the major and minor phylogenetic

243

clades (see Figure S1 in the online supplemental material).

244 245

4. Discussion

246

Previous studies of vibriosis outbreaks have shown that diseased shrimp are colonized by

247

multiple co-occurring Vibrio species (e.g. Leaño et al., 1998; Rodriguez et al., 2010; Yao et al.,

248

2018; Dai et al., 2018). Yet little is known about the diversity of co-occurring species belonging

249

to the Harveyi clade as the techniques commonly applied (e.g., phenotypic tests and SSU rRNA

11

250

sequencing) cannot accurately identify closely related members of the Harveyi clade. In this

251

study, we hypothesized that a vibriosis outbreak in P. vannamei promoted the co-occurrence of

252

multiple Harveyi clade species and we tested this hypothesis using WGS.

253 254

The identification of some Harveyi clade species can be achieved through the PCR-based

255

detection of conserved DNA sequences (e.g., Chatterjee and Haldar 2012). For example, a

256

multiplex assay, targeting the hly hemolysin gene, was developed for the simultaneous detection

257

and differentiation of V. campbellii, V. harveyi, and V. parahaemolyticus (Haldar et al., 2010).

258

Similarly, another multiplex assay targeting the topA topoisomerase, the ftsZ GTPase, and the

259

mreB ATPase, was developed for the simultaneous detection and differentiation of V. campbellii,

260

V. harveyi, V. owensii, and V. rotiferianus (Cano-Gomez et al., 2015). While these assays are

261

advantageous, the accuracy of ‘species-specific’ assays have been subject to revision; the

262

primers targeting the thermolabile hemolysis gene (tl) of V. parahaemolyticus (Bej et al., 1999)

263

have been shown to amplify sequence variants in V. alginolyticus, V. diabolicus, V. tubiashii, and

264

P. damselae (Klein et al., 2014; Yáñez et al., 2015; Turner et al., 2018).

265 266

V. parahaemolyticus was implicated as a causal agent of the vibriosis infection reported in this

267

study (Prangnell et al., 2016). Thus, we tested presumptive Vibrio strains for two V.

268

parahaemolyticus species-associated markers: the tl gene and the blaCARB-17 gene. The tl assay

269

was not accurate while the blaCARB-17 assay proved accurate and correctly identified two strains

270

as V. parahaemolyticus. The identity of one or more of the remaining strains could have been

271

determined using either of the multiplex assays (described above) designed to detect and

272

differentiate V. campbellii, V. harveyi, V. owensii, V. parahaemolyticus, and V. rotiferianus, but

12

273

to our knowledge no single assay is capable of detecting and differentiating all Harveyi clade

274

species (e.g., V. alginolyticus, V. campbellii, V. diabolicus, V. harveyi, V. jasicida, V. natriegens,

275

V. owensii, V. parahaemolyticus, and V. rotiferianus). Additionally, said assay would fail to

276

detect the presence of unknown Harveyi clade species as well as non-Harveyi clade species.

277 278

Whole-genome sequencing can accurately identify and differentiate all known Harveyi clade

279

species (Lin et al., 2009; Hoffmann et al., 2011; Urbanczyk et al., 2013; Ke et al., 2017; Turner

280

et al., 2018). It also eliminates the risk of false positives resulting from primer mishybridization

281

and it provides invaluable data at the strain and molecular level. For these reasons, combined

282

with the steady decrease in sequencing costs and the increased availability of bioinformatics

283

software and pipelines, WGS has been described as the ‘ultimate’ typing methodology that will

284

become an indispensable tool in aquaculture (Bayliss et al., 2017). In this study, the use of WGS-

285

based analyses clearly showed that moribund shrimp were colonized by at least five bacterial

286

species (i.e., P. damselae, V. alginolyticus, V. harveyi, V. parahaemolyticus, and V. rotiferianus)

287

that are known shrimp pathogens (Liu et al., 2004; Austin and Zhang, 2006; Rivas et al., 2013;

288

Zhang et al., 2014; Lee et al., 2015).

289 290

Aquaculture infections are increasingly thought to be polymicrobial (Kotob et al., 2016). Recent

291

reports have demonstrated that shrimp infections are accompanied by a reduction in bacterial

292

community diversity and the appearance of opportunistic pathogens (Dai et al., 2018; Yao et al.,

293

2018), and a recent commentary suggested that a growing number of diseases in marine

294

ecosystems may result from or be complicated by dysbiosis (Egan and Gardiner, 2016). This

295

thinking is a departure from the traditional monomicrobial understanding of infectious disease

13

296

and it is becoming increasingly clear that Koch’s postulates should be adapted to polymicrobial

297

scenarios (Nelson et al., 2012). Future efforts, such as 1) time-series data detailing the

298

progression of the infection and the appearance of primary and secondary bacteria species, and

299

2) bacteria challenge experiments that assess individual and multispecies effects, could shed light

300

on the role of infection-associated bacteria.

301 302

The comparison of infection-associated strains against global MLST databases revealed that P.

303

damselae and V. parahaemolyticus strains founded novel sequence types. The single mutation

304

difference between V. parahaemolyticus ST1931 (described in this study), ST428 (isolated

305

previously from a non-described environmental source in Vietnam), and ST1226 (isolated

306

previously from a non-described environmental source in the Philippines) serves as a reminder

307

that exotic species can be vehicles for the global transport of bacterial pathogens. Biosecurity

308

measures aim to exclude specific pathogens but the working list of pathogens, developed by the

309

United States Marine Shrimp Farming Consortium (USMSFC), is largely focused on the control

310

of viruses (Lightner, 2005). Bacteria belonging to the Vibrionaceae are more difficult to exclude

311

as they are a ubiquitous and natural component of the shrimp’s microflora.

312 313

The search for molecular features that may contribute to competitive interactions between co-

314

occurring species revealed three phylogenetically distinct T6SSs. These systems are highly

315

conserved in many Gram-negative genera and can mediate competition through the delivery of

316

toxic effector proteins (Unterweger et al., 2014; Russell et al., 2014), and T6SS-mediated

317

competition has been characterized as a form of bacteria-bacteria predation that makes available

318

resources for rapid growth and multiplication (Pukatzki and Provenzano, 2013). Here, the

14

319

maintenance of 23 systems across 13 strains suggests that interspecies competition may play a

320

role in shaping infection-associated bacterial communities. The detection of systems homologous

321

to T6SS1 of clinical V. parahaemolyticus RIMD2210633 warrants closer examination of V.

322

alginolyticus and V. harveyi strains, as a previous study detected T6SS1 in toxic strains of V.

323

parahaemolyticus and concluded that possession of this system could impart a competitive

324

advantage over co-occurring species (Li et al., 2017). The detection of an evolutionarily distinct

325

system in V. rotiferianus was novel and also warrants further examination.

326 327

5. Conclusions

328

This study revealed that multiple Harvey clade species were associated with a vibriosis outbreak

329

in Pacific white shrimp (P. vannamei). Further, co-occurring Harveyi clade species possessed

330

numerous protein secretion systems that may contribute to interspecies competition. The

331

identification and characterization of these closely related species was made possible through

332

WGS sequencing, which is increasingly viewed as a high-resolution research methodology in

333

aquaculture.

334 335

Funding

336

This work was supported by Texas A&M University-Corpus Christi.

337 338

Acknowledgments

339

We thank Lee Pinnell for assisting with bioinformatics analyses and manuscript editing. We

340

thank David Prangnell and Tzachi Samocha for providing bacteria cultures, and we thank

341

Patricia Alba and Narjol Gonzalez-Escalona for their expertise and assistance in adding our

15

342

strains to the P. damselae and V. parahaemolyticus databases. We thank Toshio Kodama and

343

Tetsuya Iida and the Pathogenic Microbes Repository Unit of the Research Institute for

344

Microbial Diseases (RIMD) for sharing the V. parahaemolyticus RIMD2210633 type strain.

345 346

Declaration of interest

347

None.

348 349

Appendix A. Supplementary data

350

The following data can be found in the online supplementary material:

351

Supplemental table S1. List of the 43 genomes included in the study.

352

Supplemental table S2. General attributes of the 13 draft genomes sequenced in this study.

353

Supplemental figure S1. Synteny of the T6SS gene clusters detected in the 13 draft genomes

354

sequenced in this study.

355 356

Data availability

357

The whole-genome shotgun projects have been deposited at DDBJ/ENA/Genbank under the

358

accession numbers listed in Table S1. The raw sequence reads were deposited in the Sequence

359

Read Archive under the BioProject accession numbers PRJNA324107, PRJNA324108, and

360

PRJNA446072. Additionally, the P. damselae and V. parahaemolyticus loci sequences and

361

MLST profiles have been deposited in the https://pubmlst.org/pdamselae/ and

362

https://pubmlst.org/vparahaemolyticus/ databases, respectively.

363 364

References

16

365

Alba, P., Caprioli, A., Cocumelli, C., Ianzano, A., Donati, V., Scholl, F., et al., 2016. A new

366

multilocus sequence typing scheme and its application for the characterization of

367

Photobacterium damselae subsp. damselae associated with mortality in cetaceans. Front.

368

Microbiol. 7, 1656. https://doi.org/10.3389/fmicb.2016.01656.

369

Austin, B., & Zhang, X.H., 2006. Vibrio harveyi: a significant pathogen of marine vertebrates

370

and invertebrates. Lett. Appl. Microbiol. 43 (2), 119-124. https://doi.org/10.1111/j.1472-

371

765X.2006.01989.x.

372

Bayliss, S.C., Verner-Jeffreys, D.W., Bartie, K.L., Aanensen, D.M., Sheppard, S.K., Adams, A.,

373

et al., 2017. The promise of whole genome pathogen sequencing for the molecular

374

epidemiology of emerging aquaculture pathogens. Front. Microbiol. 8, 121.

375

https://doi.org/10.3389/fmicb.2017.00121.

376

Bej, A.K., Patterson, D.P., Brasher, C.W., Vickery, M.C., Jones, D.D., Kaysner, C.A., 1999.

377

Detection of total and hemolysin-producing Vibrio parahaemolyticus in shellfish using

378

multiplex PCR amplification of tl, tdh and trh. J. Microbiol. Methods 36 (3), 215-225.

379

https://doi.org/10.1016/S0167-7012(99)00037-8.

380

Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., et al., 2009.

381

BLAST+: architecture and applications. BMC Bioinformatics 10 (1), 421.

382

https://doi.org/10.1186/1471-2105-10-421.

383

Cano-Gomez, A., Høj, L., Owens, L., Andreakis, N., 2011. Multilocus sequence analysis

384

provides basis for fast and reliable identification of Vibrio harveyi-related species and

385

reveals previous misidentification of important marine pathogens. Syst. Appl. Microbiol. 34

386

(8), 561-565. https://doi.org/10.1016/j.syapm.2011.09.001.

17

387

Cano-Gomez, A., Høj, L., Owens, L., Baillie, B. K., Andreakis, N., 2015. A multiplex PCR-

388

based protocol for identification and quantification of Vibrio harveyi-related species.

389

Aquaculture 437, 195-200. https://doi.org/10.1016/j.aquaculture.2014.10.050.

390

Capella-Gutiérrez, S., Silla-Martínez, J.M., Gabaldón, T., 2009. trimAl: a tool for automated

391

alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25 (15), 1972-

392

1973. https://doi.org/10.1093/bioinformatics/btp348.

393

Chatterjee, S., Haldar, S., 2012. Vibrio related diseases in aquaculture and development of rapid

394

and accurate identification methods. J. Mar. Sci. Res. Dev. 1, 2.

395

https://doi.org/10.4172/2155-9910.S1-002.

396

Chen, W.Y., Ng, T.H., Wu, J.H., Chen, J.W., Wang, H.C., 2017. Microbiome dynamics in a

397

shrimp grow-out pond with possible outbreak of acute hepatopancreatic necrosis disease. Sci

398

Rep. 7 (1), 9395. https://doi.org10.1038/s41598-017-09923-6.

399

Cock, J., Salazar, M., Rye, M., 2017. Strategies for managing diseases in non‐native shrimp

400

populations. Rev. Aquacult. 9 (3), 211-226. https://doi.org/10.1111/raq.12132.

401

Colwell, R.R., Grimes, D.J., 1984. Vibrio diseases of marine fish populations. Helgol.

402 403

Meeresunt. 37 (1), 265. https://doi.org/10.1007/BF01989311. Contreras-Moreira, B., Vinuesa, P., 2013. GET_HOMOLOGUES, a versatile software package

404

for scalable and robust microbial pangenome analysis. Appl. Environ. Microbiol., 79 (24),

405

7696-7701. https://doi.org/10.1128/AEM.02411-13.

406

Cornejo-Granados, F., Lopez-Zavala, A.A., Gallardo-Becerra, L., Mendoza-Vargas, A., Sánchez,

407

F., Vichido, R., et al., 2017. Microbiome of Pacific Whiteleg shrimp reveals differential

408

bacterial community composition between wild, aquacultured and AHPND/EMS outbreak

409

conditions. Sci. Rep. 7 (1), 11783. https://doi.org/10.1038/s41598-017-11805-w.

18

410

Dai, W., Yu, W., Xuan, L., Tao, Z., Xiong, J., 2018. Integrating molecular and ecological

411

approaches to identify potential polymicrobial pathogens over a shrimp disease progression.

412

Appl. Microbiol. Biotech. 102 (8), 3755-3764. https://doi.org/10.1007/s00253-018-8891-y.

413 414 415 416 417

Edgar, R.C., 2004. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32 (5), 1792-1797. https://doi.org/10.1093/nar/gkh340. Egan, S., Gardiner, M., 2016. Microbial dysbiosis: rethinking disease in marine ecosystems. Front. Microbiol. 7, 991. https://doi.org/10.3389/fmicb.2016.00991. Englen, M.D., Kelley, L.C., 2000. A rapid DNA isolation procedure for the identification of

418

Campylobacter jejuni by the polymerase chain reaction. Lett. Appl. Microbiol. 31 (6), 421-

419

426. https://doi.org/10.1046/j.1365-2672.2000.00841.x.

420

Giardine, B., Riemer, C., Hardison, R.C., Burhans, R., Elnitski, L., Shah, P., et al., 2005. Galaxy:

421

a platform for interactive large-scale genome analysis. Genome Res. 15 (10), 1451-1455.

422

https://doi:10.1101/gr.4086505.

423

González-Escalona, N., Martinez-Urtaza, J., Romero, J., Espejo, R.T., Jaykus, L.A., DePaola, A.,

424

2008. Determination of molecular phylogenetics of Vibrio parahaemolyticus strains by

425

multilocus sequence typing. J. Bacteriol. 190 (8), 2831-2840.

426

https://doi.org/10.1128/JB.01808-07.

427

Haldar, S., Neogi, S. B., Kogure, K., Chatterjee, S., Chowdhury, N., Hinenoya, A., 2010.

428

Development of a haemolysin gene‐based multiplex PCR for simultaneous detection of

429

Vibrio campbellii, Vibrio harveyi and Vibrio parahaemolyticus. Lett. Appl. Microbiol. 50

430

(2), 146-152. https://doi.org/10.1111/j.1472-765X.2009.02769.x.

19

431

Hoffmann, M., Monday, S. R., Fischer, M., Brown, E.W., 2012. Genetic and phylogenetic

432

evidence for misidentification of Vibrio species within the Harveyi clade. Lett. Appl.

433

Microbiol. 54 (2), 160-165. https://doi.org/10.1111/j.1472-765X.2011.03183.x.

434

Kalyaanamoorthy, S., Minh, B.Q., Wong, T.K., von Haeseler, A., Jermiin, L.S., 2017.

435

ModelFinder: fast model selection for accurate phylogenetic estimates. Nat. Methods 14 (6),

436

587. https:// doi:10.1038/nmeth.4285.

437

Ke, H.M., Prachumwat, A., Yu, C.P., Yang, Y.T., Promsri, S., Liu, K.F., et al., 2017.

438

Comparative genomics of Vibrio campbellii strains and core species of the Vibrio Harveyi

439

clade. Sci. Rep. 7, 41394. https://doi.org/10.1038/srep41394.

440

Klein, S.L., West, C.K.G., Mejia, D.M., Lovell, C.R. 2014. Genes similar to the Vibrio

441

parahaemolyticus virulence-related genes tdh, tlh, and vscC2 occur in other Vibrionaceae

442

species isolated from a pristine estuary. Appl. Environ. Microbiol., 80 (2), 595-602.

443

https://doi:10.1128/AEM.02895-13.

444

Kotob, M.H., Menanteau-Ledouble, S., Kumar, G., Abdelzaher, M., El-Matbouli, M., 2017. The

445

impact of co-infections on fish: a review. Vet. Res. 47 (1), 98. https://doi:10.1186/s13567-

446

016-0383-4.

447

Kumar, G., Engle, C.R., 2016. Technological advances that led to growth of shrimp, salmon, and

448

tilapia farming. Rev. Fish. Sci. Aquac. 24 (2), 136-152.

449

https://doi.org/10.1080/23308249.2015.1112357.

450

Lavilla-Pitogo, C.R., Baticados, M.C. L., Cruz-Lacierda, E.R., Leobert, D., 1990. Occurrence of

451

luminous bacterial disease of Penaeus monodon larvae in the Philippines. Aquaculture 91

452

(1-2), 1-13. https://doi.org/10.1016/0044-8486(90)90173-K.

20

453

Leaño, E.M., Lavilla-Pitogo, C.R., Paner, M.G., 1998. Bacterial flora in the hepatopancreas of

454

pond-reared Penaeus monodon juveniles with luminous vibriosis. Aquaculture 164 (1-4),

455

367-374. https://doi.org/10.1016/S0044-8486(98)00201-4.

456

Lee, C.T., Chen, I.T., Yang, Y.T., Ko, T.P., Huang, Y.T., Huang, J.Y., et al., 2015. The

457

opportunistic marine pathogen Vibrio parahaemolyticus becomes virulent by acquiring a

458

plasmid that expresses a deadly toxin. Proc. Natl. Acad. Sci. USA 112 (34), 10798-10803.

459

https://doi.org/10.1073/pnas.1503129112.

460

Li, P., Kinch, L.N., Ray, A., Dalia, A.B., Cong, Q., Nunan, L.M., et al., 2017. Acute

461

hepatopancreatic necrosis disease-causing Vibrio parahaemolyticus strains maintain an

462

antibacterial type VI secretion system with versatile effector repertoires. Appl. Environ.

463

Microbiol. 83 (13), e00737-17. https://doi.org/10.1128/AEM.00737-17.

464

Li, R., Chiou, J., Chan, E.W.C., Chen, S., 2016. A novel PCR-based approach for accurate

465

identification of Vibrio parahaemolyticus. Front. Microbiol. 7, 44.

466

https://doi.org/10.3389/fmicb.2016.00044.

467

Lightner, D.V., Redman, R.M., 1998. Shrimp diseases and current diagnostic methods.

468

Aquaculture 164 (1-4), 201-220. https://doi.org/10.1016/S0044-8486(98)00187-2.

469

Lightner, D.V., 2005. Biosecurity in shrimp farming: pathogen exclusion through use of SPF

470

stock and routine surveillance. J. World Aquacult. Soc. 36 (3), 229-248.

471

https://doi.org/10.1111/j.1749-7345.2005.tb00328.x.

472

Lin, B., Wang, Z., Malanoski, A.P., O'Grady, E.A., Wimpee, C.F., Vuddhakul, V., et al., 2010.

473

Comparative genomic analyses identify the Vibrio harveyi genome sequenced strains BAA‐

474

1116 and HY01 as Vibrio campbellii. Environ. Microbiol. Rep. 2(1), 81-89.

475

https://doi:10.1111/j.1758-2229.2009.00100.x.

21

476

Liu, C.H., Cheng, W., Hsu, J.P., Chen, J.C., 2004. Vibrio alginolyticus infection in the white

477

shrimp Litopenaeus vannamei confirmed by polymerase chain reaction and 16S rDNA

478

sequencing. Dis. Aquat. Org. 61 (1-2), 169-174. https:// doi:10.3354/dao061169.

479

Makino, K., Oshima, K., Kurokawa, K., Yokoyama, K., Uda, T., Tagomori, K., et al., 2003.

480

Genome sequence of Vibrio parahaemolyticus: a pathogenic mechanism distinct from that of

481

V. cholerae. Lancet 361 (9359), 743-749. https://doi.org/10.1016/S0140-6736(03)12659-1.

482

Martínez-García, P.M., Ramos, C., Rodríguez-Palenzuela, P., 2015. T346Hunter: a novel web-

483

based tool for the prediction of type III, type IV and type VI secretion systems in bacterial

484

genomes. PloS ONE 10 (4), e0119317. https://doi.org/10.1371/journal.pone.0119717.

485 486 487

Minh, B.Q., Nguyen, M.A.T., von Haeseler, A., 2013. Ultrafast approximation for phylogenetic bootstrap. Mol. Biol. Evol. 30 (5), 1188-1195. https://doi:10.1093/molbev/mst024. Moreno, E., Parks, M., Pinnell, L.J., Tallman, J.J., Turner, J.W., 2017. Draft genome sequence of

488

a Vibrio harveyi strain associated with vibriosis in Pacific white shrimp (Litopenaeus

489

vannamei). Genome Announc., 5 (7), e01662-16. https://doi.org/10.1128/genomeA.01662-

490

16.

491 492 493

Msangi, S., Kobayashi, M., Batka, M., Vannuccini, S., Dey, M.M., Anderson, J.L., 2013. Fish to 2030: prospects for fisheries and aquaculture. World Bank Report 83177 (1), 102. Nasu, H., Iida, T., Sugahara, T., Yamaichi, Y., Park, K.S., Yokoyama, K., 2000. A Filamentous

494

Phage Associated with Recent Pandemic Vibrio parahaemolyticus O3: K6 Strains. J. Clin.

495

Microbiol. 38 (6), 2156-2161.

496

Nelson, A., De Soyza, A., Perry, J.D., Sutcliffe, I.C., Cummings, S.P., 2012. Polymicrobial

497

challenges to Koch’s postulates: ecological lessons from the bacterial vaginosis and cystic

498

fibrosis microbiomes. Innate Immun. 18 (5), 774-783.

22

499

https://doi.org/10.1177/1753425912439910.

500

Nguyen, L.T., Schmidt, H.A., von Haeseler, A., Minh, B.Q., 2014. IQ-TREE: a fast and effective

501

stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32

502

(1), 268-274. https://doi:10.1093/molbev/msu300.

503

Park, K.S., Ono, T., Rokuda, M., Jang, M.H., Okada, K., Iida, T., et al., 2004. Functional

504

characterization of two type III secretion systems of Vibrio parahaemolyticus. Infect.

505

Immun. 72 (11), 6659-6665. https://doi:10.1128/IAI.72.11.6659-6665.2004.

506

Prangnell, D.I., Castro, L.F., Ali, A.S., Browdy, C.L., Zimba, P.V., Laramore, S.E., et al., 2016.

507

Some limiting factors in superintensive production of juvenile Pacific white shrimp,

508

Litopenaeus vannamei, in no‐water‐exchange, biofloc‐dominated systems. J. World

509

Aquacult. Soc. 47 (3), 396-413. https://doi.org/10.1111/jwas.12275.

510 511 512

Pukatzki, S., Provenzano, D., 2013. Vibrio cholerae as a predator: lessons from evolutionary principles. Front. Microbiol. 4, 384. https://doi.org/10.3389/fmicb.2013.00384. Richter, M., Rosselló-Móra, R., Oliver Glöckner, F., Peplies, J., 2015. JSpeciesWS: a web server

513

for prokaryotic species circumscription based on pairwise genome comparison.

514

Bioinformatics 32 (6), 929-931. https://doi:10.1093/bioinformatics/btv681.

515

Rivas, A.J., Lemos, M.L., Osorio, C.R., 2013. Photobacterium damselae subsp. damselae, a

516

bacterium pathogenic for marine animals and humans. Front. Microbiol. 4, 283.

517

https://doi.org/10.3389/fmicb.2013.00283.

518

Soto Rodriguez, S.A., Gomez Gil, B., Lozano, R., Roque, A., 2010. Density of vibrios in

519

hemolymph and hepatopancreas of diseased pacific white shrimp, Litopenaeus vannamei,

520

from Northwestern Mexico. J. World Aquacult. Soc. 41, 76-83.

521

https://doi.org/10.1111/j.1749-7345.2009.00335.x.

23

522

Rossi-Tamisier, M., Benamar, S., Raoult, D., Fournier, P.E., 2015. Cautionary tale of using 16S

523

rRNA gene sequence similarity values in identification of human-associated bacterial

524

species. Int. J. Syst. Evol. Microbiol. 65 (6), 1929-1934.

525

https://doi.org/10.1099/ijs.0.000161.

526 527 528

Russell, A.B., Peterson, S.B., Mougous, J.D., 2014. Type VI secretion system effectors: poisons with a purpose. Nat. Rev. Micro. 12 (2), 137. https://doi.org/10.1038/nrmicro3185. Shrivastava, S., Mande, S.S., 2008. Identification and functional characterization of gene

529

components of Type VI Secretion system in bacterial genomes. PloS ONE 3 (8), e2955.

530

https://doi.org/10.1371/journal.pone.0002955.

531

Soto-Rodriguez, S.A., Gomez-Gil, B., Lozano, R., 2010. ‘Bright-red’syndrome in Pacific white

532

shrimp Litopenaeus vannamei is caused by Vibrio harveyi. Dis. Aquat. Org. 92 (1), 11-19.

533

https://doi.org/10.3354/dao02274.

534

Soto-Rodriguez, S.A., Gomez-Gil, B., Lozano-Olvera, R., Betancourt-Lozano, M., Morales-

535

Covarrubias, M. S., 2015. Field and experimental evidence of Vibrio parahaemolyticus as

536

the causative agent of acute hepatopancreatic necrosis disease of cultured shrimp

537

(Litopenaeus vannamei) in Northwestern Mexico. Appl. Environ. Microbiol., 81 (5), 1689-

538

1699. https://doi.org/10.1128/AEM.03610-14.

539

Stamatakis, A., 2014. RAxML version 8: a tool for phylogenetic analysis and post-analysis of

540

large phylogenies. Bioinformatics 30 (9), 1312-1313.

541

https://doi.org/10.1093/bioinformatics/btu033.

542

Tatusova, T., DiCuccio, M., Badretdin, A., Chetvernin, V., Nawrocki, E. P., Zaslavsky, L., et al.,

543

2016. NCBI prokaryotic genome annotation pipeline. Nucleic Acids Res. 44 (14), 6614-

544

6624. https://doi.org/10.1093/nar/gkw569.

24

545

Turner, J.W., Tallman, J.J., Macias, A., Pinnell, L.J., Elledge, N.C., Nasr Azadani, D., et al.,

546

2018. Comparative genomic analysis of Vibrio diabolicus and six taxonomic synonyms: a

547

first look at the distribution and diversity of the expanded species. Front. Microbiol. 9, 1893.

548

https://doi.org/10.3389/fmicb.2018.01893.

549

Unterweger, D., Miyata, S.T., Bachmann, V., Brooks, T.M., Mullins, T., Kostiuk, B., et al.,

550

2014. The Vibrio cholerae type VI secretion system employs diverse effector modules for

551

intraspecific competition. Nat. Commun. 5, 3549. https://doi.org/10.1038/ncomms4549.

552

Urbanczyk, H., Ogura, Y., Hayashi, T., 2013. Taxonomic revision of Harveyi clade bacteria

553

(family Vibrionaceae) based on analysis of whole genome sequences. Int. J. Syst. Evol.

554

Microbiol. 63 (7), 2742-2751. https://doi.org/10.1099/ijs.0.051110-0.

555

Xiong, J., Dai, W., Li, C., 2016. Advances, challenges, and directions in shrimp disease control:

556

the guidelines from an ecological perspective. Appl. Microbiol. Biotechnol. 100 (16), 6947-

557

6954. https://doi.org/10.1007/s00253-016-7679-1.

558

Yao, Z., Yang, K., Huang, L., Huang, X., Qiuqian, L., Wang, K., et al., 2018. Disease outbreak

559

accompanies the dispersive structure of shrimp gut bacterial community with a simple core

560

microbiota. AMB Express 8 (1), 120. https://doi.org/10.1186/s13568-018-0644-x.

561

Yáñez, R., Bastías, R., Higuera, G., Salgado, O., Katharios, P., Romero, J., et al., 2015.

562

Amplification of tlh gene in other Vibrionaceae specie by specie-specific multiplex PCR of

563

Vibrio parahaemolyticus. Electron. J. Biotechnol. 18 (6), 459-463.

564

https://doi.org/10.1016/j.ejbt.2015.09.007.

565

Zhang, X.J., Yan, B.L., Bai, X.S., Bi, K.R., Gao, H., Qin, G.M., 2014. Isolation and

566

characterization of Vibrio parahaemolyticus and Vibrio rotiferianus associated with mass

567

mortality of Chinese shrimp (Fenneropenaeus chinensis). J. Shellfish Res. 33 (1), 61-69.

25

568 569

https://doi.org/10.2983/035.033.0108. Zimin, A.V., Marçais, G., Puiu, D., Roberts, M., Salzberg, S.L., Yorke, J.A., 2013. The

570

MaSuRCA genome assembler. Bioinformatics. 29 (21), 2669-2677.

571

https://doi.org/10.1093/bioinformatics/btt476.

572 573

26

574

Figure legends

575

Figure 1. Pangenome phylogeny. Maximum-likelihood phylogenetic tree based on the

576

alignment of 21 159 binary sites representing the pangenome matrix. Collapsed subclades

577

represent multiple strains. Strains shown in red font were identified/sequenced during this study.

578

Clades highlighted in blue correspond to the five species identified in this study. Branch labels

579

show bootstrap support values. Branch lengths represent the average number of substitutions per

580

site. The tree was rooted to the midpoint.

581 582

Figure 2. T6SS phylogeny. Maximum-likelihood phylogenetic tree based on the multilocus

583

sequence alignment of six T6SS proteins: ClpV, VasA, VasE, VasK, VipA and VipB. Major

584

phylogroups were characterized as T6SS1, T6SS2 and T6SS3 (highlighted in green, purple and

585

blue, respectively), and each of those phylogroups was further divided by two subgroups (A and

586

B). The two systems shown in blue font correspond to T6SS1 and T6SS2 detected in the

587

reference genome V. parahaemolyticus RIMD2210633. Branch labels show bootstrap support

588

values. Branch lengths represent the average number of substitutions per site. The tree was

589

rooted to the midpoint.

590

27

V. harveyi ATCC 43516

V. parahaemolyticus RIMD 2210633

V. rotiferianus HM-10

Strains P. damselae Hep-2a-11 P. damselae Hep-2a-14 P. damselae Hep-2a-16 P. damselae Hep-2b-22 V. alginolyticus Hep-1a-2 V. harveyi Hep-2a-10 V. parahaemolyticus Hep-2a-12 V. parahaemolyticus Hep-2b-20 V. rotiferianus Hep-1a-1 V. rotiferianus Hep-1a-3 V. rotiferianus Hep-2b-18 V. rotiferianus Hep-2b-19 Vibrio sp. Hep-1b-8

V. brasiliensis LMG 20546

seven reference genomes. Values shown in bold font exceed the 96% threshold for species delineation.

V. alginolyticus NBRC 15630

2

P.damselae subsp. piscicida 91-197

Table 1. Pairwise comparison of average nucleotide identity (ANI) values between the 13 bacteria strains sequenced in this study and

P.damselae subsp. damselae KC-Na-1

1

97.26 97.28 97.27 97.30 70.56 70.78 70.30 70.24 70.51 70.46 70.51 70.52 70.65

96.75 96.80 96.79 96.81 70.75 70.96 70.50 70.42 70.74 70.64 70.75 70.60 70.84

70.22 70.25 70.28 70.13 98.26 80.36 83.02 82.98 79.64 79.64 79.65 79.63 73.83

69.75 69.75 69.79 69.78 73.56 74.00 73.84 73.84 74.17 74.17 74.14 74.23 82.38

70.11 70.16 70.21 70.15 80.26 98.50 79.78 79.80 85.89 85.90 85.90 85.54 73.85

70.08 70.16 70.21 70.12 82.92 80.14 98.22 98.22 79.46 79.43 79.47 79.40 73.81

69.96 69.88 69.90 69.92 79.67 85.72 79.27 79.26 96.95 96.95 96.95 96.58 73.65

3 4

1

1

2 3 4

Table 2. Description of the six conserved proteins used for the multilocus sequence analysis of the 23 T6SSs detected in this study. Component Description Length (aa)a ClpV ATPase 851 | 892 VasA Inner membrane protein 582 | 607 VasE Baseplate subunit 345 | 444 VasK Inner membrane protein 1123 | 1181 VipA Contractile sheath, small subunit 148 | 174 VipB Contractile sheath, large subunit 490 | 497 a Minimum and maximum protein length b Average sequence identity in the protein alignment

Identity (%)b 52.34 38.13 38.45 34.23 48.84 53.89

1

1 2

Highlights •

3 4

Whole-genome sequencing (WGS) enabled the identification and characterization of cooccurring, infection-associated Harveyi clade species.

•

Multilocus sequence typing (MLST) revealed that P. damselae and V. parahaemolyticus

5

belonged to novel sequence types (ST56 and ST1931, respectively), although ST1931

6

was closely related to strains isolated from Southeastern Asia.

7

•

A novel multilocus sequence analysis (MLSA) of 23 type VI secretion systems (T6SSs)

8

revealed that co-occurring Harveyi clade species possessed phylogenetically distinct

9

secretion systems.

10

1

1

Conflict of interest

2

The study is original research that has not been previously published and is not under

3

consideration elsewhere. Its publication was approved by all authors, and the authors agree that it

4

will not be published elsewhere without the written consent of the copyright-holder. The authors

5

declare no conflicts of interest.