Effect of dietary energy restriction and subsequent compensatory feeding on testicular transcriptome in developing rams

Accepted Manuscript Effect of dietary energy restriction and subsequent compensatory feeding on testicular transcriptome in developing rams Y.X. Fan, Z. Wang, C.F. Ren, T.W. Ma, K.P. Deng, X. Feng, F.Z. Li, F. Wang, Y.L. Zhang PII:

S0093-691X(18)30421-7

DOI:

10.1016/j.theriogenology.2018.06.028

Reference:

THE 14610

To appear in:

Theriogenology

Received Date: 14 November 2017 Revised Date:

28 June 2018

Accepted Date: 28 June 2018

Please cite this article as: Fan YX, Wang Z, Ren CF, Ma TW, Deng KP, Feng X, Li FZ, Wang F, Zhang YL, Effect of dietary energy restriction and subsequent compensatory feeding on testicular transcriptome in developing rams, Theriogenology (2018), doi: 10.1016/j.theriogenology.2018.06.028. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

revised

ACCEPTED MANUSCRIPT 1

Effect of dietary energy restriction and subsequent compensatory feeding on testicular transcriptome in developing rams

2 3

Y.X. Fan, Z. Wang, C.F. Ren, T.W. Ma, K.P. Deng, X. Feng, F.Z. Li, F. Wang, and Y.L.

5

Zhang*

6

Jiangsu Livestock Embryo Engineering Laboratory, College of Animal Science and

7

Technology, Nanjing Agricultural University, Nanjing 210095, PR China.

9

Yanli Zhang

SC

*Corresponding author:

M AN U

8

RI PT

4

Jiangsu Livestock Embryo Engineering Laboratory

11

College of Animal Science and Technology

12

Nanjing Agricultural University

13

No.1 Weigang

14

Nanjing

15

China

16

Tel: 86-25-84395381; Fax: 86-25-84395314; E-mail address: [email protected]

18 19

EP

AC C

17

TE D

10

20 21 22 1

ACCEPTED MANUSCRIPT ABSTRACT

24

Nutritional intake and reproductive allocation are strongly associated and dietary energy

25

restriction (ER) or surpluses can affect reproductive capacity. The objective of this study was

26

to investigate the effect of energy levels on sheep testicular development. Three-month old Hu

27

sheep were assigned to four groups, and fed diets containing different levels of energy

28

(Control, maintenance energy; ER1, 85% maintenance energy; ER2, 70% maintenance energy;

29

ER3, 55% maintenance energy). Two months later, half the sheep in each group were

30

euthanized, whereas the remaining sheep were euthanized after a further 3 months feeding on

31

a compensatory energy diet. The testicular weight and reproductive hormone levels of the Hu

32

sheep were investigated. Differences in the testes of ER3 and control group sheep were

33

investigated at the transcriptional level using high-throughput sequencing. The results showed

34

that the testicular weights had decreased in the energy-restricted rams compared with the

35

controls, and that the testosterone concentration in ER3 group rams was significantly lower

36

than that in other compared groups (P < 0.05). After the period of compensatory feeding,

37

however, ER3 sheep testicular weight and testosterone concentrations were similar to those of

38

the control group sheep. In addition, the RNA sequencing results revealed that 81 genes were

39

upregulated and 180 genes were downregulated in the ER3 group compared with the control

40

group. Moreover, based on the enriched steroidogenesis, meiosis and kinases pathways, a

41

number of candidate genes potentially involved in the regulation of testicular development or

42

reproduction of Hu sheep, including CYP11A1, ALDH3B1, FDFT1, WNT2, PGR and INSR,

43

were screened. Quantitative real-time polymerase chain reaction analysis results correlated

44

well with the sequencing data. Taken together, this study provides a first insight into the

AC C

EP

TE D

M AN U

SC

RI PT

23

2

ACCEPTED MANUSCRIPT development of the testis with dietary energy restriction in sheep and shows that these

46

changes are associated with alterations in transcriptomic. The sheep testis mRNA database

47

were extended in this study will provides novel candidate regulators for future genetic and

48

molecular studies on sheep testicular development associated with energy restriction, which

49

will contribute to improving the reproductive performance of sheep.

RI PT

45

50

Key Words: Transcriptome; Testis development; Energy restriction; Compensatory feeding;

52

Hu sheep

M AN U

SC

51

53 54 55

59 60 61 62 63

EP

58

AC C

57

TE D

56

64 65 66 3

ACCEPTED MANUSCRIPT 67

1. Introduction Nutrition is an important environmental factor that affects reproductive performance and

69

breeding efficiency in sheep [1]. Reproductive performance may be severely impaired if the

70

diet is insufficient, unbalanced, or otherwise deleterious [2]. Avoiding nutritional scarcity

71

during the critical stage of mammalian reproduction presents a challenge because the supply

72

of forage varies considerably throughout the year. Thus, the correct feeding of rams is critical

73

for population growth. However, in this regard, research on the nutritional requirements for

74

ewe reproduction vastly outweighs that on rams [2].

M AN U

SC

RI PT

68

Previous studies have indicated that low nutrition can reduce the testicular weight and

76

production of spermatozoa in 1-year-old rams [3]. Larger amounts of testosterone produced

77

by the testis were found to reach the peripheral circulation in rams fed with a high diet (1 kg

78

of the basal diet plus 1.5 kg lupins /day) than in those fed with the low diet (400 g/day of the

79

basal diet) [3]. Previous studies have also shown that 72 % of the maintenance diet reduced

80

the single ejaculate volume and sperm motility rate of 2-year-old sheep [2]. However, contrast

81

with the above results, other studies have shown that nutritional treatment has no effect on the

82

levels of testosterone and testicular development [4], and Martin et al. [5] have shown that,

83

with the exception of severe underfeeding, there are only minor changes in the endocrine

84

function of the testis (testosterone production) unless season-long treatments are imposed.

85

Despite these findings, our understanding of the effects of nutrition on reproductive

86

performance remains limited. Reproduction is a complex process, and traits such as testicular

87

weight and reproductive hormone levels are affected by some genes, and thus it is very

88

important to gain an understanding of the gene regulatory mechanisms underlying the effect

AC C

EP

TE D

75

4

ACCEPTED MANUSCRIPT 89

of nutrition on sheep reproductive development. Hu sheep is an endemic species bred widely in the south of China. Due to their high

91

fecundity, Hu sheep have been introduced in northern China, where forage shortage is

92

seriously deficient during the cold winters. Thus, compensatory growth during the second

93

spring, facilitated by grazing on rich pastures, is essential to ensure the fattening/growth of

94

lambs. Although previous studies have generally focused on studying the effects of energy

95

restriction and compensation on growth performance[6], the mechanisms whereby energy

96

restriction and compensation affect sheep testicular reproductive performance remain

97

unknown. High-throughput sequencing (RNA-seq) technologies provide a unique opportunity

98

to characterize cell transcripts, to quantify transcripts and to identify di erential regulation in

99

a single experiment[7]. In recent years, RNA-seq technology has been applied to the study of

100

testicular development [8, 9], and is likely to provide valuable insights into the effects of

101

nutritional deficiency on testicular development in rams.

TE D

M AN U

SC

RI PT

90

In order to analyze the effects of different level of dietary energy restriction and

103

compensation on the testicular development of Hu sheep, we measured testicular weight and

104

levels of follicle stimulating hormone (FSH), luteinizing hormone (LH), and testosterone after

105

feeding the sheep on energy-restriction and compensatory diet. We also investigated the

106

impact of dietary energy restriction on transcriptome profile of Hu sheep testis, and analyzed

107

the molecular regulation of testicular development at different energy levels. Our results

108

provide useful information on the expression of reproductive-related genes prior to sexual

109

maturity.

110

2. Materials and methods

AC C

EP

102

5

ACCEPTED MANUSCRIPT 111

2.1 Animals and Experimental Design All protocols involving the use of animals were in accordance with the Guidelines for

113

Animal Experiments of Nanjing Agricultural University, and approved by the Animal Care

114

and Use Committee of Nanjing Agricultural University.

RI PT

112

The 3-month-old male Hu sheep (n=80, 22.02 ± 0.14 kg) used in this study were

116

