A Reduced Transcriptome Approach to Assess Environmental

Subscriber access provided by READING UNIV

Article

A Reduced Transcriptome Approach to Assess Environmental Toxicants Using Zebrafish Embryo Test Pingping Wang, Pu Xia, Jianghua Yang, Zhihao Wang, Ying Peng, Wei Shi, Daniel.L Villeneuve, Hongxia Yu, and Xiaowei Zhang Environ. Sci. Technol., Just Accepted Manuscript • DOI: 10.1021/acs.est.7b04073 • Publication Date (Web): 11 Dec 2017 Downloaded from http://pubs.acs.org on December 11, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Environmental Science & Technology is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36

Environmental Science & Technology

TOC 84x47mm (300 x 300 DPI)

ACS Paragon Plus Environment


1 2 3

A Reduced Transcriptome Approach to Assess Environmental Toxicants Using Zebrafish Embryo Test

4 5

Pingping Wang†1, Pu Xia†1, Jianghua Yang†, Zhihao Wang†, Ying Peng†, Wei Shi†,

6

Daniel L. Villeneuve ‡，Hongxia Yu†, Xiaowei Zhang†*

7

†State Key Laboratory of Pollution Control & Resource Reuse, School of the

8

Environment, Nanjing University, Nanjing, P. R. China, 210023 ‡

9

United States Environmental Protection Agency, Mid-Continent Ecology Division,

10

Duluth, MN, USA.

11

1

12

*Correspondence:

13

Xiaowei Zhang, PhD, Prof

14

School of the Environment

15

Nanjing University

16

Nanjing, 210089, China;

17

Tel.: 86-25-89680623

18

Fax: 86-25-89680623

19

E-mail: [email protected]

20

[email protected]

, these two authors contributed equally to this paper


Page 2 of 36

Page 3 of 36


21

ABSTRACT

22

Omics approaches can monitor responses and alterations of biological pathways at

23

genome-scale, which are useful to predict potential adverse effects by environmental

24

toxicants. However, high throughput application of transcriptomics in chemical

25

assessment is limited due to the high cost and lack of “standardized” toxicogenomic

26

methods. Here, a reduced zebrafish transcriptome (RZT) approach was developed to

27

represent the whole transcriptome and to profile bioactivity of chemical and

28

environmental mixtures in zebrafish embryo. RZT gene set of 1637 zebrafish Entrez

29

genes was designed to cover a wide range of biological processes, and to faithfully

30

capture gene-level and pathway-level changes by toxicants compared with the whole

31

transcriptome. Concentration-response modeling was used to calculate the effect

32

concentrations (ECs) of DEGs and corresponding molecular pathways. To validate the

33

RZT approach, quantitative analysis of gene expression by RNA-ampliseq technology

34

was used to identify differentially expressed genes (DEGs) at 32 hpf following

35

exposure to seven serial dilutions of reference chemical BPA (10~10E-5µM) or each

36

of four water samples ranging from wastewater to drinking water (relative enrichment

37

factors 10~6.4E-4). The RZT-ampliseq-embryo approach was both sensitive and able

38

to identify a wide spectrum of biological activities associated with BPA exposure.

39

Finally, water quality was benchmarked based on the sensitivity distribution curve of

40

biological pathways detected using RZT-ampliseq-embryo, and the most sensitive

41

biological pathways were identified, including those linked with adverse reproductive

42

outcomes, genotoxicity and development outcomes. RZT-ampliseq-embryo approach



43

provides an efficient and cost-effective tool to prioritize toxicants based on

44

responsiveness of biological pathways.

45 46


Page 4 of 36

Page 5 of 36


47

1. INTRODUCTION

48

A major challenge with regard to prioritizing environmental chemicals and/or

49

assessing the hazard of complex mixtures is the lack of sufficient toxicological

50

information for thousands of chemicals and endless possibility of mixtures. Toxicity

51

pathway profiling could help to predict potential apical toxicity and prioritize and

52

guide subsequent testing of the chemicals.1 Traditionally, monitoring and assessment

53

of mixtures have relied on chemistry analyses. Although high throughput targeted and

54

non-targeted analytical methods have been developed for detection of hundreds of

55

chemicals present in complex environmental samples, chemical-focused analyses

56

cannot detect contaminants with unknown structure, and cannot explain the

57

cumulative toxicity of mixtures.2 Effect-based approaches, such as high content

58

screening, can provide assessments of biological activity of environmental mixture.2

59

However, most current cell-based HTS assays are limited in their coverage of

60

biological pathways, and consequently their ability to predict a wide range of potential

61

adverse outcomes.

62

Our previous development on the reduced human transcriptome (RHT), which

63

integrated a panel of 1200 core toxicologically relevant genes and dose-response

64

modeling, has been shown to be an efficient and cost-effective approach to benchmark

65

the water quality from wastewater to drinking water using human cell-based tests.3 To

66

overcome the inherent limitations of single cell type and broaden the coverage of

67

biological pathways in assessing toxicants, here a reduced transcriptome approach

68

was developed using the zebrafish embryo as a multicellular test system. Zebrafish is



69

an excellent experimental model to study chemical induced toxicity, because of its

70

genetic similarity to humans and other vertebrates as well as its well characterized

71

genome. Zebrafish embryo in particular has received considerable attention as an

72

efficient alternative model that can be used to monitor altered molecular response

73

induced by environmental stimuli, and linked with apical endpoints such as

74

developmental and neuro toxicities. However, the utility of zebrafish model was

75

significantly constrained in omics-based toxicological research due to the high cost

76

and the lack of standardized bioinformatics protocols for dose-response modelling.

77

Integrating genomic dose-response modeling into the hazard characterization with

78

