Transcriptome Analysis Reveals Potential Antioxidant Defense

Subscriber access provided by Stony Brook University | University Libraries

Food and Beverage Chemistry/Biochemistry

Transcriptome Analysis Reveals Potential Antioxidant Defense Mechanisms in Antheraea pernyi in Response to Zinc Stress Yu Liu, Zhao-Zhe Xin, Jiao Song, Xiao-Yu Zhu, Qiuning Liu, Daizhen Zhang, Bo-Ping Tang, Chun-Lin Zhou, and Li-Shang Dai J. Agric. Food Chem., Just Accepted Manuscript • DOI: 10.1021/acs.jafc.8b01645 • Publication Date (Web): 05 Jul 2018 Downloaded from http://pubs.acs.org on July 6, 2018

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 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 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.

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

Journal of Agricultural and Food Chemistry

1

Transcriptome Analysis Reveals Potential Antioxidant Defense Mechanisms in

2

Antheraea pernyi in Response to Zinc Stress

3

Yu Liua,b,d,#, Zhao-Zhe Xina,d,#, Jiao Songc,#, Xiao-Yu Zhua, Qiu-Ning Liua,*, Dai-Zhen

4

Zhanga, Bo-Ping Tanga,*, Chun-Lin Zhoua, Li-Shang Daib,*

5 6

a

7

Innovation Center for Coastal Bio-agriculture, Jiangsu Provincial Key Laboratory of

8

Coastal Wetland Bioresources and Environmental Protection, School of Ocean and

9

Biological Engineering, Yancheng Teachers University, Yancheng 224051, People's

Jiangsu Key Laboratory for Bioresources of Saline Soils, Jiangsu Synthetic

10

Republic of China

11

b

12

People's Republic of China

13

c

14

Republic of China

15

d

16

Technology, Nanjing 210009, People's Republic of China

17

* Corresponding author: Bo-Ping Tang, Li-Shang Dai, and Qiu-Ning Liu

18

Tel/fax: +86 515 88233991

19

E-mail: [email protected] (BP Tang), [email protected] (LS Dai),

20

[email protected] (QN Liu).

21

Address: Jiangsu Key Laboratory for Bioresources of Saline Soils, Jiangsu Synthetic

22

Innovation Center for Coastal Bio-agriculture, Jiangsu Provincial Key Laboratory of

23

Coastal Wetland Bioresources and Environmental Protection, School of Ocean and

24

Biological Engineering, Yancheng Teachers University, Yancheng 224051, People's

25

Republic of China.

26

#

School of Pharmaceutical Sciences, Wenzhou Medical University, Wenzhou 325035,

College of Life Science, Anhui Agricultural University, Hefei 230036, People's

College of Biotechnology and Pharmaceutical Engineering, Nanjing University of

These authors contributed equally to this work.

27 28

ACS Paragon Plus Environment


29

ABSTRACT

30

The growth and development of the Chinese oak silkworm, Antheraea pernyi (A.

31

pernyi), is strongly influenced by environmental conditions, including heavy metal

32

pollution. An excess of heavy metals causes cellular damage through the production

33

of free-radicals reactive oxygen species. In this study, transcriptome analysis was

34

performed to investigate global gene expressions when A. pernyi was exposed to zinc

35

infection. With RNA sequencing (RNA-Seq), a total of 25,795,510 and 38,158,855

36

clean reads were obtained from zinc-treated and control fat bodies libraries,

37

respectively. We identified 2,399 differentially expression genes (DEGs) (1,845

38

up-regulated and 544 down-regulated genes) in zinc-treated library. Further, Kyoto

39

Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed that these

40

DEGs were related to peroxisome pathway that was associated with antioxidant

41

defense. Our results suggest that fat bodies of A. pernyi constitute a strong antioxidant

42

defense against heavy metal contamination.

43

KEYWORDS: Antheraea pernyi; antioxidant mechanisms; heavy metals; zinc;

44

transcriptome

45

INTRODUCTION

46

A. pernyi, one of the most important wild silkworm species in economy, produce

47

thousands of tons of tussah silk per year and thus play a vital role in the textile

48

industry.1,2 A. pernyi diverged from dipterans 240 million years ago and was

49

domesticated in China in the 18th Century.3,4 To date, its breeding ground has covered

50

more than 700,000 acres, mainly distributing in China, Japan, India, and southeast


Page 2 of 36

Page 3 of 36


51

Asian countries.5 It has been reported that 90% of the world’s tussah cocoons and

52

silks originate from China, which provides the suitable habitat for silkworm survival.6

53

A. pernyi is an completely metamorphic insect with a diapause state and its life

54

