Comparative Transcriptome Analysis of Two Ipomoea aquatica Forsk

Jun 6, 2016 - To our knowledge, this is the first attempt to use RNA-Seq to detect Cd-induced transcriptional differences between two water spinach cu...
0 downloads 0 Views 1MB Size
Subscriber access provided by UNIV OF NEBRASKA - LINCOLN

Article

Comparative transcriptome analysis of two Ipomoea aquatica Forsk. cultivars targeted to explore possible mechanism of genotype dependent accumulation of cadmium Ying-Ying Huang, Chuang Shen, Jing-Xin Chen, Chun-Tao He, Qian Zhou, Xiao Tan, Jiangang Yuan, and Zhongyi Yang J. Agric. Food Chem., Just Accepted Manuscript • DOI: 10.1021/acs.jafc.6b01267 • Publication Date (Web): 06 Jun 2016 Downloaded from http://pubs.acs.org on June 8, 2016

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.

Journal of Agricultural and Food Chemistry 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 37

Journal of Agricultural and Food Chemistry

1

Comparative transcriptome analysis of two Ipomoea aquatica Forsk. cultivars

2

targeted to explore possible mechanism of genotype dependent accumulation of

3

cadmium

4 5

Ying-Ying Huang, Chuang Shen, Jing-Xin Chen, Chun-Tao He, Qian Zhou, Xiao Tan,

6

Jian-Gang Yuan and Zhong-Yi Yang*

7 8

State Key Laboratory for Biocontrol, School of Life Sciences, Sun Yat-Sen University,

9

Xingang Xi Road 135, Guangzhou, 510275, China

10 11

*

12

Mail address: Xingang Xi Road 135, Guangzhou, 510275, China.

13

Email: [email protected]

14

Tel: +86 2084113220

Corresponding Author

15

ACS Paragon Plus Environment

Journal of Agricultural and Food Chemistry

16

Page 2 of 37

ABSTRACT

17

A low-shoot-Cd (QLQ) and a high-shoot-Cd cultivar (T308) of water spinach

18

(Ipomoea aquatica Forsk.) were used to investigate molecular mechanism of the

19

genotype difference in cadmium (Cd) accumulation. RNA-Seq under 9 h and 72 h

20

cadmium exposures (5 mg L-1) were undertaken to explore Cd induced genotype

21

differences in molecular processes. In total, 253,747,540 clean reads were assembled

22

into 57,524 unigenes. Among them, 6,136 and 10,064 unigenes were differentially

23

expressed in QLQ and T308, respectively. Cell wall biosynthesis genes such as GAUT

24

and laccase and three Cd efflux genes (Nramp5, MATE9 and YSL7) had higher

25

expression levels in QLQ, while the genes in sulfur and glutathione metabolism

26

pathway e.g. sulfate transporter and cysteine synthase showed higher expression

27

levels in T308. These findings would be useful for further understanding of the

28

mechanisms related to genotype-dependent Cd accumulation and developing the

29

molecular assisted screening and breeding of low-shoot-Cd cultivars for water

30

spinach.

31

Low-shoot-Cd

accumulation,

Water

32

KEYWORDS:

33

aquatica Forsk.), Transcriptome, DEGs, Molecular mechanisms

34

ACS Paragon Plus Environment

spinach