randomly assigned to four groups: a normal energy group [Control ， 100% energy

117

requirements based on feeding standards for meat-producing sheep and goats (NY/T

118

816-2004, China); n = 20], a 85% maintenance energy group (ER1; n = 20), a 70%

119

maintenance energy group (ER2; n = 20), and a 55% maintenance energy group (ER3; n = 20).

120

The diet energy levels during the initial period of the study were 11.64 (Control), 9.73 (ER1),

121

8.40 (ER2), and 6.84 (ER3) MJ/kg. After 2 months of feeding, half the lambs in each group

122

(10 lambs per group) were euthanized, whereas during the following 3-month compensatory

123

period, the remaining lambs in all groups were fed the same diet containing 12.44 MJ/kg

124

energy, with water provided to the sheep ad libitum. The compensatory groups were

125

designated as follows: the control group (Control-C, 100% energy requirements prior to the

126

compensation period, n = 10), the 85% energy group with energy compensation (ER1-C; n =

127

10), the 70% energy group with energy compensation (ER2-C; n = 10), and the 55% energy

128

group with energy compensation (ER3-C; n = 10). The amount of dietary energy provided to

129

the compensated animals corresponded to the feeding standards for meat-producing sheep and

130

goats (NY/T 816-2004, China). The compositions of the experimental diets are shown Table

131

1.

132

2.2 Measurement of Testes Weight and Blood levels of Reproductive Hormones

AC C

EP

TE D

M AN U

SC

115

6

ACCEPTED MANUSCRIPT After the dietary restriction and compensatory periods, all the lambs were weighed

134

before slaughter. Blood samples were collected from the jugular vein by venipuncture at 0500,

135

0700, and 1100 h at the end of the dietary restriction and compensatory periods. The collected

136

blood was immediately placed in heparinized tubes, maintained on ice, separated by

137

centrifugation for 10 min at 3000 r/min at 4°C, and stored at –20°C for subsequent

138

reproductive hormone (FSH, LH and testosterone) measurements. For each lamb, a “testes

139

ratio” (i.e., testes to body weight ratio) was calculated. On the basis of significant differences

140

between the control and ER3 group, three testicle samples were selected randomly from each

141

group for RNA-Sequence analysis.

SC

M AN U

142

RI PT

133

2.3 Collection of Testicular Tissue and RNA Extraction

All the testicular tissues were collected from the same area, placed in Sample Protector

144

for RNA/DNA (Takara, Dalian, China), snap-frozen in liquid nitrogen, and stored in -80°C

145

cryogenic refrigerator for preservation. Total RNA was extracted from the testicular tissues of

146

both control and ER3 group sheep using a high purity total RNA extraction kit (Bioteke,

147

Beijing, China). The purified RNA was monitored on 1% agarose gels to ensure no genomic

148

DNA contamination. RNA purity was checked using a NanoPhotometer® spectrophotometer

149

(IMPLEN, CA, USA), and RNA concentration was measured using a Qubit® RNA Assay Kit

150

in Qubit® 2.0 Fluorometer (Life Technologies, CA, USA). RNA integrity was assessed using

151

an RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 System (Agilent Technologies,

152

Santa Clara, CA, USA).

153

2.4 Construction of Library and Illumina Sequencing

154

AC C

EP

TE D

143

Equal amounts of RNA samples (3 µg RNA per sample) were used to prepare a cDNA 7

ACCEPTED MANUSCRIPT library. In brief, mRNA was purified from total RNA using poly-T oligo-attached magnetic

156

beads and fragmented. First-strand cDNA was generated using random hexamer primers and

157

M-MuLV Reverse Transcriptase (RNase H-). Second-strand cDNA synthesis was

158

subsequently performed using DNA Polymerase I and RNase H. Following the addition of a

159

single A base and adapters, the blunt-end cDNA fragments were subsequently ligated.

160

Thereafter, the fragments were amplified by PCR to create the final cDNA library.

RI PT

155

After library construction, Qubit2.0 was used to perform a preliminary analysis and the

162

library concentration was diluted to 1 ng/µL. An Agilent Bioanalyzer 2100 was used to

163

examine the length of inserted fragments in the library and after reaching the expected insert

164

size, accurate quantification of the effective concentration of the library (library effective

165

concentration > 2 nM) was carried out using qRT-PCR. The clustering of index-coded

166

samples was performed using a cBot Cluster Generation System in conjunction with a TruSeq

167

PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After

168

cluster generation, the library preparations were sequenced using the Illumina Hiseq X Ten

169

platform and 125 bp/150 bp paired-end reads were generated.

170

2.5 Quality Evaluation of Sequencing Data

AC C

EP

TE D

M AN U

SC

161

171

Base sequencing error rates were determined from Phred numerical score (Phred score,

172

Qphred) and by formula [formula1: Qphred = -10 log10 (e)] transformation. The GC content

173

was then determined to ensure that AT and GC were separated. After removing low-quality

174

reads from the raw data, clean reads were obtained by removing reads containing adapters and

175

ploy-N. These clean reads were used for subsequent analyses.

176

In this study, we used TopHat2 software to select a gene model for genome positioning 8

ACCEPTED MANUSCRIPT analysis of the filtered sequences. Reads from different regions of the reference genome

178

(exons, introns, and intergenic regions) were analyzed. The read positions on chromosomes,

179

length of chromosomes, the relation of the length of chromosome and the total number of

180

reads were visualized. Pearson correlation analysis was used to perform RNA - Sequence

181

correlation tests and R language was used to calculate the Pearson correlation coefficients.

182

2.6 Alternative Splicing Analysis and Prediction of Novel Transcripts Alternative

splicing

events

were

analyzed

using

rMATS

SC

183

RI PT

177

analysis

software