history can be divided into four stages: egg, larva, pupa, and moth. The pupae of A.

55

pernyi contains all the essential amino acids and are considered to be a source of

56

high-quality dietary protein.7 Furthermore, A. pernyi is now used as a model organism

57

for a variety of research purposes because of its ease of rearing and experimental

58

manipulation compared with other lepidopteran insects. However, the health and

59

vigor of A. pernyi are threatened by numerous parasites and pathogens, including

60

viruses,8 mites,9 bacteria10 and protoza.11 In addition, the growth and development of

61

larva are also strongly affected by environmental conditions, including heavy metal

62

pollution, climate, foliage and diseases.12

63

Heavy metals are natural compounds found in soil, rock, air, water and living

64

organisms and are essential micronutrients for growth and development of organisms,

65

such as copper and zinc. Heavy metals, however, are also serious pollutants and have

66

led to a worldwide environmental problem.13,14 These substances pose serious threats

67

to human health and global ecosystems via various ways, including biogeochemical

68

cycle, microbial activity and human activities.15,16 An excess of heavy metals in cells

69

cause oxidative stress because of the formation of free radicals and reactive oxygen

70

species (ROS).17,18 ROS can be produced in all aerobic organisms due to the

71

mandatory dependence of oxidative metabolism on the ground-state of molecular

72

oxygen.19 Insects have evolved an array of antioxidant defense systems to keep



Page 4 of 36

73

themselves from oxidative damages.20 Superoxide dismutase (SOD), catalase (CAT),

74

glutathione peroxide (GSH-PX) and glutathione reductase (GR) play vital roles in

75

protecting organisms from endogenous ROS, and they can clean ROS to avoid

76

oxidative damage.19 In addition, glutathione S-transferases (GSTs) are involved in the

77

metabolism

78

antioxidant.19,21 These antioxidant enzymes are also thought to be essential

79

components of insect antioxidant defense systems.22 In recent years, there have been

80

several reports of antioxidant enzymes present in different species of Lepidopteran

81

insect larvae. 23-25 However, to our knowledge, antioxidant enzymes involved in the

82

response to zinc challenge have not been reported. Zinc, the second most abundant

83

metal in organisms, has been identified as a contaminant of concern in Japan and

84

selected as a model toxicant.26

of

polycyclic

hydrocarbons,

pesticides

and

the

endogenous

85

Transcriptome sequencing is widely used for genome analysis and functional

86

gene identification, aiding our understanding of genetic responses of hosts to heavy

87

metals and the molecular mechanisms of antioxidant defense systems. In the present

88

study, transcriptome analysis was conducted on A. pernyi fat bodies using an Illumina

89

sequencing platform to identify functional genes. As an important organ, insect fat

90

body is equivalent to the liver and adipose of vertebrates, and plays key roles in

91

metabolic activities including nutrient supply and intermediary metabolism.27-29 The

92

differentially expressed genes (DEGs) were identified in different functional

93

databases and the results were annotated to investigate their potential functions. The

94

findings provide a further understanding of antioxidant mechanisms in A. pernyi.


Page 5 of 36


95

MATERIALS AND METHODS

96

Insect collection and challenge experiments

97

Individuals of A. pernyi were collected from the Sericultural Research Institute of

98

Henan and reared on oak leaves until pupation. The pupae were kept at room

99

temperature in advance, and then were divided into two groups (n = 3 per group): the

100

zinc exposure group and the control group. The former was treated with 10 µL ZnCl2

101

solution while the latter was treated with 10 µL phosphate buffered saline (PBS)

102

solution. After 24 h infection, the fat bodies of these two groups were collected,

103

labeled and then stored in liquid nitrogen for RNA extraction.

104

Total RNA isolation, cDNA library construction, and sequencing

105

Library construction and RNA-Seq were performed at Beijing BioMarker

106

Technologies (Beijing, China). Total RNA from fat bodies tissue was isolated with

107

TRIzol Reagent (Invitrogen, USA) according to the manufacturer’s instructions. After

108

extraction, RNA degradation and contamination were monitored on 1% agarose gels.

109

In addition, RNA purity, concentration and integrity were evaluated with NanoDrop

110

(Implen, Westlake Village, CA, U.S.A.), Qubit 2.0 (Life Technologies, Carlsbad, CA,

111

U.S.A.), and Agilent 2100 (Agilent Technologies, Santa Clara, CA, U.S.A)

112

instruments, respectively.

113

Briefly, mRNA was purified from total RNA using poly-T oligo-attached

114

magnetic beads, and mRNA fragmentation was carried out using divalent cations

115

under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×).

116