(Ipomoea

Page 3 of 37

Journal of Agricultural and Food Chemistry

35

INTRODUCTION

36

Cadmium (Cd) is a non-essential element with high toxicity and can be

37

accumulated in the human body over time from the digestion of food containing Cd,

38

being a threat to human health with excessive intake 1, 2. Cd contamination has been a

39

worldwide environmental problem 3. Therefore, take efficient technologies to reduce

40

Cd content in crops is crucially important. Wide variations among cultivars in the

41

concentration of Cd have been reported in many species, including rice (Oryza sativa)

42

4

43

spinach (Ipomoea aquatic Forsk.) 8. It has been indicated that low-Cd cultivars can be

44

used to reduce the movement of Cd into the human diet 1.

, asparagus bean (Dolichos sinensis) 5, hot pepper (Capsicum annuum) 6, 7 and water

Water spinach is an important leafy vegetable in southern China, it can

45

8, 9

46

accumulate Cd easily when growing in Cd-contaminated soils

. A high-shoot-Cd

47

cultivar (T308) and a low-shoot-Cd cultivar (QLQ) were identified in our previous

48

studies

49

the differences between shoot concentrations of the two cultivars were 3.1~3.4 fold 9.

50

A reciprocal grafting experiment of QLQ and T308 showed that the plants of T308

51

scion with QLQ rootstock accumulated less Cd than the shoot of non-grafted T308

52

and the plants of QLQ scion with T308 rootstock accumulated higher Cd than the

53

shoots of non-grafted QLQ. These observations indicated that the genotype-dependent

54

Cd accumulation of water spinach largely depends on bio-processes occurring in roots

55

11

56

and T308) had been observed by suppression subtractive hybridization (SSH) in our

8, 10

. When growing in Cd contaminated soil (0.593 mg kg-1dry weight, DW),