(http://rnaseq-mats.soueceforge.net/index.html). Different types of alternative splicing events

185

were analyzed using the quantitative methods of junction count only (JC only) and reads on

186

target and junction counts (JC + reads on target), where an FDR < 0.05 was used as the

187

screening criterion. Reads were then mapped against ovine genome assembly v.3.1 (Oar_v3.1).

188

The Cufflinks v2.1.1 Reference Annotation Based Transcript (RABT) assembly method was

189

used to construct and identify both known and novel transcripts from the TopHat alignment

190

results.

191

2.7 Single Nucleotide Polymorphisms (SNP) and Insertion and Deletion (InDel) Analyses

EP

TE D

M AN U

184

Picard-tools v1.96 and samtools v0.1.18 were used to sort, mark duplicated reads and

193

reorder the bam alignment results for each sample. GATK2 (v3.2) software was used to

194

perform SNP and InDel calling. Finally, analysis results were obtained after filtering.

195

2.8 Quantification of Gene Expression Level

AC C

192

196

Gene expression levels were estimated in terms of fragments per kilo base of transcript

197

sequence per millions of base pairs sequenced (FPKM). HTSeq v0.6.1 software and a union

198

model were used to analyze the expression levels of genes in each sample. 9

ACCEPTED MANUSCRIPT Gene Ontology（GO）and Kyoto Encyclopedia of Genes and Genomes（KEGG）

199

2.9

200

Enrichment Analysis of Differentially Expressed Genes To identify the biological pathways and functions of differentially expressed genes

202

(DEGs), GO and KEGG analyses were performed on DEGs. We obtained GO annotation

203

terms and performed GO functional enrichment using GOseq software based on the Wallenius

204

non-central hyper-geometric distribution [10]. KEGG pathway analysis was performed using

205

KOBAS software (http://kobas.cbi.pku.edu.cn/home.do) against the KEGG database.

206

2.10 Validation of DEGs

M AN U

SC

RI PT

201

The total RNA was extracted from testis samples using a Baiteke RNA extraction reagent

208

kit (Baiteke, Beijing, China) and cDNA synthesis was performed using an RNA reverse

209

transcription reagent kit (Takara, Shanghai, China). Six genes related to reproduction were

210

selected to confirm the sequencing data by qRT-PCR. The sequences and NCBI Reference

211

Sequence of the primers used for amplification of the target genes are shown in Table S1. The

212

expression level of the six DEGs were compared between the following pairs of groups:

213

Control vs ER3, Control-C vs ER3-C, Control vs Control-C, and ER3 vs ER3-C. The

214

qRT-PCR was performed using the Step One Plus qRT-PCR System (Applied BioSystems

215

Carlsbad, USA) and fluorescence was detected using SYBR Green (No.04913850001; Roche,

216

Mannheim, Germany) in a reaction volume of 20 µL. The comparative quantification of each

217

RNA was normalized to glyceraldehyde 3-phosphate dehydrogenase (GAPDH) using the

218

2-∆∆Ct method and all reactions were carried out three times.

219

2.11 Statistical Analyses

220

AC C

EP

TE D

207

All data were analyzed using SPSS 19.0 (SPSS Inc. Chicago, IL, USA) and are presented 10

ACCEPTED MANUSCRIPT as mean values ± SEM from 10 different individuals (n=10) with three biological replicates.

222

The testicular weight, slaughter weight, testes ratio, and blood levels of reproductive

223

hormones were compared using single factor analysis of variance, and qRT-PCR results were

224

analyzed using Student’s t-test. Values of P < 0.05 were considered to indicate a statistically

225

significant difference and those of P < 0.01 were considered to indicate a statistically

226

extremely significant difference.

227

3. Results

228

3.1 Effects of Restricted and Compensatory Feeding on Testicular Development and

229

Reproductive Hormones in the Blood

M AN U

SC

RI PT

221

After a period of restricted feeding, the body and testicular weights of lambs decreased

231

with an increase in energy restriction (P < 0.05). The testes ratio in the control group was

232

significantly higher than that in the ER2 and ER3 groups (P < 0.05; Table 2). Although there

233

were no significant differences in the LH and FSH levels among the different groups (P >

234

0.05), testosterone levels in the control, ER1, and ER2 groups were significantly higher than

235

those in the ER3 group (P < 0.05) (Table 3). After compensatory feeding, there was no

236

significant difference in the above indices of testicular development and reproductive

237

hormones among the groups (P > 0.05) (Tables 4-5).

238

3.2 Analysis of RNA-sequence Data

AC C

EP

TE D

230

239

In order to obtain an overview of the testis transcriptome of sheep and identify the genes

240

involved in testis development under dietary energy restriction, six cDNA libraries were

241

prepared from RNA extracted from testes samples collected from control and ER3 group rams

242

and sequenced using the Illumina Hiseq platform. In order to ensure that the RNA sequencing 11

ACCEPTED MANUSCRIPT data satisfied the criteria for genome-wide transcriptomic analysis, we conducted standard

244

analyses of a quality control. A total 181,786,368 (26.36 G) and 157,273,776 (22.81 G) raw

245

reads were obtained for the control group and ER3 group, respectively (Table S2). Of these,

246

60.9%–69.7% of the total reads were mapped to exons in genomic regions. After removing

247

low-quality reads with a Qpred quality of less than 20, a total of 175,771,456 and 152,109,294

248

clean reads were filtered to guarantee the sequencing data. The average error rate of base

249

sequencing was approximately 0.01 %, indicating a high sequencing quality. Total clean reads

250

were mapped to the reference genome, with approximately 77.63%–79.26% of the clean reads

251

being aligned successfully (Table S3). Under normal circumstances, the greater the length of

252

the chromosome, the greater the number of reads located within the chromosome [11]. In the

253

present study, we found that testicular development-related reads were enriched in

254

chromosomes 1, 2, 3 and X (Figures S1a-S1f). The correlation of gene expression in different

255

samples is an important index for assessing the reliability of the experiment and the selection

256

of samples. Here, the correlation coefficient within a group was close to 1, indicating that the

257

expression pattern between samples was highly similar. All the squares of the Pearson

258

correlation coefficients between samples were greater than 0.8 (Figure S2), meeting the

259

standard requirements.

AC C

EP

TE D

M AN U

SC

RI PT

243

260

FPKM provides a measure of gene expression levels. On the basis of the determined

261

FPKM values, gene expression levels were classified into five categories (Table S4). Among

262

these, the largest proportion (33.84%–36.86%) of genes exhibited low expression (0 < FPKM

263

< 1), followed by the second largest group (21.71%–23.99%) with FPKM values greater than

264

3 and less than 15. Only a small proportion of genes were highly expressed. 12

ACCEPTED MANUSCRIPT 3.3 Alternative Splicing, SNP, and InDel Analyses

266

The degree of alternative splicing indicates the enrichment of protein complexity. In this

267

study, five primary types of alternative splicing event were detected (Figure 1): skipped exon,

268

mutually exclusive exon, alternative 5′ splice site, alternative 3′ splice site, and retained intron.

269

Junction count only (JC only) and reads on target and junction counts (JC + read on target)

270

were used to detect the type of alternative splicing, and they yielded similar results: the vast

271

majority of alternative splicing events (<10,000) were skipped exon events, with mutually

272

exclusive exon events being the second most abundant events. Alternative 5′splice site,

273

alternative 3′ splice site, and retained intron events made up less than 100 of the alternative

274

splicing events. Among the 430,000–480,000 SNPs, 37,000–40,000 InDels were from one

275

sample in this transcriptome dataset (Table S5).

276

3.4 DEGs between Control and ER3 Groups

TE D

M AN U

SC

RI PT

265

The FPKM distribution of all genes in the ER3 and Control groups were similar (Figure

278

2a). In the six examined samples, 29,307 genes were detected, and 261 DEGs with fold

279

changes greater than two were obtained (P < 0.05, FPKM > 1). Among the DEGs, 81 genes

280

were upregulated, whereas the remaining 180 genes were downregulated (Figure 2b) in the

281

ER3 group compared with the control group. Genes with similar expression patterns or

282

enriched in the same metabolic process or cellular pathway may have similar function, and

283

thus the function of unknown genes or new functions of known genes can be predicted by

284

clustering the genes into classes. In this study, a hierarchical cluster analysis was performed to

285

classify the DEGs (Figure 2c).

286

3.5 GO and KEGG Analysis of Differentially Expressed Genes

AC C

EP

277

13

ACCEPTED MANUSCRIPT The GO and KEGG databases were used for gene annotation (Tables S3 and S4). The top

288

30 significantly enriched GO terms (ER3 vs Control groups) are shown in Figure 3a. A total

289

of 18,134 (61.88%) genes were assigned to 162 cellular components, 251 molecular functions

290

and 1392 biological processes. Molecular function (GO: 0003674), biological process (GO:

291

0008150), and binding (GO: 0005488) were the three terms enriched with the most genes.

RI PT

287

KEGG annotation is the process of mapping interesting genes to metabolic pathways. In

293

our study, 261 DEGs were mapped to 92 KEGG pathways. The top 20 significantly enriched

294

KEGG pathways are shown in Figure 3b, and among these pathways, we identified two

295

signaling pathways that have been well-documented to be essential in testes development,

296

namely, the ovarian steroidogenesis pathway, and the drug metabolism-cytochrome P450

297

pathway. In addition, a further 11 important pathways related to the regulation of male sexual

298

function were identified, including the mitogen-activated protein kinase pathway,

299

phosphatidylinositol-3-kinase protein kinase B (PI3K-Akt) pathway, AMPK pathway and Wnt

300

pathway, and so on (Table 6). The significantly enriched GO and KEGG results suggested that

301

some DEGs were associated with testes development in Hu sheep. Thus, it can be speculated

302

that energy restriction may affect male reproduction by altering the transcriptional profile of

303

the testes.

304

3.6 Verification of DEGs by qRT-PCR

AC C

EP

TE D

M AN U

SC

292

305

To validate the RNA-Sequence results, the expression levels of six DEGs primarily

306

involved in ram reproduction were analyzed using qRT-PCR for the following group pairs:

307

Control vs. ER3, Control-C vs. ER3-C, Control vs. Control-C, and ER3 vs. ER3-C (Figure 4).

308

The expression level of the six genes showed a pattern that correlated well with that of the 14

ACCEPTED MANUSCRIPT sequencing data (Figure 4a and 4b). The expression level of INSR was apparently upregulated

310

in sheep fed the 55% energy (Figure 4a) and compensation diets (Figure 4c), and also

311

upregulated in the Control-C group compared with the Control group, and the ER3-C group

312

compared with the ER3 group. The CYP11A1, FDFT1, and Wnt2 genes showed similar

313

expression tendencies, wherein all were downregulated after feed restriction and

314

compensation, but were upregulated following compensatory feeding in the Control-C and

315

ER3-C groups compared with the Control and ER3 groups, respectively (Figure 4d and 4e).

316

Although expression of the ALDH3B1 gene recovered to normal levels after low expression

317

during restricted feeding, it was downregulated in the Control-C and ER3-C groups compared

318

with the Control and ER3 groups, respectively, thereby indicating that it might play a role in

319

the development of the testes prior to sexual maturity. The compensatory growth phase

320

significantly increased the low expression levels measured during restricted feeding. Under

321

restricted feeding and normal feeding conditions, expression levels of the PGR gene were

322

downregulated in the ER3 and Control-C groups compared with the Control group,

323

respectively, whereas after compensatory feeding, PGR expression was significantly

324

upregulated.

325

4

AC C

EP

TE D

M AN U

SC

RI PT

309

Discussion

326

In this study, we measured testicular weight and blood levels of reproductive hormones

327

in Hu rams after restricted and compensatory feeding. A calorie-restricted diet affects

328

reproduction by altering the metabolism in young mammals [5]. Previous studies have

329

indicated that energy restriction decreases the secretion of pituitary peptide hormones, thereby

330

reducing the transport of cholesterol from intracellular stores to the mitochondria via the 15

ACCEPTED MANUSCRIPT cyclic adenosine monophosphate signal pathway. As result, testosterone etc synthesis

332

decreased. Reduced synthesis of the steroid hormone testosterone can influence sperm count

333

and motility by acting on Sertoli cells or spermatogonia. The results of the present study,

334

showing that testicular weight and ratio decreased with increasing dietary restriction are

335

consistent with this decreased secretion of pituitary hormones. Unexpectedly, however, we

336

observed that blood levels of LH and FSH also had a slight tendency to decline. After

337

compensatory feeding, the levels of testosterone returned to normal, as did testicular weight

338

and testes ratio. These findings indicate that compensatory feeding can compensate for the

339

decline in reproductive performance due to energy restriction during the early stages of

340

growth. These findings suggest that energy-compensated feeding can restore the pre-puberty

341

reproductive performance of lambs to normal levels, and that early restricted feeding would

342

not adversely affect the later breeding performance of adult rams.

TE D

M AN U

SC

RI PT

331

In recent years, high-throughput sequencing has been used in bioinformatics analysis of

344

humans, mice, maize, cattle, and other species [12, 13]; however, analysis of the sheep testis

345

transcriptome has rarely been carried out. To further study the effects of feed restriction on

346

reproduction, RNA-sequencing was used to examine the transcriptional differences between

347

the testes of lambs fed energy-limited or maintenance energy diets. From six samples we

348

identified 261 genes that were differentially expressed between the ER3 group and the control

349

group, with 81 being upregulated and 180 downregulated in the former compared with the

350

latter. Some of these DEGs, including CYP11A1, INSR, ALDH3B1, are related to reproduction,

351

and our finding indicate that these DEGs might play important roles in the regulation of

352

reproductive hormone secretion, testicular development and spermatogenesis.

AC C

EP

343

16

ACCEPTED MANUSCRIPT 353

KEGG analysis indicated that, the aforementioned DEGs are enriched in multiple

354

reproduction-related pathways. Ovarian steroidogenesis and steroid hormone biosynthesis

355

pathways play important roles in the regulation of mammalian reproduction, and can regulate

356

the expression of related proteins and promote the synthesis of testosterone, which

357

can promote testes development [14]. The oocyte meiosis pathway and adherens junction

358

pathway may be related to meiosis and mitosis, which affect testicular size. Furthermore, the

359

Wnt signaling pathway can promote self-renewal of spermatogonial stem cells,

360

spermatogenesis and proliferation of progenitor cell populations in mammalian testes [15].

361

Disrupting this pathway in humans can lead to infertility and testicular cancer [15]. The

362

mitogen-activated protein kinase pathway is related to the development of germ line stem

363

cells [16]. Previous studies have shown that a downstream transcriptional factor of p38

364

mitogen-activated protein kinase is a key mediator of extracellular stimulus-induced testicular

365

injury [17]. There are also many functional genes involved in reproductive hormone

366

production in the steroid hormone biosynthetic pathway and the drug metabolism-cytochrome

367

P450 pathway. The steroidogenic hormone biosynthesis pathway can influence the secretion

368

of testosterone in Leydig cells by affecting the synthesis of testosterone precursors [18], and

369

cytochrome P450 enzymes can metabolize multiple sterols in vitro to establish novel

370

branching points in cholesterol synthesis for conversion to testosterone [19]. These findings

371

emphasize the importance of reproductive hormone production and testicular size as well as

372

their corresponding genes in the development of testes.

AC C

EP

TE D

M AN U

SC

RI PT

in turn

373

Alternative splicing is a source of protein diversity and biological complexity [20] , and

374

among the six samples collected in the present study, we detected five primary types of 17

ACCEPTED MANUSCRIPT alternative splicing event and a number of genes containing one or more alternative splices.

376

Among these, the skipped exon and mutually exclusive exon events comprised the largest

377

proportion. Alternative splicing is also considered as a standard for detecting the complexity

378

of cell and function in eukaryotes [21]. Most of the genomic variation is attributed to SNPs

379

and InDels, which may also increase functional complexity and offer the highest resolution

380

for tracking disease-related genes and population history [22]. In the present study, we

381

detected 430,000–480,000 SNPs, and 37,000–40,000 InDels among the different samples,

382

which are useful for the identification of candidate genes and biomarkers for sheep fertility

383

[23], and should be useful as molecular markers for functional genomics and breeding

384

research in this species and other closely related species [24]. Therefore, the large number of

385

SNPs and InDels that we detected in the testes warrant further analysis

M AN U

SC

RI PT

375

The qRT-PCR was performed to validate the expression of reproduction-associated genes

387

selected from among the DEGs. The results were consistent with RNA sequencing data. In

388

this study, CYP11A1 was enriched in the steroidogenesis and steroid hormone biosynthesis

389

pathways. P450sec the cholesterol side chain cleavage enzyme encoded by the CYP11A1 gene

390

is instrumental in the synthesis of sex hormones, and is known to catalyze the initial step of

391

steroidogenesis

392

monophosphate–mediated signals [25]. Additionally, desmosterol can be converted to

393

pregnenolone by CYP11A1 in another branch of the cholesterol pathway leading toward

394

steroid hormones [26]. Therefore, we hypothesized that the low expression of CYP11A1 in

395

pubertal Hu sheep after energy-restricted feeding may have resulted in the low levels of

396

testosterone recorded in the present study. The enzyme encoded by FDFT1 is involved in the

AC C

EP

TE D

386

in

the

adrenals

and

18

gonads

through

cyclic

adenosine

ACCEPTED MANUSCRIPT production of C-30 squalene in a reaction that is also considered to be a branch of the

398

isoprenoid pathway. Previous studies have shown that a mutation in the FDFT1 gene of

399

heterozygous individuals can lead to a 50% reduction in enzyme activity in the testes. A

400

further gene, PGR, plays a vital role in the biosynthetic pathway of sex hormones [27]. Thus,

401

although the genes CYP11A1, FDFT1 and PGR are capable of regulating the synthesis of

402

steroid hormones , they were all expressed at low levels under the condition of energy

403

restriction imposed in the present study, indicating that they may cooperatively cause low

404

testosterone levels.

M AN U

SC

RI PT

397

Wnt family members play roles in the development of sexual dimorphism in the

406

mammalian gonads[28]. Among these, the Wnt2 gene plays an important role in the formation

407

of precursors in the early stage of gonad development, particularly during the embryonic

408

period [29, 30], and this may contribute to the development of spermatogonial stem cells,

409

Sertoli cells or mesenchymal cells. Thus, a decrease in Wnt2 gene expression in the ER3 rams

410

may have contributed to the reduced testicular weight of pubertal Hu sheep.

TE D

405

ALDH3B1 plays an important defensive role against oxidative and aldehyde stressors

412

that may impair spermatogenesis. Therefore, expression of the ALDH3B1 gene is indicative of

413

sperm production capacity [31]. Although, lamb semen was not collected in the present study,

414

low expression of the ALDH3B1 gene might indicate a decrease in testicular sperm

415

production. INSR gene expression affects the insulin/IGF-1 signaling pathway, which is

416

required for FSH-mediated spermatogenic cells proliferation, and may play an important role

417

in regulating the final number of spermatogenic cells, testicular size, and daily sperm output

418

[32]. In the present study, FSH levels in the blood were unchanged, whereas the expression of

AC C

EP

411

19

ACCEPTED MANUSCRIPT INSR was upregulated. Although we expected that this would promote an increase in

420

testicular size, no such increase was observed. Potentially, this finding may be attributed to

421

severe energy restriction having a negative effect on the molecular mechanism of the

422

signaling pathway; however, the specific reasons remain unknown and require further study.

423

After compensatory feeding, expression of the ALDH3B1 gene returned to normal levels,

424

whereas that of the PGR gene was significantly increased. In contrast, the other four DEGs

425

examined showed expression patterns similar to those observed following restricted feeding.

426

This indicates that 55% energy-restricted feeding for2 months had a significant effect on the

427

expression of the INSR, CYP11A1, FDFT1, PGR and Wnt2 genes, and that the expression of

428

these genes did not return to normal levels after 3 months of compensatory feeding. At

429

different stages of animal development, gene expression levels may vary widely. Male Hu

430

sheep reach puberty at 120 days[33], which is earlier than that of other breeds[34]. Therefore,

431

in the present study, the lambs were sexually mature after the restricted period and

432

compensatory feeding. We found that expression of the CYP11A1, FDFT1, ALDH3B1 and

433

Wnt2 genes in the control and ER3 groups showed the same patterns at different

434

developmental stages at 5 months and 8 months of age. This suggests that previous energy

435

restriction did not change the expression of these four genes prior to maturation. The

436

expression of PGR was significantly increased at 8 -months of age after restricted and

437

compensatory feeding, which was independent of an increase in age, as PGR expression

438

decreased with age in the control group. Therefore, we speculate that the PGR gene may

439

contribute to the development of testes during the phase of growth prior to sexual maturation.

440

However, the underlying molecular mechanisms require further study.

AC C

EP

TE D

M AN U

SC

RI PT

419

20

ACCEPTED MANUSCRIPT 441

5

Conclusions Overall, this study represents the first utilization of Illumina sequencing technology to

443

conduct a transcriptome analysis of sheep testis associated with dietary energy restriction.

444

RNA-seq enabled us to obtain 261 DEGs, among which many genes potentially involved in

445

the regulation of reproduction were identified. The expression patterns of these genes could

446

serve to explain the di erences in reproductive characteristics (testicular weight and

447

testosterone concentration) caused by energy restriction and compensation. The finding of this

448

study shows that energy restriction induced significant transcriptomic alterations in Hu sheep

449

testis, and provides a valuable transcriptomic resource for Hu sheep and a basis for future

450

research on the molecular mechanisms underlying the relationships between nutrition and

451

reproduction.

452

Supplementary Materials

TE D

M AN U

SC

RI PT

442

Figure S1: The number of reads mapped to the chromosomes of Hu sheep fed on 55%

454

and 100% (control) energy requirement diets. Figure S2: Pearson correlation analysis between

455

the testes samples of Hu sheep fed on 55% and100% (control) energy requirement diets. Table

456

S1: Details of primer sequences, expected product sizes and GenBank accession numbers of

457

genes used for qRT-PCR. Table S2: Sequencing data quality assessment for Hu sheep fed on

458

55 % and 100% (control) energy requirement diets. Table S3: Reads compared with the

459

reference genome of Hu sheep fed on 55% and 100% (control) energy requirement diets.

460

Table S4: Number of genes with different expression levels in Hu sheep fed on 55% and

461

100% (control) energy requirement diets. The RNA-Seq data from this study have been

462

submitted to the GEO database with the accession number is GSE111453.

AC C

EP

453

21

ACCEPTED MANUSCRIPT 463

Acknowledgements We are grateful to all members of the F. Wang laboratory for their contributions to

465

sample determination. The research was supported by the Special Fund for Agro-scientific

466

Research in the Public Interest (No.201303144) and National Industrial Technology System

467

of Sheep & Goat (CARS-38).

468

References

469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501

[1] Luo F, Jia R, Ying S, Wang Z, Wang F. Analysis of genes that influence sheep follicular development by

RI PT

464

SC

different nutrition levels during the luteal phase using expression profiling. Animal Genetics. 2016;47:354-64. [2] Guan YJ, Malecki IA, Hawken PAR, Linden MD, Martin GB. Under-nutrition reduces spermatogenic efficiency and sperm velocity, and increases sperm DNA damage in sexually mature male sheep. Animal reproduction

M AN U

science. 2014;149:163-72.

[3] Hotzel MJ, Markey CM, Walkden-Brown SW, Blackberry MA, Martin GB. Morphometric and endocrine analyses of the effects of nutrition on the testis of mature Merino rams. J Reprod Fertil. 1998;113:217-30. [4] Martin GB, Tjondronegoro S, Blackberry MA. Effects of nutrition on testicular size and the concentrations of gonadotrophins, testosterone and inhibin in plasma of mature male sheep. J Reprod Fertil. 1994;101:121-8. [5] Martin GB, Blache D, Miller DW, Vercoe PE. Interactions between nutrition and reproduction in the management of the mature male ruminant. animal. 2010;4:1214-26.

[6] Fontes CAD, Guimaraes RFM, de Almeida MLV, de Campos OF, de Almeida FQ, Sant'Ana ND. Evaluation of 2007;36:698-708.

TE D

compensatory growth in crossbred Holstein-Gyr steers: intake and performance. Rev Bras Zootecn. [7] Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10:57-63.

[8] Meng XL, Liu P, Jia FL, Li J, Gao BQ. De novo Transcriptome Analysis of Portunus trituberculatus Ovary and

EP

Testis by RNA-Seq: Identification of Genes Involved in Gonadal Development. PLoS One. 2015;10:e0128659. [9] Wu J, Xiong S, Jing J, Chen X, Wang W, Gui JF, et al. Comparative Transcriptome Analysis of Differentially Expressed Genes and Signaling Pathways between XY and YY Testis in Yellow Catfish. PLoS One.

AC C

2015;10:e0134626.

[10] Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology. 2010;11:: R14. [11] Marquez Y, Brown JW, Simpson C, Barta A, Kalyna M. Transcriptome survey reveals increased complexity of the alternative splicing landscape in Arabidopsis. Genome Research. 2012;22:1184-95. [12] Chao T, Wang G, Wang J, Liu Z, Ji Z, Hou L, et al. Identification and Classification of New Transcripts in Dorper and Small-Tailed Han Sheep Skeletal Muscle Transcriptomes. Plos One. 2016;11:e0159638. [13] Venugopal S, Martinez-Arguelles DB, Chebbi S, Hullin-Matsuda F, Kobayashi T, Papadopoulos V. Plasma membrane origin of the steroidogenic pool of cholesterol used in hormone-induced acute steroid formation in Leydig cells. Journal of Biological Chemistry. 2016;291. [14] Chand D, Colacci M, Dixon K, Kollara A, Brown TJ, Lovejoy DA. C-terminal region of teneurin-1 co-localizes with the dystroglycan complex in adult mouse testes and regulates testicular size and testosterone production. Histochemistry and Cell Biology. 2014;141:191-211. 22

ACCEPTED MANUSCRIPT [15] Kumar M, Atkins J, Cairns M, Ali A, Tanwar PS. Germ cell-specific sustained activation of Wnt signalling perturbs spermatogenesis in aged mice, possibly through non-coding RNAs. Oncotarget. 2016. [16] Amoyel M, Anderson J, Suisse A, Glasner J, Bach EA. Socs36E Controls Niche Competition by Repressing MAPK Signaling in theDrosophilaTestis. Plos Genetics. 2016;12:e1005815. [17] Banerjee B, Nandi P, Chakraborty S, Raha S, Sen PC, Jana K. Resveratrol ameliorates benzo(a)pyrene-induced testicular dysfunction and apoptosis: involvement of p38 MAPK/ATF2/iNOS signaling. Journal of Nutritional Biochemistry. 2016;34:17-29.

RI PT

[18] Mei H, Guan X, Wei L, Li P, Yang M, Liu C. Di-(2-ethylhexyl) phthalate inhibits testosterone level through disturbed hypothalamic–pituitary–testis axis and ERK-mediated 5α-Reductase 2. Science of the Total Environment. 2016;563-564:566-75.

[19] Motohashi M, Wempe MF, Mutou T, Okayama Y, Kansaku N, Takahashi H, et al. In utero-exposed di(n-butyl) phthalate induce dose dependent, age-related changes of morphology and testosterone-biosynthesis enzymes/associated proteins of Leydig cell mitochondria in rats. Journal of Toxicological Sciences. 2016;41. read sequencing rate. BMC Bioinformatics. 2015;16:1-13.

SC

[20] Liu X, Shi X, Chen C, Li Z. Improving RNA-Seq expression estimation by modeling isoform- and exon-specific [21] Pan Q, Shai O, Blencowe BJ, Lee LJ, Frey BJ. Deep surveying of alternative splicing complexity in the human

M AN U

transcriptome by high-throughput sequencing. Nature Genetics. 2009;40:1413-5.

[22] Altshuler D, Pollara VJ, Cowles CR, Van Etten WJ, Baldwin J, Linton L, et al. An SNP map of the human genome generated by reduced representation shotgun sequencing. Nature. 2000;407:513-6. [23] Fischer D, Laiho A, Gyenesei A, Sironen A. Identification of Reproduction-Related Gene Polymorphisms Using Whole Transcriptome Sequencing in the Large White Pig Population. G3 (Bethesda, Md). 2015;5. [24] Ding H, Luo Y, Liu M, Huang J, Xu D. Histological and transcriptome analyses of testes from Duroc and Meishan boars. Scientific Reports. 2016;6:20758.

TE D

[25] Reddy KR, Deepika MLN, Supriya K, Latha KP, Rao SSL, Rani VU, et al. CYP11A1 microsatellite (tttta)n polymorphism in PCOS women from South India. Journal of Assisted Reproduction and Genetics. 2014;31:857-63.

[26] Jure A, Sandeep G, Rok K, Marko G, Martina P, Ales B, et al. Cytochrome P450 metabolism of the post-lanosterol intermediates explains enigmas of cholesterol synthesis. Scientific Reports. 2016;6:28462.

EP

[27] Fischer D, Laiho A, Gyenesei A, Sironen A. Identification of Reproduction-Related Gene Polymorphisms Using Whole Transcriptome Sequencing in the Large White Pig Population. G3-Genes Genom Genet. 2015;5:1351-60.

[28] Dong WL, Tan FQ, Yang WX. Wnt signaling in testis development: Unnecessary or essential? Gene.

AC C

502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545

2015;565:155-65.

[29] Defalco T, Camara N, Bras SL, Doren MV. Nonautonomous sex determination controls sexually dimorphic development of the Drosophila gonad. Developmental Cell. 2008;14:275-86. [30] Franke L, Werner R, Richter-Unruh A, Schwab KO, Hiort O, Holterhus PM. Potential roles of wingless-type-MMTV integration site family, member 2 (WNT2) in human male external genital development. Exp Clin Endocr Diab. 2007;115:S90-S. [31] Marchitti SA, Brocker C, Orlicky DJ, Vasiliou V. Molecular characterization, expression analysis, and role of ALDH3B1 in the cellular protection against oxidative stress. Free Radical Biology & Medicine. 2010;49:1432-43. [32] Pitetti JL, Calvel P, Zimmermann C, Conne B, Papaioannou MD, Aubry F, et al. An essential role for insulin and IGF1 receptors in regulating sertoli cell proliferation, testis size, and FSH action in mice. Molecular Endocrinology. 2013;27:814-27. [33] Yue GH. Reproductive characteristics of Chinese Hu sheep. Animal Reproduction Science. 1996;44:223-30. 23

ACCEPTED MANUSCRIPT [34] Dyrmundsson OR. Puberty and early reproductive performance in sheep rams lambs. Animal Breeding

EP

TE D

M AN U

SC

RI PT

Abstracts. 1973.

AC C

546 547 548 549 550 551 552 553 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

24

ACCEPTED MANUSCRIPT Figure legends： ：

591

Fig. 1. Detection of alternative splicing events in the groups of Hu sheep fed on diets

592

containing 55% (ER3) and 100% (control) of energy requirements (n = 3).

593

Alternative splicing events were detected using the junction counts and reads on target

594

methods. The numbers of five alternative splicing types showed similar percentages of the

595

total alternative splicing events using these two methods. SE, skipped exon; RI, retained

596

intron; MXE, mutually exclusive exon; A5SS, alternative 5′ splice site; A3SS, alternative 3′

597

splice site; JC only, junction count only; JC + read on target, reads on target and junction

598

counts.

599

Fig. 2. Gene expression and dynamic changes in differentially expressed genes (DEGs)

600

between Hu sheep fed on diets containing 55% (ER3) and 100% (control) of energy

601

requirements (n = 3).

602

(a) Genes expressed at different levels. Different-colored areas indicate gene expression of

603

different groups: red indicates the ER3 group and, green indicates the control group. (b)

604

Volcano plot of DEGs in the ER3 group compared with the control group. Different colors

605

indicate different levels of expression: green indicates upregulation and red indicates

606

downregulation. (c) Hierarchical clustering analysis of DEGs. Each column represents a

607

group sample (ER3 group or control group), and each row represents a gene. Changes in color

608

represent changes in gene expression, where red represents increased gene expression and

609

blue represents decreased gene expression. Darker colors indicate a greater difference in

610

expression.

611

Fig. 3. Gene ontology (GO) classification and Kyoto encyclopedia of genes and genomes

AC C

EP

TE D

M AN U

SC

RI PT

590

25

ACCEPTED MANUSCRIPT (KEGG) enrichment pathways of the groups of Hu sheep fed on diets containing 55%

613

(ER3) and 100% (control) of energy requirements (n = 3).

614

(a) DEGs were enriched in three categories, shown as different colors. Red indicates

615

biological process; blue indicates cellular component; green indicates molecular function. (b)

616

A KEGG pathway enrichment analysis between the DEGs of groups A and D. The different

617

colors represent different enrichment factors. The size of the plot corresponds to the number

618

of genes analyzed. DEGs, differentially expressed genes.

619

Fig.4. Validation of the expression levels of differentially expressed genes (DEGs) related

620

to the reproduction of Hu sheep by qRT-PCR (n = 10) and RNA-seq (n = 3).

M AN U

SC

RI PT

612

(a) The relative expression level of DEGs, INSR (Insulin receptor), CYP11A (cytochrome

622

P450 family 11 subfamily A member 1), ALDH3B1 (aldehyde dehydrogenase 3 family

623

member B1), FDFT1 (farnesyl-diphosphate farnesyltransferase 1), Wnt2 (Wnt family member

624

2), and PGR (progesterone receptor) in ER3 and Control groups determined by qRT-PCR. (b)

625

The relative expression level of DEGs in ER3 and Control groups determined by RNA-seq. (c)

626

The relative expression level of DEGs in Control-C and ER3-C groups determined by

627

qRT-PCR. (d) The relative expression level of DEGs in Control and Control-C groups

628

determined by qRT-PCR. (e) The relative expression level of DEGs in ER3 and ER3-C groups

629

determined by qRT-PCR. The relative expression levels of DEGs in testes were determined by

630

qRT-PCR and normalized to the expression of glyceraldehyde 3-phosphate dehydrogenase

631

(GAPDH). Control, 100% energy requirements; ER3, 55% energy group; Control-C, 100%

632

energy requirements during the compensation period; ER3-C, energy compensation after 55%

633

restricted feeding. * P < 0.05. ** P < 0.01.

AC C

EP

TE D

621

26

Revised

ACCEPTED MANUSCRIPT

Table. 1. Dietary composition and nutrient levels in Hu sheep（ （air-dry basis） ） Energy restriction period

Energy

Items compensatory period

Control

ER1

ER2

ER3

Barley husk

0

0

0

0

Rice straw

0

26.00

42.00

64.00

0

Soybean straw

30.00

16.00

10.50

0

0

Cassava residue

12.00

8.00

9.50

4.00

21.00

Corn

28.00

17.00

5.00

0

18.00

Soybean meal

18.00

19.00.

21.00

23.00

23.00

Corn bran

8.00

10.00

8.00

5.00

7.00

Premix1

4.00

4.00

4.00

4.00

5.00

Total

100.00

100.00

100

100

100.00

DE (MJ/ kg)

11.64

9.73

8.40

6.84

12.44

CP (%)

15.19

15.22

15.19

15.17

13.82

EP

TE D

M AN U

SC

RI PT

26

Control, maintenance energy; ER1, 85% maintenance energy; ER2, 70% maintenance energy;

1

AC C

ER3, 55% maintenance energy.

The premix provided the following per kg of diets: Fe 56 mg, Cu 15 mg, Mn 30 mg, Zn 40 mg, I

1.5 mg, Se 0.5 mg, Co 0.25 mg, S 3.2 g, VA 2150 IU, VD 170 IU, VE 131 IU, super concentrated YuanKangbao 2.7 g, 2% monensin l.6 g, sodium sulfate 10.1 g. All ingredients are calculated according to the raw material. Except the DE value, the rest of the ingredients are displayed as percentage.

ACCEPTED MANUSCRIPT Table. 2. Effect of restricted feeding on Hu sheep body weight and testes (n = 10) Control

Items

Mean

SEM

ER1

Mean

SEM

ER2

Mean

ER3

SEM

Mean

P-Value

0.60

<0.01

slaughter weight/kg

31.60a

0.67

25.83b

0.38

25.59b

0.42

Testicular weight/g

250.83a

23.61

181.43b

9.91

145.83bc

2.34

97.77c

21.80

<0.01

Testes ratio /%

0.792a

0.067

0.704a

0.047

0.571b

0.016

0.436b

0.095

0.025

RI PT

22.29c

SEM

SC

Control, maintenance energy; ER1, 85% maintenance energy; ER2, 70% maintenance energy;

a,b,c

M AN U

ER3, 55% maintenance energy.

Mean values with different superscript letters within a row are significantly different (P <

AC C

EP

TE D

0.05)

ACCEPTED MANUSCRIPT

Table. 3. Effect of energy restricted feeding on reproduction hormone in Hu sheep blood (n=10) ER1

Mean

SEM

SEM

Mean

LH (mIU/mL)

45.09

5.04

46.89

3.83

43.49

FSH (mIU/mL)

78.64

3.95

75.95

5.61

76.94

Testosterone (nmol/L)

15.61a

0.51

15.85a

0.50

16.24a

SEM

Mean

SEM

42.13

5.29

0.383

3.95

72.5

6.12

0.547

1.46

12.14b

0.56

0.020

Control, maintenance energy; ER1, 85% maintenance energy; ER2, 70% maintenance energy; ER3, 55% maintenance energy; FSH, follicle stimulating hormone; LH, luteinizing hormone.

EP

TE D

Mean values with different superscript letters within a row are significantly different (P < 0.05)

AC C

a,b

P-Value

2.59

SC

Items

ER3

M AN U

Mean

ER2

RI PT

Control

ACCEPTED MANUSCRIPT

Table. 4. Effect of energy compensation on Hu sheep body weight and testes (n = 10)

Items

Mean

SEM

ER1-C

Mean

SEM

ER2-C

Mean

ER3-C

RI PT

Control-C

SEM

Mean

SEM

P-Value

44.28

0.65

44.01

0.96

43.20

1.72

42.00

1.22

0.557

Testicular weight/g

265.15

27.39

287.45

14.98

330.20

27.85

346.13

27.85

0.184

Testes ratio /%

0.60

0.065

0.65

0.044

0.76

0.066

0.83

0.102

0.157

M AN U

SC

slaughter weight/kg

Control-C, 100% energy requirements prior to the compensation period; ER1-C, the 85% energy group with energy compensation; ER2-C, the 70% energy group with energy compensation;

AC C

EP

TE D

ER3-C, the 55% energy group with energy compensation.

ACCEPTED MANUSCRIPT

Table. 5. Effect of energy compensation on reproduction hormone in Hu sheep blood (n =10) ER1-C

Mean

SEM

SEM

Mean

LH (mIU/mL)

43.50

3.31

45.62

4.22

42.72

FSH (mIU/mL)

76.91

8.52

76.16

6.45

78.69

Testosterone (nmol/L)

18.42

1.46

18.55

0.62

17.47

SEM

Mean

SEM

42.05

3.60

0.934

5.79

74.86

8.86

0.718

0.20

16.20

1.95

0.329

Control-C, 100% energy requirements prior to the compensation period; ER1-C, the 85% energy group with energy compensation; ER2-C, the 70% energy group with energy compensation; ER3-C, the 55% energy group with energy compensation; FSH, follicle stimulating hormone; LH,

AC C

EP

TE D

luteinizing hormone.

P-Value

5.16

SC

Items

ER3-C

M AN U

Mean

ER2-C

RI PT

Control-C

ACCEPTED MANUSCRIPT

Table. 6. The DEGs related to Hu sheep testes development through KEGG classification.

KEGG pathway

CYP11A1

Ovarian steroidogenesis

KEGG ID

RI PT

Reproductive related gene

oas04913

Steroid hormone biosynthesis

oas00140

Drug metabolism - cytochrome P450

oas00982

DUSP5

Mitogen-activated protein kinase signaling

oas04010

M AN U

SC

ALDH3B1

pathway

PI3K-Akt signaling pathway

oas04151

Adherens junction

oas04520

Ovarian steroidogenesis

oas04913

TE D

INSR

AMPK signaling pathway

oas04152

Wnt signaling pathway

oas04310

Estrogen signaling pathway

oas04915

Oxytocin signaling pathway

oas04921

FDFT1

Steroid biosynthesis

oas00100

GUCY1A3

Oxytocin signaling pathway

oas04921

ATP1A2

Cyclic adenosine monophosphate

oas04024

WNT2

AC C

EP

KCNJ6(GIRK2)

signaling pathway PGR (NR3C3)

Progesterone-mediated oocyte maturation

oas04914

DEGs, differentially expressed genes; KEGG, kyoto encyclopedia of genes and genomes;

ACCEPTED MANUSCRIPT CYP11A1, cytochrome P450 family 11 subfamily A member 1; ALDH3B1, aldehyde dehydrogenase 3 family member B1; DUSP5, dual specificity phosphatase 5; INSR, insulin receptor ; Wnt2, Wnt family member 2; KCNJ6(GIRK2), potassium voltage-gated channel

RI PT

subfamily J member 6; FDFT1, farnesyl-diphosphate farnesyltransferase 1; GUCY1A3, guanylate cyclase 1 soluble subunit alpha; ATP1A2, ATPase Na+/K+ transporting subunit alpha 2;

AC C

EP

TE D

M AN U

SC

PGR(NR3C3), progesterone receptor.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT Revised

Highlights 1. Dietary energy levels affect sheep testicular weight and reproductive hormone levels.

RI PT

2. Dietary energy restriction affect sheep testicular transcriptional levels.

3. 81 up-regulated and 180-downregulated genes were identified in ER group compared

AC C

EP

TE D

M AN U

SC

to the control.

Revised

ACCEPTED MANUSCRIPT

Transcriptome profiling of Hu sheep testicular development associated with dietary energy restriction and during subsequent compensatory feeding

TE D

M AN U

SC

RI PT

Y.X. Fan, Z. Wang, T.W. Ma, K.P. Deng, X. Feng, F.Z. Li, F. Wang, and Y.L. Zhang*

AC C

EP

Supplementary Figure S1 The number of reads mapped to the chromosome of Hu sheep fed on diets containing 55% (ER3) and 100% (Control) of energy requirements. The figure a-f represent control group (A13, A15, A16), ER3 (D9, D10, D11) group, respectively. The abscissa is the length of the chromosome information (in units of Mb), the ordinate is the number of reads mapped to the chromosome (in units of K).

ACCEPTED MANUSCRIPT Supplementary Figure S2. Pearson correlation between the testicles samples of Hu sheep fed on diets containing 55% (ER3) and 100% (control) of energy requirements. High-quality samples should have values higher than 0.800: the higher the value, the stronger is the correlation. The Control group include three replicas: A13, A15, A16 and ER3 group include three duplications: D9, D10, D11.

Gene

NCBI Reference Sequence

INSR

XM_012177948.2

CYP11A1

NM_001093789.1

ALDH3B1

XM_012117488.2 GAAI01002542.1

FDFT1

XM_012129046.2

Wnt2

NM_001195319.1

PGR

XM_012169085.2

GAPDH

NM_001190390.1

RI PT

Supplementary Table S1 Details of primer sequences, expected product sizes and GenBank accession numbers of genes used for qRT– PCR. Primer sequence (5’-3’)

Product size (bp)

F-GTACCCTGGAGAAGTGTGCC R-CCCGGAAGAGCAGCAAGTAA F-CCATACCAGTGTCCCTCTGC R-TGGCCTTGACATCCTCCAAG F-TCCGAGATCAGCAAGAGCAC R-GCCCGCCCAGCACCA F-GTGCTATTCCACAGATTTATCACAG R-TAGATGGGCGAGTAGTGGCT F-TGACTTCAGGAAAACGGGCA R-CTTCCGGGTAATGTGGGAGG F-AACTACCTGAGGCCCGATTC R-TCCGAAAACCTCCAAGAACCA F-GTCAAGGCAGAGAACGGGAA R-GGTTCACGCCCATCACAAAC

SC

191 230

M AN U

91 139 289 288

TE D

232

EP

INSR, Insulin receptor ; CYP11A, cytochrome P450 family 11 subfamily A member 1; ALDH3B1, aldehyde dehydrogenase 3 family member B1; FDFT1, farnesyl-diphosphate farnesyltransferase 1; Wnt2, Wnt family member 2; PGR, progesterone receptor; GAPDH, glyceraldehyde 3-phosphate dehydrogenase

Group

Raw reads1

A13

54425742

Clean Reads2 52670316

A15

66943166

64746080

9.71G

0.01

97.52

94.12

49.94

A16

60417460

58355060

8.75G

0.01

97.45

93.99

50.19

D9

51866144

50163418

7.52G

0.01

97.46

93.98

49.85

D10

54903088

53090676

7.96 G

0.01

97.50

94.09

49.11

D11

50504544

48855200

7.33 G

0.01

97.51

94.12

49.57

Sample

Control

ER3

AC C

Supplementary Table S2 Sequencing data quality assessment of Hu sheep fed on diets containing 55% (ER3) and 100% (Control) of energy requirements Clean Base3 7.90G

Error rate (%)4 0.01

Q20(%)5

Q30(%)5

97.55

94.19

GC content(%)6 50.11

The Control group include three replicas: A13, A15, A16 and ER3 group include three duplications: D9, D10, D11. 1 Raw reads: Counts the original sequence data, counting the number of sequencing sequences for each file in four rows.

ACCEPTED MANUSCRIPT 2

RI PT

Clean reads: the calculation method with Raw Reads, but the statistical data for the filtered sequencing data. Subsequent bioinformatics analysis is based on Clean reads. 3 Clean bases: The number of sequencing sequences multiplied by the length of the sequencing sequence and converted to the unit of G. 4 Error rate: The average base sequence error rate. 5 Q20, Q30: Calculate the percentage of the base that the Phred value is greater than 20, 30, respectively, of the total base. 6 GC content: The sum of the total number of bases G and C in the total number of bases. Supplementary Table S3 Reads compared with reference genome of Hu sheep fed on diets containing 55% (ER3) and 100% (Control) of energy requirements

Uniquely mapped4 Read-1 Read-2 Reads map to '+'5 Reads map to '-'5 Non-splice reads

52670316

64746080

58355060

41436772 (78.67%) 975201 (1.85%) 40461571 (76.82%) 20837313 (39.56%) 19624258 (37.26%) 20265469 (38.48%) 20196102 (38.34%) 24415317 (46.35%) 16046254 (30.47%)

50725179 (78.34%) 1092046 (1.69%) 49633133 (76.66%) 25611871 (39.56%) 24021262 (37.1%) 24888259 (38.44%) 24744874 (38.22%) 28564390 (44.12%) 21068743 (32.54%)

45301748 (77.63%) 999695 (1.71%) 44302053 (75.92%) 22890964 (39.23%) 21411089 (36.69%) 22191581 (38.03%) 22110472 (37.89%) 25555135 (43.79%) 18746918 (32.13%)

AC C

Splice reads6

A16

D9

D10

D11

50163418

53090676

48855200

39760251 (79.26%) 832313 (1.66%) 38927938 (77.6%) 20119900 (40.11%) 18808038 (37.49%) 19510573 (38.89%) 19417365 (38.71%) 22929110 (45.71%) 15998828 (31.89%)

41973212 (79.06%) 1012044 (1.91%) 40961168 (77.15%) 21147373 (39.83%) 19813795 (37.32%) 20507957 (38.63%) 20453211 (38.53%) 25813371 (48.62%) 15147797 (28.53%)

38349911 (78.5%) 858434 (1.76%) 37491477 (76.74%) 19351813 (39.61%) 18139664 (37.13%) 18743151 (38.36%) 18748326 (38.38%) 23057342 (47.2%) 14434135 (29.54%)

M AN U

Multiple mapped3

A15

TE D

Total mapped

2

A13

EP

Total reads

1

ER3

SC

Control

Sample name

The Control group include three replicas: A13, A15, A16 and ER3 group include three duplications: D9, D10, D11. 1 Total reads: The number of sorted data after the sequencing of the data. 2 Total mapped: the number of sequencing sequences that can be positioned on the genome; 3 Multiple mapped: The number of sequencing sequences with multiple alignment positions on the reference sequence; 4 Uniquely mapped: The number of sequencing sequences with unique alignment positions on the reference sequence. 5 Reads map to '+', Reads map to '-': The sequencing sequence compares to the positive and negative chains of the genome. 6 Splice reads: The statistics of the sequencing sequence (also known as Junction reads) on the two exons are compared in total mapped, and the non-splice reads are the statistics of the

ACCEPTED MANUSCRIPT entire sequence of the sequences to the exon. The percentage of Splice reads depends on the length of the sequencing clip.

Supplementary Table S4 Number of genes with different expression levels of Hu sheep fed on diets containing 55% (ER3) and 100% (Control) of energy requirements

3-15 16-60 >60

A13

A15

A16

D9

10177 (34.73%) 4535 (15.47%) 6676 (22.78%) 5390 (18.39%) 2529 (8.63%)

10803 (36.86%) 5057 (17.26%) 6721 (22.93%) 4206 (14.35%) 2520 (8.60%)

10079 (34.39%) 5132 (17.51%) 6996 (23.87%) 4563 (15.57%) 2537 (8.66%)

10084 (34.41%) 5135 (17.52%) 7032 (23.99%) 4530 (15.46%) 2526 (8.62%)

D10

D11

9918 (33.84%) 4770 (16.28%) 6717 (22.92%) 5400 (18.43%) 2502 (8.54%)

10200 (34.80%) 4765 (16.26%) 6364 (21.71%) 5419 (18.49%) 2559 (8.73%)

RI PT

1-3

ER3

SC

0-1

Control

M AN U

FPKM Interval

name SNP InDel

A13 435160 38188

Control A15 457986 39183

TE D

The Control group include three replicas: A13, A15, A16 and ER3 group include three replicas: D9, D10, D11. The number of genes in different expression levels was calculated by counting the number of genes at different expression levels and the level of expression of individual genes. Supplementary Table S5 Number of SNP and InDel of Hu sheep fed on diets containing 55% (ER3) and 100% (Control) of energy requirements A16 447849 38482

Average 446998 38618

D10 478411 41982

ER3 D11 452898 39622

D9 432531 37326

Average 454613 39643

AC C

EP

The Control group include three replicas: A13, A15, A16 and ER3 group include three replicas: D9, D10, D11; SNP, Single Nucleotide Polymorphisms; InDel, Insert and Deletion.