First-strand cDNA was synthesized using random hexamer primers and M-MuLV



117

reverse transcriptase (RNase H-). Subsequently, second-strand cDNA was synthesized

118

using DNA polymerase I and RNase H. These double-stranded cDNA fragments

119

underwent end repair, addition of a single ‘A’ base, and ligation of adapters, and

120

adaptor-modified fragments were isolated by gel purification and amplified by PCR to

121

create the final cDNA library.

122

Transcriptome Sequencing of the cDNA libraries derived from the zinc exposure

123

and control groups was performed on Illumina HiSeq 2500 platform, which is based

124

on synthesis technology. Before proceeding with the analysis, it is necessary to

125

confirm that the quality of the reads was sufficient to ensure the accuracy of the

126

sequence assembly and subsequent analysis. The de novo assembly of RNA-Seq was

127

performed using Trinity software.30

128

Functional unigenes annotation and classification

129

The functional unigenes were annotated based on the following databases: NCBI

130

non-redundant protein/nucleotide sequences (Nr/Nt, http://www.ncbi.nlm.nih.gov/, E

131

values cutoff ≤ 1e-5),31 Protein family (Pfam, https://pfam.sanger.ac.uk/, E values

132

cutoff ≤ 1e-5),32 euKaryotic Ortholog Groups/Clusters of Orthologous Groups of

133