. Cd-induced gene expression differences of the two water spinach cultivars (QLQ

ACS Paragon Plus Environment

Journal of Agricultural and Food Chemistry

Page 4 of 37

12

57

previous study

. Although there are some studies focus on the molecular

58

mechanisms of Cd hyperaccumulation in several species

59

mechanisms of cultivar dependent difference in Cd accumulation of crops are poorly

60

understood

61

mechanisms of low-shoot-Cd accumulation of water spinach at transcriptional level.

62

For water spinach, which has no reference genome , de novo assembly using

63

high-throughput deep sequencing technology is an effective approach for

64

transcriptome analysis. RNA-Seq has been successfully used in studies of Cd stressed

65

Noccaea caerulescens and Sedum alfredii which has also no reference genome 16, 17.

13, 14

, the molecular

1, 15

. The objective of this study is to uncover the possible molecular

66

Accordingly, RNA-Seq is employed to explore the differences in root

67

transcriptional responses to Cd stress between the two cultivars after short- (9 h) and

68

long-term (72 h) Cd treatment. To our knowledge, this is the first attempt to use

69

RNA-Seq to detect Cd-induced transcriptional differences between two water spinach

70

cultivars with different Cd accumulation abilities. It is expected that the results will

71

help to clarify molecular mechanisms that are responsible for low-shoot-Cd

72

accumulation in water spinach and will be useful to build up breeding methodology of

73

low Cd accumulation cultivars based on molecular assistant breeding techniques.

74

METHODS AND MATERIALS

75

Plant Materials and Treatments

76

Two water spinach cultivars, QLQ (low-shoot-Cd cultivar) and T308

77

(high-shoot-Cd cultivar) were used in this study. The two cultivars were selected from

78

30 cultivars in our previous study, and the differences between shoot Cd

ACS Paragon Plus Environment

Page 5 of 37

Journal of Agricultural and Food Chemistry

79

concentrations of the two cultivars were 3.1~3.4 fold when grew in Cd contaminated

80

soil (0.593 mg kg-1 dry weight, DW) 9. When grew in Hoagland solution containing 6

81

mg L-1 CdCl2 for 24 h, the differences between shoot Cd concentrations was 1.6 fold

82

12

83

. Seeds of the cultivars were surface disinfected with 1‰ carbendazim (m/v)

for

84

0.5 h, then thoroughly washed with deionized water, and sown into pots filled with

85

sand (3 plants per pot) and then were irrigated with a half strength Hoagland nutrient

86

solution for 4 weeks. The plants were cultured at constant temperature condition (25±

87

1°C) with a photoperiod of 14 h (12000 lx). Three Cd exposure lengths were then

88

carried out as treatment by irrigating with the half strength Hoagland solution

89

containing 5 mg L-1 CdCl2 for 0 h, 9 h and 72 h, respectively. The three treatment

90

groups were denoted as E0, E9 and E72, respectively, and each treatment repeated 3

91

times.

92

The shoots and roots of the three plants for each treatment were harvested

93

separately and then washed with deionized water. Thereafter one plant is randomly

94

selected from each pot for further study.

95

Determination of Cd Concentration

96

The roots were washed thorough with tap water and soaked in 25 mM CaCl2

97

twice (5 min each time) to remove heavy metal ions from the root surface. All

98

samples were thoroughly washed with deionized water, dried at 105°C for 30 min,

99

and then dried at 70°C to constant weight. The dried samples were digested (HNO3:

100

H2O2, 5:1) and then analyzed by atomic absorption spectroscopy (Hitachi Z-5300,

ACS Paragon Plus Environment

Journal of Agricultural and Food Chemistry

101

Japan). A certified reference material (CRM) of plant (GBW-07603, provided by the

102

National Research Center for CRM, China) with a Cd concentration of 0.38 mg kg-1

103

was used to control the precision of the analytical procedures. Data statistical analyses

104

including independent samples t-test and one-way ANOVA with the least significant

105

difference (LSD) test were performed using excel 2013 and SPSS 18.0.

106

RNA Extraction

107

Total RNA were extracted separately from the root samples of the three plants as

108

replicates for each treatment with the aid of the RNA Easyspin Isolation System

109

(Aidlab, China) according to the manufacturer’s instruction, and the genomic DNA

110

residues were removed during RNA extraction. RNA concentration and quality were

111

tested in an Agilent 2100 Bioanalyzer (Agilent, USA). RNA quality was also

112

confirmed by RNase free agarose gel electrophoresis. Equal amounts of total RNA

113

extracted from the three replicate plants of each treatment were pooled to construct

114

the cDNA library, so so 6 cDNA libraries (2 cultivars × 3 treatments) were

115

constructed in the present study. Pooling is a cost-effective strategy to identify gene

116

expression profiles and has been well-justified based on statistical and practical

117

considerations 18.

118

cDNA Library Construction and Sequencing

119

Poly (A) mRNA was isolated by using oligo-dT beads (Qiagen, USA) and

120

broken down into short fragments (200 nt) by adding fragmentation buffer.

121

First-strand cDNA was generated using random hexamer-primed reverse transcription

122

and second-strand cDNA was synthesized using RNase H and DNA polymerase I.

ACS Paragon Plus Environment

Page 6 of 37

Page 7 of 37

Journal of Agricultural and Food Chemistry

123

The cDNA fragments were purified with a QIAquick PCR purification kit (Qiagen,

124

USA) and resolved with EB buffer for end reparation, poly (A) addition and ligated to

125

sequencing adapters. After agarose gel electrophoresis, suitable fragments (about 200

126

nt) were used to construct the final cDNA library and was sequenced using Illumina

127

HiSeq™ 2000 (pair-end sequencing, length 125 bp).

128

De Novo Transcriptome Assembly and Annotation

129

Illumina adapter sequences were cut, and low quality reads that did not meet the

130

quality threshold were removed by a Perl program (version 5.18.4). Raw reads

131

containing adapter, reads with more than 10% unknown nucleotides, and low quality

132

reads with over 50% low Q-value (≤5) base were removed as well. The clean reads

133

were de novo assembled using Trinity 2.0 19. The resulting transcripts were annotated

134

with BLASTX (version 2.2.29), a tool widely used to search sequence similarities, by

135

comparing against the database of NCBI non-redundant protein database (Nr,

136

http://www.ncbi.nlm.nih.gov/), Cluster of Orthologous Groups of proteins database

137

(COG, http://www.ncbi.nlm.nih.gov/COG/), Kyoto Encyclopedia of genes and

138

genomes

139

(http://www.expasy.ch/sprot), with an e-value cut off of 10-5.

140

Identification of Differentially Expressed Genes

(KEGG,

http://www.genome.jp/kegg/)

and

Swiss-Prot

141

Clean reads were used to map reference sequence by SOAPaligner/soap2

142

(version 2.21), a tool designed for short sequence alignment. The expression levels of

143

all genes were determined using the RPKM (reads per kb per million reads) method ,

144

according to the formula of RPKM (A)=106C/(NL 10-3). In which, RPKM stands for

ACS Paragon Plus Environment

Journal of Agricultural and Food Chemistry

Page 8 of 37

145

the expression of gene A; C stands for the number of reads uniquely mapped to gene

146

A; N stands for the number of reads uniquely mapped to all genes; L stands for the

147

total length of gene A. For genes with more than one alternative transcript, the longest

148

transcript was used to calculate the RPKM

149

influence of different gene lengths and sequencing discrepancies on the gene

150

expression calculation. The false discovery rate (FDR) was calculated to adjust the

151

threshold of p value using the method described by Benjamini et al

152

between each two samples were considered differentially expressed when the false

153

discovery (FDR) ≤ 0.001 and the absolute value of log2Ratio≥ 1.

154

GO and Pathway Enrichment Analyses of DEGs

20

. The RPKM method can eliminate the

21

. Unigenes

155

To determine the main biological functions of the differentially expressed genes,

156

Gene Ontology (GO, http://www.geneontology.org/) terms were assigned by using

157

Blast2GO (version 2.8) to search the Nr database 22. GO classification was performed

158

in WEGO (http://wego.genomics.org.cn/cgi-bin/wego/index.pl)

159

categorization results were expressed as 3 independent hierarchies for molecular

160

function, biological process and cellular component. The KEGG pathway annotation

161

was performed using BLASTX software against the Kyoto Encyclopedia of Genes

162

and Genomes (KEGG) database according to the method described by Kanehisa et al

163

24

164

correction, and a corrected p-value≤ 0.05 was chosen as the threshold value to

165

determine significantly enriched GO terms. For KEGG enrichment analysis, the

166

pathway with a FDR value≤ 0.05 was considered as an enriched pathway. Pearson

23

, and the GO

. For GO enrichment analysis, all p-values were adjusted with the Bonferroni

ACS Paragon Plus Environment

Page 9 of 37

Journal of Agricultural and Food Chemistry

167

correlation coefficient analysis and principal component analysis (PCA) on all

168

samples were performed by R package (version 3.1.3, http://www.r-project.org/) using

169

the log transformed RPKM values of all genes.

170

Quantitative Real-Time Reverse Transcription PCR

171

The expression levels of 20 genes were determined by quantitative real-time

172

reverse transcription PCR (qRT-PCR). First strand cDNA was generated by purified

173

total RNA using a PrimeScriptTM RT regent kit (Takara, Japan). The water spinach

174

actin (unigene0033718) was selected as a reference gene. All primers for qRT-PCR

175

were shown in Table S1 (suporting information). qRT-PCR was performed on a

176

LightCycler® 480 Real-Time PCR System (Roche, Germany) with SYBR GreenII

177

PCR Master Mix (Takara, Japan). The qPCR cycling program was run as follows:

178

95°C for 3 min, followed by 45 cycles of 95°C for 15 sec and 60°C for 1 min. Each

179

qRT-PCR analysis was performed in triplicate. Results were analyzed with the

180

integrated LightCycler® 480 service software. Expression levels of the tested genes

181

were determined by CT values and calculated by the ∆∆Ct method.

182

RESULTS

183

Cd Concentration of Water Spinach

184

Root dry weight in both two cultivars at E9 and E72 did not show significant

185

differences compared with those at E0 (Figure S1). Shoot and root Cd concentrations

186

of water spinach cultivars QLQ and T308 are shown in Figure 1. Concentrations in

187

both cultivars increased with duration of Cd treatment. T308 always accumulated

188

higher Cd in shoots than QLQ, but only the difference at E72 was significant

ACS Paragon Plus Environment

Journal of Agricultural and Food Chemistry

189

(p0.05).

194

De Novo Sequence Assembly and Functional Annotation

-1

(dry weight basis), respectively, indicating that the genotype

195

RNA-Seq analyses of the two cultivars were performed on the Illumina system to

196

explore the differences that may contribute to their metal-related traits at

197

transcriptional level in roots. As mentioned before that there is no appropriate

198

reference genome sequence available for water spinach, therefore all the 253,747,540

199

clean reads were de novo assembled using Trinity to construct unique consensus

200

sequence. The assembly yielded a total of 57,524 unigenes with an N50 unigene size

201

of 1,627 bp, an average length of 936.74 bp, a minimum length of 201 bp and a

202

maximum length of 17,439 bp. A summary of the unigenes size distribution is given

203

in Figure S2. The raw reads of RNA sequencing were deposited in the Sequence Read

204

Archive (SRA) database of NCBI with the accession number of SRP071620.

205

Only 55.81% (32104) unigenes of all the assembled unigenes (57,524) could be

206

annotated to the 4 databases (Nr, Swiss-Prot, COG and KEGG). GO analysis was

207

performed to classify the function of annotated unigenes. GO analysis categorized

208

16,394 unigenes belonging to cellular component, biological process and molecular

209

function in which cell (8,049), metabolic process (8,734), and catalytic activity (7,981)

210

was the dominant groups, respectively (Figure S3). In addition, 9,166 unigenes were

ACS Paragon Plus Environment

Page 10 of 37

Page 11 of 37

Journal of Agricultural and Food Chemistry

211

mapped to 125 KEGG pathways. The dominant pathway with most unigenes was

212

metabolic pathways, followed by biosynthesis of secondary metabolite, ribosome and

213

plant hormone signal transduction (Table S2).

214

The Pearson correlation between QLQ-72 h and QLQ- 0 h was lower than that

215

between QLQ-9 h and QLQ-0 h, while it was reversed in T308. The correlations

216

between two cultivars decreased with duration of Cd treatment (Table S3). Principal

217

component analysis (PCA) showed that 68.8% of the transcriptional changes

218

associated with Cd stress were explained by the first two PCA components. The first

219

and second component did not cluster the samples by Cd treatment time in both two

220

cultivars, but the transcriptional response to Cd stress of QLQ and T308 could be

221

separated by PCA analyses. For the first principal component, assignment of QLQ-72

222

h was similar to that of T308-9 h (Figure S4).

223

Differential Expressed Genes (DEGs)

224

Unigenes were commonly considered differentially expressed when the FDR ≤

225

0.001 and the absolute value of log2Ratio ≥ 1. RNA-Seq data showed that the gene

226

expression profiles of the both cultivars changed significantly after exposing to 5 mg

227

L-1 CdCl2 at E9 and E72 as compared to E0. The number of DEGs between two

228

cultivars increased over time, and this result was comparable with the increment of

229

the differences in shoot Cd concentrations between two cultivars (Table S4). DEGs

230

between QLQ and T308 at E9 and E72 obtained from Venn analysis are shown in

231

Figure 2. A total of 79 up-regulated and 98 down-regulated unigenes were shared by

232

QLQ and T308 at different treatment time (Table S5).

ACS Paragon Plus Environment

Journal of Agricultural and Food Chemistry

233

Trend Analysis of DEGs

234

To examine the expression profile of the DEGs in QLQ and T308, trend analysis

235

were performed by Short Time-series Expression Miner (STEM) Software. A total of

236

6,136 DEGs in QLQ were clustered into 8 profiles, and only 2 profiles recorded with

237

a significance p