wide-range doses, particularly in lower, environmentally-relevant doses, has been

79

seen to be valuable in risk assessment4. Genomic studies on a wide-range doses could

80

help to determine new biomarkers and to derive point of departure for chemical risk

81

assessment. For instance, certain endocrine disrupting chemicals have been reported

82

to alter gene expression in a non-monotonic manner at low dose, which indicate

83

potential novel molecular mechanism.5 In addition, application of multiple doses with

84

single replicate using zebrafish embryo has been shown to effectively identify

85

androgen-responsive genes.6 Concentration-dependent bioactivity of water samples

86

could indicate potential early responses.7 Pathway analysis based on the active values

87

of differentially concentration-dependent genes implicate the potential bioactivity of

88

samples, which can be used in diagnostic analysis of chemical profiles.3 However,

89

utilization of biological pathway responses derived from dose-dependent genomic

90

data is still limited in hazard characterization.6, 8, 9


Page 6 of 36

Page 7 of 36


91

The objectives of this study were three fold. The first was to curate a reduced gene

92

list from zebrafish transcriptome (RZT) that can comprehensively represent biological

93

pathways and toxicologically relevant processes, and be quantified by Ion Ampliseq

94

Technology (RZT-ampliseq) (Figure 1). Second, we aimed to develop a chemical test

95

protocol integrating RZT-ampliseq and dose-response modeling in zebrafish embryo

96

(RZT-ampliseq-embryo). Bisphenol A (BPA), a well-studied endocrine disruptor

97

frequently detected in water samples, was selected as a reference chemical. Finally,

98

we wanted to evaluate the performance of RZT-ampliseq-embryo for use in hazard

99

assessment of environmental mixtures. The mixture samples tested in this study were

100

a set of water extracts which have been previously characterized by a battery of in

101

vitro assays10 and reduced human transcriptome (RHT) method3.

102 103

2. MATERIALS &METHODS

104

2.1. Design of a gene set for reduced zebrafish transcriptome. The RZT gene set

105

was selected to represent the key biological pathways and toxicologically relevant

106

processes in zebrafish (Danio rerio) genome (Figure 1).

107

associated with key biological pathways (in Entrez ID formats) was curated from

108

three databases, including Kyoto Encyclopedia of Genes and Genomes (KEGG),11

109

zebrafish orthologs of L1000 landmark genes12 and zebrafish orthologs of pathway

110

reporter genes (Table 1).13 The centrality values of genes were calculated using

111

CentiScaPe in Cytoscape software.14 Centrality values are node parameters

112

demonstrating the relevant position of nodes in a whole network. Higher centrality


First, a list of genes


Page 8 of 36

113

value suggests more central roles of a gene in biological pathways. Then the numbers

114

of significantly enriched KEGG pathways and GO terms (adjusted p-value 0, Reads > 20) and the coefficient of

204

variation (CV) (SD/Mean) of each gene’s expression abundance. For each gene, the

205

relationship between CV and sequencing depth was fitted with loess model and then

206

the minimum sequencing depth of ensuring CV < 15 % was calculated. To evaluate

207

the mRNA profiling performance of ampliseq on the RZT gene set, the number of

208

detected and undetected genes, as well as the each gene expression abundance

209

measured by RZT-ampliseq were compared to that by microarray platform

210

(GSE43186) and RNA-seq30 platform on 36 hpf zebrafish embryo. Correlations of the

211

gene expression abundance between different technologies were calculated using

212

number of reads per amplicon for RZT-ampliseq, RPKM (reads per kilo base per

213

million reads) values for RNA-seq and signal intensity values for microarray

214

technology. Finally, to evaluate the repeatability of RZT-ampliseq, the CV of RZT

215

gene set in zebrafish embryos of 0.1% DMSO (n=3) from six batches were analyzed

216

using the edge package.27

217

Pathway-level and biological process validation. For single dose experiment,

218

functional enrichment analysis of identified DEGs was performed using a one-sided

219

Fisher’s exact test on GO of Biological Process (BP), and KEGG pathways with RZT

220

gene list (Table S5) as background. For full dose experiment, the EC values of GO


Page 12 of 36

Page 13 of 36


221

terms and KEGG pathways were calculated as the geometric mean of EC values of

222

matched DEGs. Only GO terms or KEGG pathways matched by at least three genes

223

were included in EC calculation and further analysis. Finally, to analyze overall

224

biological potency of each sample, the proportionally ranked distribution of GO and

225

KEGG of EC values was fitted with a four-parameter dose-response curve using

226

GraphPad Prism 5.0 software (San Diego, CA, U.S.A.).

227

The molecular responses profiling (DEGs, KEGG pathways of DEGs) of 0.1 µM

228

BPA treatment by RZT-ampliseq were compared with whole transcriptome analysis

229

of BPA archived in NCBI.31 To compare RZT-ampliseq-embryo approach with

230

existing Toxcast high throughput in vitro assays with regard to biological activities

231

associated with BPA exposure, the responsive gene endpoints and molecular

232

pathways (KEGG, GO BP terms) identified by the both methods were evaluated. The

233

responsive molecular gene endpoints were DEGs captured by dose-response model

234

analysis of RZT- ampliseq-embryo. The responsive genes of Toxcast in vitro assay

235

were

236