proteins (KOG/COG, https://www.ncbi.nlm.nih.gov/COG/, E values cutoff 1e-3),33

134

Swiss-Prot (http://www.ebi.ac.uk/uniprot, E values cutoff ≤ 1e-5),34 Kyoto

135

Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/, E values

136

cutoff ≤ 1e-10),35 and Gene Ontology (GO, http://www.geneontology.org/, E values

137

cutoff ≤ 1e-6). The Blast2GO program36 was used to assign GO annotations to three

138

main categories (molecular function, cellular component and biological process)


Page 6 of 36

Page 7 of 36


139

based on the Nr annotation (E values cutoff ≤ 1e-6). Then, GO functional classification

140

was performed with WEGO software37 and the KEGG pathways of unigenes were

141

analyzed based on BLASTX hits against the KEGG database.

142

Quantification of gene expression levels and DEG screening

143

The gene expression levels were estimated by RNA-Seq by Expectation

144

Maximization (RSEM) for each sample.38 Clean data were mapped back onto the

145

assembled transcriptome using TopHat software (version 2.1.1), and the read count

146

for each gene was obtained from the mapping results. Calculation of unigene

147

expression was measured by fragment per kilobase of exon per million fragments

148

mapped (FPKM) based on the number of uniquely mapped reads.39 The longest

149

transcript was selected to calculate the FPKM when a gene has more than one

150

alternative transcripts.40 Differential expression analysis of two groups was performed

151

using the DESeq R package to screen out the DEGs. To correct for multiple testing,

152

the false discovery rate (FDR) was calculated to adjust the threshold of p-value. DEGs

153

were defined as genes with FDR values < 0.01, and |log2(foldchange)| > 1 (minimal

154

2-fold difference in expression).41

155

GO and KEGG enrichment analyses

156

GO enrichment analysis of DEGs was implemented using the topGO R packages

157

based on the Kolmogorov–Smirnov test.42 KEGG, as a database resource, was used to

158

interpret high-level functions and utilities of the cell, the organism and the ecosystem

159

based on molecular-level information.43 Furthermore, KOBAS software was used to

160

test the statistical significance of DEG enriched in KEGG pathways.44



161

RESULTS AND DISCUSSION

162

Transcriptome sequence and assembly

163

The transcriptome of A. pernyi was assembled de novo using paired-end raw reads

164

generated by the Illumina HiSeq2500. After filtering out redundant and short reads,

165

25,795,510 clean reads were obtained in the control (PBS-treated) group and

166

38,158,855 in the zinc exposure group (Table 1). The Illumina guidelines were

167

applied for data quality assessment and the result showed that more than 92% of the

168

sequencing data for each sample were found to have a quality score of Q30. The GC

169

counts for the control group and the zinc exposure group were 43.12% and 43.07%,

170

respectively. Besides, 20,423,257 (79.17%) and 29,020,615 (76.05%) clean reads

171

were respectively obtained from two groups, which successfully matched to the A.

172

pernyi genome. According to the statistics, 65.80% of the reads were uniquely

173

mapped to the genome for the control group and 71.38% for the zinc exposure group.

174

Furthermore, 16.89% and 28.62% of the reads, respectively, were multiply mapped to

175

the genome for the control group and the zinc exposure group. Most unigenes are

176

distributed at 200-300 bp, 300-500 bp, and 500-1000 bp (Fig. 1). These results

177

demonstrated that the sequencing quality was relatively high, which means

178

the unigenes were suitable for subsequent annotation analysis.

179

Functional annotation and classification

180

To analyze the gene function information, unigenes were annotated using seven

181

databases (Nr, Swiss-Prot, COG, KOG, eggNOG, GO and KEGG). In total, 22,620 of

182

unigenes were matched in Nr. In Nr annotation, 22,543 unigenes were matched to


Page 8 of 36

Page 9 of 36


183

multiple species genomes, including Bombyx mori (39.48%), Danaus plexippus

184

(16.92%), Tribolium castaneum (11.29%), Ceratitis capitata (1.66%), Papilio xuthus

185

(1.30%), Dendroctonus ponderosae (1.20%), Manduca sexta (0.84%), Acanthamoeba

186

castellanii (0.73%), Acyrthosiphon pisum (0.71%), Nematostella vectensis (0.65%),

187

and others (25.02%) (Table 2).

188

GO classification is a standardized system for categorizing genes and gene

189

products across species and it was performed using the Blast2GO program. On the

190

basis of the functional annotation, 12,602 unigenes were classified into three main GO

191

categories: biological process, cellular component and molecular function. These

192

unigenes were further divided into 57 subcategories as shown in Figure 2 and Table 3.

193

Among these subcategories, 20, 19, and 18 of them were clustered in the biological

194

process, cellular component and molecular function categories, respectively. In the

195

biological process category, the most abundant groups were “metabolic process” (GO:

196

0008152) with 7,861 unigenes and “cellular process” (GO: 0009987) with 6,377

197

unigenes. Within the cellular component category, most of the unigenes were assigned

198

to the terms “cell” (GO 0005623) and “cell part” (GO 0044464), which involved in

199

the basic structural and functional unit of all organisms.45 However, among sequences

200

categorized as molecular function, 6,542 unigenes were predicted to have “catalytic

201

activity” (GO:0003824) and 6,412 unigenes were predicted to have “binding” (GO:

202

0005488) activity. These results indicate that the annotated unigenes were most often

203

associated with various types of biological process. According to the GO analysis,

204

immune-related genes were enriched in the biological process categories, with many



205

unigenes assigned to the subcategories “immune system process” (GO:0002376) and

206

“response to stimulus” (GO:0050896). A total of 1,475 and 138 unigenes were

207

assigned to the subcategories “immune system process” and “response to stimulus”,

208

respectively. In addition, 70 unigenes were assigned to “antioxidant activity”

209

(GO:0016209).

210

The KOG database is used to evaluate the completeness of transcriptomes and the

211

effectiveness of the annotation processes. In this study, 14,961 unigenes were

212

annotated and classified into 25 KOG categories (Fig. 3 and Table 4). Among the

213

functional classes, the category with the greatest number of unigenes was “general

214

function prediction only” (R, 3,115, 20.82%), followed by “signal transduction

215

mechanisms” (T, 1,692, 11.31%), and “posttranslational modification, protein

216

turnover, chaperones” (O, 1,517, 10.14%). The three categories “cell motility”,

217

“nuclear structure”, and “extracellular structures” contained the lowest number of

218

unigenes, accounting for 36, 69, and 124 unigenes, respectively. In total, 169 unigenes

219

were assigned to the cluster of “defense mechanisms”, indicating that these unigenes

220

might be involved in antioxidant defense in A. pernyi. Additionally, the functions of

221

890 unigenes were unknown.

222

Identification and analysis of DEGs

223

To identify genes in A. pernyi whether or not expressed differently in response to zinc

224

exposure, the numbers of clean reads from the transcriptomes of the control and zinc

225

exposure groups were mapped back, the read counts compared, and the gene

226

expression levels estimated. Pearson correlation analysis revealed a significant


Page 10 of 36

Page 11 of 36


227

difference between the two groups (r2 = 0.4813; Fig. 4). Differential expression

228

analysis revealed 2,399 significant DEGs (1,845 upregulated and 554 downregulated)

229

in the zinc exposure group compared with the control group, with FDR 1 (Fig. 5).

231

GO and KEGG enrichment analysis of DEGs

232

To further investigate the function of DEGs, unigenes were annotated using different

233

functional databases: 1,004 in COG, 655 in GO, 851 in KEGG, 1,211 in KOG, 1,336

234

in Pfam, 968 in Swiss-Prot, 1,462 in eggNOG and 1,397 in Nr. According to GO

235

enrichment analysis using Blast2GO with a correct p-value

Transcriptome Analysis Reveals Potential Antioxidant Defense

Recommend Documents