(https://www.epa.gov/chemical-research/toxicity-forecaster-toxcasttm-data).

237

responsive molecular endpoints were converted to zebrafish orthologous genes.

download

from The

238

2.5 Comparison of RZT with in vitro bioassays & RHT method on mixtures. A

239

supervised approach was used to assess the RZT representation of the previous in

240

vitro bioassays. First, gene sets associated with cellular toxicity pathways tested by in

241

vitro bioassays were manually curated from Wiki Pathways and Gene Ontology,

242

KEGG (Table S6). Then the EC of each pathway was calculated by the geometric



243

mean of the ECs of matched DEGs. Pathway patterns identified by the RZT approach

244

was shown by heatmap using gplot package.32 The hierarchical clusters of water

245

samples identified by RZT analysis were compared with the results of in vitro

246

bioassays.

247

To evaluate the sensitivity and specificity of RZT-ampliseq-embryo in

248

identification of bioactivity of mixtures, the results of RZT-ampliseq-embryo were

249

compared to that by RHT-ampliseq using human HepG2 and MCF7 cells3 on the

250

same sample set. Briefly, the sensitivity of 50% biological potency of water samples

251

identified by RZT were compared with those identified by RHT in HepG2 and MCF7

252

cells in terms of KEGG or GO. In addition, linear regression was conducted on values

253

of 50% biological potency identified by RZT and RHT. Finally, the coverage of most

254

sensitive pathways (top 20 sensitive KEGG pathways) of Eff2, the sample with

255

potential highest and broadest bioactivity were compared between RZT and RHT

256

approaches.

257

3. RESULTS AND DISCUSSION

258

3.1. Development of RZT testing method using zebrafish embryo

259

In silico validation of RZT gene set. The developed RZT gene set consists of 1637

260

zebrafish Entrez ID genes, including a list of 1000 genes with greatest pathway

261

centrality scores and a list of 724 toxicology-relevant genes. The 1000

262

pathway-central genes were shown to be the minimum number of genes representing

263

the maximum biological pathways in terms of GO BP terms and KEGG pathways

264

(Figure 2A). Toxicology-relevant genes (n=724) were selected to provide linkages


Page 14 of 36

Page 15 of 36


265

between molecular mechanism and apical endpoints (Table 1). Then 44 genes were

266

removed by the online designer either because their background expression was too

267

high or too low, or because effective multiplexed primers could not be designed

268

(Table S7). This resulted in 1637 genes as the final RZT gene set (Table S5).

269

The RZT gene set showed a broad coverage of biological pathways, where 95%

270

KEGG pathways and 94% GO BP terms were represented by at least one gene in RZT

271

gene set (Figure 2B and C). The uncovered pathways were mainly associated with

272

basic metabolic processes (Table S8). Furthermore, the RZT gene set of 1637 genes

273

were significantly enriched in 29 KEGG pathways and 839 GO BP terms (adjusted

274

p0) were detected in all 32hpf embryo RNA samples (n=18)

310

under sequencing depth of one million (Figure S5A). Of the 179 genes not detected

311

(count=0) in 32 hpf embryo by RZT-ampliseq, 47 genes were consistently undetected

312

in 36 hpf embryo in a previous RNA-seq study

313

detected by RNA-seq instead of RZT-ampliseq platform, 86% (93 out of 108) had

314

normalized reads less or equal than 10 counts (RPKM), among which 72% (79 out of

315

108 genes) had a sequence read below 1 count. Nine genes weren’t detected by

316

RNA-seq but were detected by RZT-ampliseq platform, among which 5 genes had a

317

normalized read below 10. The mRNA expression abundances quantified by

318

RZT-ampliseq and RNA-seq platform showed a nearly linear relation (R= 0.81)

319

between two platforms for the 1270 gene transcripts detected by both (counts>0)

320

(Figure S5B). Moreover, out of 82 transcripts that were not detected by RZT-ampliseq

321

but were detected by microarray, only 29 were above 8 on log2 scale (Figure S5C, D).

322

A similar linear relationship (R=0.77) was observed for the 1254 common genes

323

between the RZT-ampliseq approach and a microarray platform (Figure S5D).

30

. Among the other 108 genes

324

3.2. RZT assessment of a classical chemical: BPA.

325

The RZT approach showed good repeatability for quantifying transcriptional

326

response to chemical by zebrafish embryo. Common CV of 32 hpf embryo mRNA

327

samples exposed to 0.1% DMSO from 8-32 hpf were 13% (biological replication

328

within

329

RZT-ampliseq-embryo (Figure S6). This variation was acceptably low when

330

compared with other RNA profiling technology, such as qPCR (CV: 1% ~ 15%),

one

batch),

14%

(biological

replication


between

6

batch)

for


331

microarray (CV:5% ~ 15%) or RNA-seq (CV:10% ~ 15%).33 After exposure of two

332

independent batches of embryos to a single dose of 10 µM BPA, 67, 45 DEGs

333

(ANOVA, p=3) were identified in two platforms, but these

349

were only involved in fundamental apoptosis process (FoxO signaling pathway, p53

350

signaling pathway) and regulation of actin cytoskeleton (Figure S8). However, 98

351

DEGs were identified by dose-response analysis of the embryo exposed to the seven,

352

10-fold, dilutions of BPA (10 E-5 ~10 µM). Three and five of the DEGs identified by


Page 18 of 36

Page 19 of 36


353

dose-response analysis were also detected as DEGs in embryos exposed to 10 and 0.1

354

µ M, respectively (Figure S7A).

355

One significant advantage of dose-response analysis by RZT-ampliseq is the

356

sensitivity analysis of genes and biological pathways in response to chemicals, which

357

could aid inference regarding the potentially sensitive apical endpoint effects. The

358

responsive DEGs were mainly fitted with U-shaped models, which suggest that the

359

mode of hormesis dominates the low dose response of transcriptome (Figure S7C).

360

However, there are alternative interpretations other than true hormesis. For example,

361

the time-course of dynamic transcriptional response may be different at different

362

doses, such that at higher doses, the transcript abundance may have peaked earlier but

363

has fallen by 32 h.

364

higher doses, associated with triggering more and more AOPs in the organism,

365

thereby causing more and more disruption of normal development. What is effectively

366

a monotonic response to the chemical may produce a non-monotonic dose-response

367

for a given snap-shot in time. The response genes (ECRW >DW). The EC values of the most sensitive

415

KEGG or GO pathways of DW samples were 1-2 orders of magnitude higher than

416

those of effluent samples, suggesting relatively weak biological effects were induced

417

by DW. (Figure S14)



418

The enriched pathways in RZT analysis could be used to prioritize potential

419

biological endpoints for future assessment. Specifically, the most sensitive KEGG or

420

GO BP pathways may be linked with adverse outcome. For instance, the top 20%

421

sensitive KEGG or GO BP pathways of Eff2 were mainly associated with adverse

422

reproductive outcomes, genotoxicity and development outcomes (Figure S10). The

423

most sensitive RZT profile of MF samples was associated with the adverse outcome

424

of endocrine system disturbance and development. For DW, the most sensitive

425

responses were also associated with genotoxicity, which might be due to the fact that

426

the samples of metropolitan drinking water treatment plant contained toxic

427

“byproduct” chemicals from chlorination9. All development relevant pathways (GO

428

terms, each covered at least 3 DEGs suggested Eff2 and MF samples might induce

429

potential development toxicity while RW and DW samples may not (Figure 4A). The

430

predicted adverse outcomes were corroborated by zebrafish embryo 48hpf lethality

431

and 120hpf sub-lethal development experiment 9(Figure 4A).

432

The RZT-ampliseq-embryo method provided more sensitive pathway profiles for

433

the four water samples comparable to the previous 103 in vitro bioassays. (Figure

434

S11) First, the RZT profiling analysis revealed similar pattern although with greater

435

sensitivity than the in vitro bioassays. Similar to the in vitro cell-based bioassays,

436

RZT-ampliseq-embryo approach could clearly distinguish wastewater (Eff2, MF)

437

from reclaimed or clean water samples (DW, RW). Biological responses associated

438

with genotoxicity and oxidative stress responses, xenobiotic metabolism, PR,

439

RAR/RXR were identified using both RZT-ampliseq-embryo and in vitro bioassays.


Page 22 of 36

Page 23 of 36


440

However, the RZT method identified weak TR and hypoxia-related responses (RPE:

441

0-10), which were not detected in the in vitro assays. The sensitive detection on the

442

thyroid hormone signaling pathways by the RZT-ampliseq-embryo method might be

443

contributed by the active involvement of TR pathway during zebrafish embryo

444

development. However, the RZT profile showed low sensitivity for ER and GR,

445

which comes as no surprise, because ER and GR would show high sensitivity at 48hpf

446

or later phase zebrafish embryo37.

447

The RZT approach provided a broader biological coverage than in vitro cell-based

448

bioassays which focused on preselected biomarkers assessment. Thus, the RZT

449

approach could support more comprehensive assessment of biological effects of

450

environmental samples. Beside the responsive endpoints in vitro bioassays, the RZT

451

approach also identify other biological responses such as development and

452

reproduction-related pathways, some of which might be linked to AOPs of regulatory

453

concerns. For example, profiling analysis on the enriched KEGG pathways (Figure

454

S12) indicated that Eff2, MF samples other than RW, DW samples exhibited potential

455

development toxicity by inducing several relevant pathway responses, including

456

VEGF signaling pathway, Hedgehog signaling pathway and Vascular smooth muscle

457

contraction. An established AOPs “disruption of VEGFR signaling leading to

458

developmental defects” ( https://aopwiki.org/wiki/index.php/Aop:43) might be used to

459

guide the prediction of developmental risk by the observed molecular event.

460

Additionally, Hedgehog pathway was verified to play an essential role in

461

cardiomyocyte differentiation and heart morphogenesis in model species including



462

mouse38 and zebrafish39. Furthermore, wastewater Eff2 and MF samples showed

463

prominent pathway effects involved in the reproductive axis including oocyte meiosis,

464

TGF-beta signaling pathway, progesterone-mediated oocyte maturation, GnRH

465

signaling pathway, ErbB signaling pathway. However, inlet and outlet samples of

466

metropolitan drinking water treatment plant (RW, DW) showed very weak

467

reproduction-related effects.

468 469

3.4 Comparison between RZT-embryo and RHT cells profiles of water samples.

470

Although only four common water samples were tested by RZT and RHT

471

approaches, the water quality of four water samples ranked by 50% biological potency

472

identified by RZT-embryo were consistent with those by RHT in MCF7 and HepG2

473

cells (Figure S13,S14). The values of 50% biological potency by RZT-embryo

474

correlated well with those by RHT in MCF7 (R2=0.95) in terms of GO (Figure S13),

475

and correlated well with those by RHT in HepG2 (R2=0.95) in terms of KEGG

476

(R2=0.99) (Figure S14).

477

RZT-embryo also provided different profiles of altered genes & pathways of the

478

four water samples from that by RHT approach, which might be due to the greater

479

biological complexity represented by a fish embryo, compared to a single cell type.

480

The most sensitive pathways identified by RZT following exposure to the water

481

samples was distinct with those by RHT in HepG2 and MCF7. Take Eff2 for example

482

(Figure 4B), only four KEGG pathways were overlapped between the 20 most

483

sensitive pathways (with lowest EC values) identified by RZT-embryo and RHT in


Page 24 of 36

Page 25 of 36


484

HepG2 and MCF7 cells. Nine of the 20 most sensitive pathways uniquely identified

485

by RZT-embryo were associated with basic biological processes, which may suggest

486

that rapidly developing and differentiating zebrafish embryos were more sensitive to

487

alterations of basic processes, such as oxidative phosphorylation, than the single

488

cell-type in vitro system. Moreover, the most sensitive KEGG pathways identified by

489

RHT in HepG2 and RHT in MCF7 showed cell-type responses, such as pathways

490

involved in immune response and cellular communication, which were not among in

491

the most sensitive KEGG pathways by RZT-embryo assay. However, some cell-type

492

specific responses, including endocrine response in MCF7 and metabolism response

493

in HepG2, were also identified by RZT as sensitive KEGG pathways responding to

494

Eff2.

495

The RZT-embryo approach was more sensitive than RHT in HepG2 and MCF7

496

cells in terms of estimating 50% biological potency of water samples. The value of 50%

497

biological potency of each sample identified by RZT was generally one order of

498

magnitude less than that identified by RHT in either HepG2 or MCF7 (Figure S13,14,

499

Table S3A, B). Furthermore, a distinct pathway sensitivity distribution in response to

500

the MF sample was identified by RZT-embryo compared to RHT in HepG2 and

501

MCF7 cells. Although the potency of the median sensitive pathway (50% biological

502

potency) following exposure to MF was lower than that of RW and DW in

503

RZT-embryo, MF was more potent than that of RW and DW at the most sensitive

504

pathways which were profiled by the RZT-embryo. (Figure 13, 14) These highly

505

sensitive biological responses induced by MF were primarily related to embryo



(eg.

heart

jogging,

embryo

pattern

specification,

Page 26 of 36

506

development

notochord

507

development, determination of left/right symmetry) (Figure 4A), which might be

508

related to developmentally toxic pollutants present in the MF sample. The MF sample

509

was water taken after microfiltration using filters disinfected by chlorination to avoid

510

biofouling in a water reclamation plant, in which micro-pollutants of such as

511

carbamazepine, a teratogen, with the highest detected concentrations (1.9 µg / L) out

512

of ten sample24 were present. Carbamazepine has been reported to disturb embryonic

513

development with increasing hatching rate, body length, swim bladder appearance and

514

yolk sac absorption rate at 1µg/L.40 However, knowledge gaps associated with

515

unknown chemicals present in the mixtures and their potential combined effects still

516

exist in toxicological assessment of these environmental samples. An effect-directed

517

analysis (EDA) integrating extract fractionation and instrument analysis with the

518

sensitive RZT approach may be used to identify the chemicals responsible for the

519

observed effect in future.

520

In conclusion, we developed a RZT approach by integrating reduced transcriptome,

521

RNA-ampliseq technology and a zebrafish embryo test as a novel approach to assess

522

environmental toxicant. First, the RZT approach by multiple dose-response analysis

523

could identify early molecular response and molecular mechanism of single chemical

524

which would help to predict apical effect. Additionally, the RZT approach has

525

potential to be used to evaluate and prioritize chemicals for further testing and

526

potentially to predict adverse outcomes. Second, RZT-ampliseq-embryo approach

527

effectively discriminated relatively clean from more polluted environmental samples,


Page 27 of 36


528

and the responsive pathways revealed by RZT analysis could help to prioritize

529

targeted biological end points for future assessment. Future studies should explore the

530

utility of other stages, such as 48hpf (sensitivity for endocrine disrupting effects), in

531

the identification of early molecular response for assessing the hazard potencies of

532

environmental toxicants in early developmental stage.

533

Acknowledgements. For support, we thank National Natural Science Foundation of

534

China (Grant No. 21322704), Environmental Protection Foundation of Jiangsu

535

(ZX2015009) and the European Union Seventh Framework Programme (The

536

SOLUTIONS project, grant 603437). P.X. was supported by Program B for

537

Outstanding Ph.D. Candidates of Nanjing University (No. 201701B018). X.Z. was

538

supported by the Fundamental Research Funds for the Central Universities. Thanks

539

Professor Beate Escher and Eutox, UQ for provision of the water samples and for

540

helpful discussion. Mention of trade names or commercial products does not

541

constitute endorsement or recommendation for use by the U.S. Environmental

542

Protection Agency.

543 544

Supporting Information

545 546 547

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.est.xxxxx. Description of methods for estimating transcriptional point of departure (PODt); S2,

548

CV change with different reads count of all genes; fourteen figures and eight tables

549

(PDF, XLSX)

550 551 552 553 ACS Paragon Plus Environment


554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597

REFERENCE 1.

Kavlock, R.; Chandler, K.; Houck, K.; Hunter, S.; Judson, R.; Kleinstreuer, N.; Knudsen, T.; Martin,

M.; Padilla, S.; Reif, D.; Richard, A.; Rotroff, D.; Sipes, N.; Dix, D., Update on EPA's ToxCast Program: Providing High Throughput Decision Support Tools for Chemical Risk Management. Chem. Res. Toxicol. 2012, 25, (7), 1287-1302. 2.

Altenburger, R.; Ait-Aissa, S.; Antczak, P.; Backhaus, T.; Barcelo, D.; Seiler, T.-B.; Brion, F.; Busch,

W.; Chipman, K.; Lopez de Alda, M.; de Aragao Umbuzeiro, G.; Escher, B. I.; Falciani, F.; Faust, M.; Focks, A.; Hilscherova, K.; Hollender, J.; Hollert, H.; Jaeger, F.; Jahnke, A.; Kortenkamp, A.; Krauss, M.; Lemkine, G. F.; Munthe, J.; Neumann, S.; Schymanski, E. L.; Scrimshaw, M.; Segner, H.; Slobodnik, J.; Smedes, F.; Kughathas, S.; Teodorovic, I.; Tindall, A. J.; Tollefsen, K. E.; Walz, K.-H.; Williams, T. D.; Van den Brink, P. J.; van Gils, J.; Vrana, B.; Zhang, X.; Brack, W., Future water quality monitoring - Adapting tools to deal with mixtures of pollutants in water resource management. Sci. Total Environ. 2015, 512, 540-551. 3.

Xia, P.; Zhang, X.; Zhang, H.; Wang, P.; Tian, M.; Yu, H., Benchmarking Water Quality from

Wastewater to Drinking Waters Using Reduced Transcriptome of Human Cells. Environ. Sci. Technol. 2017, 51, (16), 9318-9326. 4.

Beausoleil, C.; Ormsby, J.-N.; Gies, A.; Hass, U.; Heindel, J. J.; Holmer, M. L.; Nielsen, P. J.; Munn,

S.; Schoenfelder, G., Low dose effects and non-monotonic dose responses for endocrine active chemicals: Science to practice workshop: Workshop summary. Chemosphere 2013, 93, (6), 847-856. 5.

Faulk, C.; Kim, J. H.; Jones, T. R.; McEachin, R. C.; Nahar, M. S.; Dolinoy, D. C.; Sartor, M. A.,

Bisphenol A-associated alterations in genome-wide DNA methylation and gene expression patterns reveal sequence-dependent and non-monotonic effects in human fetal liver. Environ. Epigenet. 2015, 1, (1). 6.

Fetter, E.; Smetanova, S.; Baldauf, L.; Lidzba, A.; Altenburger, R.; Schuettler, A.; Scholz, S.,

Identification and Characterization of Androgen-Responsive Genes in Zebrafish Embryos. Environ. Sci. Technol. 2015, 49, (19), 11789-11798.


Page 28 of 36

Page 29 of 36


598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640

7.

Berninger, J. P.; Martinovic-Weigelt, D.; Garcia-Reyero, N.; Escalon, L.; Perkins, E. J.; Ankley, G. T.;

Villeneuve, D. L., Using Transcriptomic Tools to Evaluate Biological Effects Across Effluent Gradients at a Diverse Set of Study Sites in Minnesota, USA. Environ. Sci. Technol. 2014, 48, (4), 2404-2412. 8.

Thomas, R. S.; Allen, B. C.; Nong, A.; Yang, L.; Bermudez, E.; Clewell, H. J., III; Andersen, M. E., A

method to integrate benchmark dose estimates with genomic data to assess the functional effects of chemical exposure. Toxicol. Sci. 2007, 98, (1), 240-248. 9.

Smetanova, S.; Riedl, J.; Zitzkat, D.; Altenburger, R.; Busch, W., High-throughput

concentration-response analysis for omics datasets. Environ. Toxicol. Chem. 2015, 34, (9), 2167-2180. 10. Escher, B. I.; Allinson, M.; Altenburger, R.; Bain, P. A.; Balaguer, P.; Busch, W.; Crago, J.; Denslow, N. D.; Dopp, E.; Hilscherova, K.; Humpage, A. R.; Kumar, A.; Grimaldi, M.; Jayasinghe, B. S.; Jarosova, B.; Jia, A.; Makarov, S.; Maruya, K. A.; Medvedev, A.; Mehinto, A. C.; Mendez, J. E.; Poulsen, A.; Prochazka, E.; Richard, J.; Schifferli, A.; Schlenk, D.; Scholz, S.; Shiraish, F.; Snyder, S.; Su, G.; Tang, J. Y. M.; van der Burg, B.; van der Linden, S. C.; Werner, I.; Westerheide, S. D.; Wong, C. K. C.; Yang, M.; Yeung, B. H. Y.; Zhang, X.; Leusch, F. D. L., Benchmarking Organic Micropollutants in Wastewater, Recycled Water and Drinking Water with In Vitro Bioassays. Environ. Sci. Technol. 2014, 48, (3), 1940-1956. 11. Collins, F. S.; Gray, G. M.; Bucher, J. R., Toxicology. Transforming environmental health protection. Science 2008, 319, (5865), 906-7. 12. Duan Q, R., SP, Clark NR, Wang Z, Fernandez NF, Rouillard AD, Readhead B, Tritsch SR, Hodos R, Hafner M, Niepel M, Sorger PK, Dudley JT, Bavari S, Panchal RG, Ma’ayan A, L1000CDS2: LINCS L1000 characteristic direction signatures search engine. NPJ. Syst. Biol. Appl. 2016, 2, (16015). 13. Zhang, J. D.; Küng, E.; Boess, F.; Certa, U.; Ebeling, M., Pathway reporter genes define molecular phenotypes of human cells. BMC Genomics 2015, 16, (1), 342. 14. Scardoni, G.; Petterlini, M.; Laudanna, C., Analyzing biological network parameters with CentiScaPe. Bioinformatics 2009, 25, (21), 2857. 15. Yu, G.; Wang, L. G.; Han, Y.; He, Q. Y., clusterProfiler: an R package for comparing biological themes among gene clusters. Omics : a journal of integrative biology 2012, 16, (5), 284-7. 16. Villeneuve, D. L.; Garcia-Reyero, N.; Martinović-Weigelt, D.; Li, Z.; Watanabe, K. H.; Orlando, E. F.; Lalone, C. A.; Edwards, S. W.; Burgoon, L. D.; Denslow, N. D., A graphical systems model and tissue-specific functional gene sets to aid transcriptomic analysis of chemical impacts on the female teleost reproductive axis. Mutat. Res. 2011, 746, (2), 151-162. 17. Jiang, J.; Wu, S.; Wu, C.; An, X.; Cai, L.; Zhao, X., Embryonic exposure to carbendazim induces the transcription of genes related to apoptosis, immunotoxicity and endocrine disruption in zebrafish (Danio rerio). Fish Shellfish Immunol. 2014, 41, (2), 493-500. 18. Li, Y.; Qi, X.; Yang, Y.-W.; Pan, Y.; Bian, H.-M., Toxic effects of strychnine and strychnine N-oxide on zebrafish embryos. Chin. J. Nat. Med. 2014, 12, (10), 760-767. 19. Guiu, J.; Bergen, D. J. M.; De Pater, E.; Islam, A. B. M. M. K.; Ayllon, V.; Gama-Norton, L.; Ruiz-Herguido, C.; Gonzalez, J.; Lopez-Bigas, N.; Menendez, P.; Dzierzak, E.; Espinosa, L.; Bigas, A., Identification of Cdca7 as a novel Notch transcriptional target involved in hematopoietic stem cell emergence. J. Exp. Med. 2014, 211, (12), 2411-2423. 20. Verleyen, D.; Luyten, F. P.; Tylzanowski, P., Orphan G-Protein Coupled Receptor 22 (Gpr22) Regulates Cilia Length and Structure in the Zebrafish Kupffer's Vesicle. PLoS One 2014, 9, (10) :e110484.



641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683

21. Wanglar, C.; Takahashi, J.; Yabe, T.; Takada, S., Tbx Protein Level Critical for Clock-Mediated Somite Positioning Is Regulated through Interaction between Tbx and Ripply. PLoS One 2014, 9, (9) :e107928. 22. Xu, M.; Liu, D.; Dong, Z.; Wang, X.; Wang, X.; Liu, Y.; Baas, P. W.; Liu, M., Kinesin-12 Influences Axonal Growth During Zebrafish Neural Development. Cytoskeleton 2014, 71, (10), 555-563. 23. Xu, M.; Liu, D.; Dong, Z.; Wang, X.; Wang, X.; Liu, Y.; Baas, P. W.; Liu, M., Kinesin-12 influences axonal growth during zebrafish neural development. Cytoskeleton 2014, 71, (10), 555-63. 24. Bluhm, K.; Otte, J. C.; Yang, L.; Zinsmeister, C.; Legradi, J.; Keiter, S.; Kosmehl, T.; Braunbeck, T.; Straehle, U.; Hollert, H., Impacts of Different Exposure Scenarios on Transcript Abundances in Danio rerio Embryos when Investigating the Toxicological Burden of Riverine Sediments. PLoS One 2014, 9, (9) :e106523. 25. Robinson, M. D.; McCarthy, D. J.; Smyth, G. K., edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, (1), 139-40. 26. Farmahin, R.; Williams, A.; Kuo, B.; Chepelev, N. L.; Thomas, R. S.; Barton-Maclaren, T. S.; Curran, I. H.; Nong, A.; Wade, M. G.; Yauk, C. L., Recommended approaches in the application of toxicogenomics to derive points of departure for chemical risk assessment. Arch. Toxicol. 2017, 91, (5), 2045-2065. 27. Robinson, M. D.; McCarthy, D. J.; Smyth, G. K., edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, (1), 139-140. 28. Ritz, C.; Streibig, J. C.; Leeuw, J. D.; Zeileis, A., Bioassay Analysis Using R. J. STAT. SOFTW. 2005, 12, (i05), 1-22. 29. Bornkamp, B.; Pinheiro, J.; Bretz, F., DoseFinding: Planning and Analyzing Dose Finding experiments. 2016. 30. Yang, H.; Zhou, Y.; Gu, J.; Xie, S.; Xu, Y.; Zhu, G.; Wang, L.; Huang, J.; Ma, H.; Yao, J., Deep mRNA Sequencing Analysis to Capture the Transcriptome Landscape of Zebrafish Embryos and Larvae. PLoS One 2013, 8, (6) :e64058. 31. Saili, K. S.; Tilton, S. C.; Waters, K. M.; Tanguay, R. L., Global gene expression analysis reveals pathway differences between teratogenic and non-teratogenic exposure concentrations of bisphenol A and 17 beta-estradiol in embryonic zebrafish. Reprod. Toxicol. 2013, 38, 89-101. 32. Warnes, G. R.; Bolker, B.; Bonebakker, L.; Gentleman, R.; Huber, W.; Liaw, A.; Lumley, T.; Maechler, M.; Magnusson, A.; Moeller, S., gplots: various R programming tools for plotting data. 2016. 33. Zhang, J. D.; Schindler, T.; Kueng, E.; Ebeling, M.; Certa, U., Highly sensitive amplicon-based transcript quantification by semiconductor sequencing. BMC Genomics 2014, 15 :565. 34. Kinch, C. D.; Ibhazehiebo, K.; Jeong, J.-H.; Habibi, H. R.; Kurrasch, D. M., Low-dose exposure to bisphenol A and replacement bisphenol S induces precocious hypothalamic neurogenesis in embryonic zebrafish. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, (5), 1475-1480. 35. Wang, C.; Zhang, J.; Li, Q.; Zhang, T.; Deng, Z.; Lian, J.; Jia, D.; Li, R.; Zheng, T.; Ding, X.; Yang, F.; Ma, C.; Wang, R.; Zhang, W.; Wen, J. G., Low concentration of BPA induces mice spermatocytes apoptosis via GPR30. Oncotarget 2017, 8(30):49005-49015. 36. Wu, M.; Pan, C.; Chen, Z.; Jiang, L.; Lei, P.; Yang, M., Bioconcentration pattern and induced apoptosis of bisphenol A in zebrafish embryos at environmentally relevant concentrations. Environ. Sci. Pollut. Res. Int. 2017, 24, (7), 6611-6621.


Page 30 of 36

Page 31 of 36


684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708

37. Schiller, V.; Wichmann, A.; Kriehuber, R.; Schaefers, C.; Fischer, R.; Fenske, M., Transcriptome alterations in zebrafish embryos after exposure to environmental estrogens and anti-androgens can reveal endocrine disruption. Reprod. Toxicol. 2013, 42, 210-223. 38. Zhang, X. M.; Ramalho-Santos, M.; McMahon, A. P., Smoothened mutants reveal redundant roles for Shh and Ihh signaling including regulation of L/R asymmetry by the mouse node. Cell 2001, 105, (6), 781-792. 39. Hami, D.; Grimes, A. C.; Tsai, H.-J.; Kirby, M. L., Zebrafish cardiac development requires a conserved secondary heart field. Development 2011, 138, (11), 2389-2398. 40. Qiang, L.; Cheng, J.; Yi, J.; Rotchell, J. M.; Zhu, X.; Zhou, J., Environmental concentration of carbamazepine accelerates fish embryonic development and disturbs larvae behavior. Ecotoxicology 2016, 25, (7), 1426-1437.

Table 1. Sources and selection criteria of reduced zebrafish transcriptome (RZT) gene panel. 1000 core genes were selected by their central roles in the network of genes collected from public databases. A final list was union of the 1000 core genes and toxicology-relevant genes. Category of genes

Public databases Core genes

Toxicology-relevant genes RZT gene panel

Number of zebrafish orthologs 4260 1022 1019 1000 326 173 176 152 1637

Sources KEGG database18 L1000 landmark genes19 Pathway reporter genes20 Selected by their central roles in the gene network ToxCast AOP wiki Graphical gene model23 Retrieved from references24-31 Core genes + toxicology-relevant genes

709 710 711



Page 32 of 36

712 713 714 715 716 717 718 719 720 721 722 723 724 725 726

Table 2. Treatments and exposure design for the validation of reduced zebrafish transcriptome (RZT) approach using embryo test with exposure time of 8-32hpf. Type

Treatment

Solvent

DMSO

Single Chemical

Exposure level

Dimethyl Sulfoxide

0.1% DMSO (6 batches, n=3)

BPA

Bisphenol A

Eff2

Secondary sewage effluent Microfiltration treated effluent River water

MF Mixture

Description

25

RW DW a

10~10E-5 (µM) with 10-fold dilution, 0.1, 10 µM (n=3) b, 10 µM (n=3 of second batch)

10~6.4E-4 (REFa) with 5-fold dilution

Drinking water b

REF: relative enrichment factor. If not otherwise stated, all treatment groups were with a single replicate in addition to three vehicle controls (0.1% DMSO)

727


Page 33 of 36


Figure 1. Design and workflow of the reduced zebrafish transcriptome (RZT) approach. RZT-ampliseq, RZT genes were quantified by Ion Ampliseq Technology; RZT-ampliseq-embryo, a chemical test protocol integrating RZT-ampliseq and dose-response modeling in zebrafish embryo; EC, effect concentration; DEGs, differentially expressed genes. 83x156mm (300 x 300 DPI)



Figure 2. Selection and in silico validation of RZT gene set list. (A) Investigation of minimum number of candidate genes for representing the maximal biological pathways including KEGG pathways and GO BP terms. The number of significantly enriched biological pathways was according to the number of candidate genes ranked by their centrality scores. The curves were fitted into Gaussian model. The red dash line means the cut-off of 1000, where the number of top ranked candidate genes may be low enough for representing maximal biological pathways. The percentage of biological pathways coverage of (B) KEGG pathways and (C) GO BP terms by 1637 genes from RZT gene set. (D, E) Comparison on the point of departure (PODt) estimated by RZT gene and the whole transcriptome in previously published studies (EMTAB-832 and GSE55618, respectively). The distributions of PODt estimated by ten approaches using RZT gene set (blue) and whole zebrafish genome (yellow) were showed in boxplot. The black bold lines represent PODt. The number represents the ratio of PODt between by RZT gene set and whole genome (larger value to smaller value). In plot (D), the solid lines represent LOAEL (13.5 µM) for pericardial edema (green line), and LOAEL (28 µM) for malformed heart (red line) induced by flusilazole in zebrafish embryo at 24 hpf. The


Page 34 of 36

Page 35 of 36


dashed lines represent 3-fold ranges of corresponding LOAELs. In plot (E), the solid and dashed red lines represent 1/3 and 1/10 values of LOAEL (8 mM) for liver damage induced by isoniazid in zebrafish embryo reported by Zhang12.

83x139mm (300 x 300 DPI)



Figure 3. Concentration-dependent network of differentially expressed genes (DEGs) (p=3 were included in this analysis. Common molecular endpoints were labeled by red triangle. Pathway scores (EC or AC50 value) were the geometric mean of the effect concentrations (ECs) (RZT-ampliseq-embryo) or AC50 (Toxcast) values of the relevant genes. The horizontal distance of 50% biological potency between RZT-ampliseq-embryo and Toxcast in vitro assays was labeled in red.

72x260mm (300 x 300 DPI)


Page 36 of 36

Page 37 of 36




Figure 4. (A) Development-related biological processes affected by four water samples from RZT-ampliseqembryo in 32 hpf. Plotted were log 10 EC values with a unit of REF (relative enrichment factor). (B) Venn diagram of top 20 sensitive KEGG pathways ranked by EC values identified by RZT, RHT in HepG2 and RHT in MCF7, respectively, for Eff2 sample. The labels of “zebrafish” in blue, “HepG2” in yellow and “MCF7” in green stand for approaches of RZT, RHT in HepG2 and RHT in MCF7, respectively. The labels in red stand for the main function of KEGG pathways. 84x196mm (300 x 300 DPI)


Page 38 of 36

A Reduced Transcriptome Approach to Assess Environmental

Recommend Documents