MicroRNA-106b Regulates Milk Fat Metabolism via ATP Binding

Subscriber access provided by ALBRIGHT COLLEGE

Biotechnology and Biological Transformations

microRNA-106b regulates milk fat metabolism via ATP binding cassette subfamily A member 1 (ABCA1) in bovine mammary epithelial cells Zhi Chen, Shuangfeng Chu, Xiaolong Wang, Yongliang Fan, Tiayin Zhan, Abdelaziz Adam Idriss Arbab, Mingxun Li, Huimin Zhang, Yongjiang Mao, Juan J. Loor, and Zhangping Yang J. Agric. Food Chem., Just Accepted Manuscript • DOI: 10.1021/acs.jafc.9b00622 • Publication Date (Web): 20 Mar 2019 Downloaded from http://pubs.acs.org on March 21, 2019

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 28

Journal of Agricultural and Food Chemistry

1

Full title: microRNA-106b regulates milk fat metabolism

2

via ATP binding cassette subfamily A member 1 (ABCA1) in

3

bovine mammary epithelial cells

4

Authors: Zhi Chen,

5

Zhan,§ Abdelaziz Adam Idriss Arbab,

6

Mao,

7

Corresponding Author: Zhangping Yang

8

Author to whom correspondence should be addressed: E-Mail: [email protected];

9

Tel:+86-051S44-87979269

†‡

†‡

Shuangfeng Chu, †‡

†‡

Xiaolong Wang, Mingxun Li,

†‡

†‡

Yongliang Fan,

Huimin Zhang,

†‡

†‡

Tiayin

Yongjiang

Juan J. Loor, Zhangping Yang*†‡ ⊥

10

Affiliation: 1 †College of Animal Science and Technology, Yangzhou University,

11

Yangzhou 225009, PR China;

12

‡

13

Ministry of Education, Yangzhou University, Yangzhou 225009, China

14

§

15

Science and Technology, Northwest A&F University, Yangling, Shaanxi 712100, PR

16

China

17

⊥

18

Division of Nutritional Sciences, University of Illinois, Urbana, IL 61801, USA

Joint International Research Laboratory of Agriculture & Agri-Product Safety,

Shanxi Key Laboratory of Molecular Biology for Agriculture, College of Animal

Mammalian Nutrition Physiology Genomics, Department of Animal Sciences and

19 20 21 22 23 24 25 26 27 1

ACS Paragon Plus Environment


Page 2 of 28

28

Abstract

29

Research on the mechanisms regulating milk fat synthesis in dairy cows is essential in

30

the context of identifying potential molecular targets that in the long-term can help

31

develop appropriate molecular breeding programs. Although some studies have

32

revealed that microRNA (miRNA) affect lipid metabolism by targeting specific genes,

33

joint analysis of miRNA and target mRNA data from bovine mammary tissue has

34

revealed few clues regarding the underlying mechanisms controlling milk fat synthesis.

35

The objective of the present study was to use high-throughput sequencing and

36

bioinformatics analysis to identify miRNA and mRNA pairs and explore further their

37

potential roles in regulating milk fat synthesis. A total of 233 pairs of negatively-

38

associated miRNA and mRNA pairs were detected. Among those, there were 162 pairs

39

in which the miRNAs were down-regulated and the target mRNAs up-regulated.

40

Among the identified miRNA, miR-106b can bind the 3′-UTR of ATP binding cassette

41

subfamily A member 1 (ABCA1), a gene previously identified as having a positive

42

association with bovine milk fat synthesis. Overexpressing miR-106b in bovine

43

mammary epithelial cells decreased while inhibiting miR-106b increased triglyceride

44

and cholesterol content, confirming its role in lipid metabolism. The present study

45

allowed for constructing a miR-106b-ABCA1 regulatory network map, hence,

46

providing a theoretical basis for targeting this gene in molecular breeding of dairy cows.

47

Key Words: miR-106b, milk fat metabolism, ABCA1, Bovine mammary epithelial

48

cells

49 50

INTRODUCTION

51

Bovine milk contains a variety of biologically-active fatty acids that have the potential

52

to affect nutrient metabolism during growth and development in humans1. Therefore,

53

studying the regulation of milk fat synthesis and fatty acid composition remains an

54

active area of research2. In the past few decades, studies on bovine milk fat synthesis

55

regulation have focused on single group analysis or gene function verification without

56

a more comprehensive approach using “multiomics” tools that can accelerate

57

understanding of the underlying molecular regulatory mechanisms3.

58

Following the completion of the human genome project, the function and regulatory

59

mechanisms of non-coding RNA began to be examined4. miRNAs, a type of non-coding

60

RNA, are now well-known to play an important role in regulating transcription of

61

protein-coding genes5. Mammary fatty acid metabolism involves the expression, 2


Page 3 of 28


62

network regulation, and signal transduction of multiple genes (including miRNAs).

63

Although the first gene regulatory networks controlling milk fat synthesis were reported

64

10 years ago, the recent application of high-throughput sequencing technology has

65

underscored that a large number of regulatory factors and mechanisms remain to be

66

identified6.

67

The main objective of the present study was to use high-throughput sequencing data

68

from bovine mammary tissue in the non-lactating (“dry” period) and lactating periods

69

to determine expression profiles of microRNA (miRNA) and mRNA. Those data were

70

used to 1) establish a miRNA and target mRNA regulatory network map, 2) establish a

71

network of milk-related genes, and 3) elucidate milk fat regulation pathways relevant

72

to dairy cattle lactation. In vitro studies with isolated bovine mammary epithelial cells

73

allowed for more detailed mechanistic evaluation miRNA and target mRNA

74

interactions in the context of lipid metabolism. Among the key findings, the analysis

75

highlighted the existence of a regulatory pathway between miR-106b and ATP binding

76

cassette transporter A1 (ABCA1).

77 78

MATERIALS AND METHODS

79

Animals, sampling, and RNA extraction

80

Holstein cows used in this research were from the elite herd from the experimental farm

81

at Yangzhou University, Yangzhou, Jiangsu Province, China. Three Holstein cows of

82

similar body weight were biopsied during the non-lactating (“dry” period) and early

83

lactation. Surgical methods were used to collect 1–2 g of mammary tissue per cow.

84

Tissue was washed with diethyl pyrocarbonate-treated water and stored immediately in

85

liquid nitrogen. Total RNA was extracted from tissue using Trizol reagent (Invitrogen,

86

Carlsbad, CA). The RNA integrity number method was used to detect RNA quality,

87

and only good-quality RNA was used for the study7, 8.

88

miRNA sequencing analysis

89

Construction and sequencing of small RNA libraries

90

The RNA libraries began to be constructed from RNA and included the addition of a 3′

91

and 5′ linker, RT-PCR enrichment, cDNA library construction, and purification by

92

polyacrylamide gel electrophoresis. Library quantification (using the Qubit machine

93

quantification miRNA library) and sequencing (using an Illumina HiSeq. 2500) were

94

conducted as described previously9.

95

Differential expression analysis of miRNA 3



Page 4 of 28

96

miRNA expression was calculated as transcripts per million (TPM). The TPM formula

97

was as follows: (number of reads per miRNA alignment)/(number of reads from the

98

total sample alignment)  106. The number of miRNA expressed per million alignments

99

was used as an indicator. The DEseq approach was used to calculate differentially

100

expressed miRNAs in each sample. Differential miRNAs were defined as those with a

101

p-value < 0.05 and a difference ≥ 2-fold10.

102

Transcriptome sequencing analysis

103

Digital gene expression library preparation

104

After total RNA extraction, mRNA was enriched and purified using magnetic beads

105

coated with Oligo (dT), followed by random fragmentation of the mRNA using a

106

fragmentation reagent. The experiment was completed in the preparation of the entire

107

library. The constructed library was sequenced with an Illumina Hiseq 2500 at a

108

sequence read length of 1 × 36 bp11.

109

Data quality control

110

After raw data of the second-generation sequencing was obtained, the quality of the

111

original

112

(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

113

Data processing

114

Unique sequence data were obtained and the corresponding copy number determined.

115

A chi-square test was used to screen differentially expressed genes using the differential

116

test method. Genes with p < 0.05 and false discovery rate (FDR) < 0.05 were defined

117

as differentially expressed genes. The selected differentially expressed genes were

118

further analyzed for clustering, gene ontology (GO) enrichment, and Kyoto

119

Encyclopedia of Genes and Genomes (KEGG) enrichment.

120

Construction of miRNA and target mRNA regulatory networks

121

The miRanda software was used to conduct analysis of differentially expressed miRNA

122

and mRNA. From this analysis, we focused on negative regulatory relationships

123

between miRNA and specific transcripts followed by construction of a negative

124

regulation microRNA network map12.

125

Cell culture and transfection

126

Resuscitated dairy bovine mammary epithelial cells were cultured at 37°C in 5% CO2

127

and appropriate humidity conditions. After the cells were grown to approximately 80%

128

in culture dishes, they were infected with small RNA chemical synthesis reagents, and

reads

was

examined

using

the

FastQC

online

tool

4


Page 5 of 28


129

the cells were collected after 48 h for subsequent index detection. Three biological

130

replicates were conducted for each treatment.

131

Triglyceride and cholesterol content assays

132

Six-well cell culture plates were used to grow and process cells; analysis were

133

performed in triplicate. After adding 150 μL of lysate and resting at 4°C for 15 min, the

134

cells and the lysate were collected in 1.5-mL tubes. Cells were then homogenized while

135

on ice, and the supernatant collected after centrifugation. A total of 10 μL of supernatant

136

was used for protein concentration using the bicinchoninic acid method13, 14.

137

RT-qPCR and western blotting

138

The mature miRNA expression level was determined using the S-Poly (T) assay using

139

the S3 primer as specific reverse primer. The 18S ribosomal RNA (rRNA) was used as

140

internal control. Briefly, reverse transcription was performed as follows: the 10 µL

141

reaction included 0.2 µg total RNA, 2.5 µL of 4×reaction buffer, 1 µL of poly A/RT

142

enzyme mix, 1 µL of 0.5 µM RT primer. The reaction was performed at 37 °C for 30

143

min, followed by 42 °C for 30 min, then 75 °C for 5 min.

144

The mRNA expression level was determined using RT-qPCR. A total of 0.5 µg of RNA

145

was used for cDNA synthesis using the Prime Script® RT Reagent Kit (Perfect Real

146

time, Takara Bio Inc., Kusatsu, Japan). The expression of target genes was normalized

147

to β -actin. The primer sequences are listed in Supplemental file S4. Relative

148

expression was calculated using the 2-△△Ct method15-18.

149

Proteins extracted from cells were separated by sodium dodecyl sulfate polyacrylamide

150

gel electrophoresis, transferred to a nitrocellulose membrane (Millipore, Burlington,

151

MA) and probed with primary monoclonal rabbit anti-ABCA1 (#ab18180, Abcam,

152

Cambridge, UK), rabbit anti-β-casein (51067-2-AP, Proteintech Group, Wuhan, China)

153

and monoclonal mouse anti-β-actin (66009-1-IG, Proteintech Group) 19-21.

154

Luciferase reporter assay

155

Primers for the 3′-untranslated region (UTR) of the ABCA1 gene containing the miR-

156

106b site were used for amplification. The bovine ABCA1 gene 3′-UTR was used as

157

template. The 3′-UTR fragment of the ABCA1 gene and the psiCHECK-2 luciferase

158

vector obtained by PCR amplification were digested by XhoI and NotI (S5) and ligated

159

overnight at 4°C using T4 ligase. The 3′-UTR of the ABCA1 gene was cloned into the

160

psiCHECK-2 luciferase vector and verified by sequencing in the next step. Bovine

161

mammary epithelial cells (BMECs) were plated in 48-well plates and transfection was 5



Page 6 of 28

162

performed using the Roche Transfection Reagent X-treme GENE HP DNA

163

Transfection Reagent (Roche, Mannheim, Germany) when the cell density reached

164

approximately 50,000 cells per well. Fluorescence intensity measurements were

165

conducted using the Dual-Glo luciferase assay system kit22.

166

Statistical analysis

167

SPSS 18.0 software (SPSS Inc., Chicago, IL) was used for statistical analysis. The mean

168

value of the continuous variable application was expressed as standard deviation. One-

169

way analysis of variance was performed to determine significance at p < 0.05.

170 171

RESULTS

172

miRNA sequencing analysis

173

Rfam comparison filtering

174

The clean read sequences were compared with data in the Rfam (version 10.0) database

175

using blastn software. Results with an E-value ≤ 0.01 were extracted, and the sequences

176

of rRNA, small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), and transfer

177

RNA (tRNA) annotated. We attempted to find and remove possible rRNAs, small

178

conditional RNAs, snoRNAs, snRNAs, and tRNAs and to perform small RNA removal

179

(Table 1).

180

Repeat sequence filtering

181

Using RepeatMasker (version 4.0.1) software, the filtered mRNA sequences of

182

degradation fragments were aligned with the repeat database to identify possible repeat

183

sequences. Repeat sequences were not removed prior to performing the alignments of

184

known miRNA and new miRNA predictions (Fig. S1).

185

Known miRNA base preference statistics

186

Development of a miRNA from a precursor to a mature molecule is performed by

187

digestion with Dicer. In general, the specificity of the restriction site allows the first

188

base of the mature miRNA sequence to have a strong bias toward the base U. There are

189

also some known statistical data for other sites, such as that numbers 2–4 are generally

190

missing U, and that number 10 is biased towards A (Fig. S2).

191

miRNA expression analysis

192

We estimated the expression level of miRNAs using the small RNA sequencing

193

analysis. This involved determining the mature sequence of the species and counting

194

the newly predicted miRNA sequence (Figs. 1, S1).

195

Transcriptome sequencing analysis 6


Page 7 of 28


196

Differential expression analysis

197

When using RNA sequencing data to compare and analyze whether there is differential

198

expression of the same gene in two samples, two criteria can be used: fold-change,

199

which is the multiple of the change of the expression level of the same gene in two

200

samples, and the p-value or False discovery rate (FDR adjusted p-value). The

201

calculation of FDR values requires that the p-value of each gene be calculated first. The

202

default screening criteria were p < 0.05 and a fold-change greater than 2-fold. We

203

performed mapping analysis for the differentially expressed genes that were screened

204

out (Figs. 2, S2).

205

Gene ontology (GO) enrichment analysis of differentially expressed genes

206

After differentially expressed genes were obtained, we performed GO enrichment

207

analysis and tabulated their functions (Fig. 3). The number of differentially expressed

208

genes included in each GO item was counted, and the significance of enrichment in

209

each GO term was calculated using the hypergeometric distribution test method. The

210

GO functional enrichment analysis method was as follows: the whole protein coding

211

gene/transcript was the background list, and the differential protein coding

212

gene/transcript list was the candidate list. The hypergeometric distribution test was used

213

to calculate a p-value that represented whether the GO function set was significantly

214

enriched among the differentially expressed protein-encoding gene/transcript list. The

215

P value was corrected by the Benjamini & Hochberg multiple test to obtain the FDR.

216

The calculation returned an enriched p-value; a small p-value indicating that the

217

differentially expressed gene was enriched in the GO term.

218

The hypergeometric distribution test calculates the p-value as follows:

219

220 221

The enrichment score calculation formula is:

222 223

where N is the number of genes with GO annotation, n is the number of GO annotated

224

genes among differentially expressed genes in N, M is the number of genes annotated 7



Page 8 of 28

225

to a particular GO term, and m is the number of differentially expressed genes annotated

226

to a particular GO term.

227

Analysis of differential gene KEGG enrichment

228

The pathway analysis of differentially expressed genes was conducted using the KEGG

229

database, a public database of manually-curated biological pathways. The significance

230

of differential gene enrichment in each pathway entry was calculated using a

231

hypergeometric distribution test. The result of the calculation returned a p-value

232

denoting the significance of the enrichment, with a small p-value indicating that the

233

differentially expressed gene was enriched in the pathway (Fig. 4).

234

Association analysis of miRNA and target mRNA expression

235

Correlation analysis of mRNA and miRNA sequencing results resulted in 233 of

236

negatively correlated miRNA and target mRNA pairs. Among them, miRNAs were

237

down-regulated and mRNAs were up-regulated in 162 pairs. The most-enriched

238

miRNA were miR-382 and miR-424-5p with 12 target genes. Ten differentially

239

expressed miRNAs (miR-223, miR-184, miR-132, miR-1246, miR-130, miR-196a,

240

miR-205, miR-200b, miR-31, and miR-145) were selected for association with 48

241

target genes (Figs. 5, S3). These differential miRNAs and target mRNA may play an

242

important role in milk fat synthesis regulation in dairy cows.

243

miR-106b specific targeting of ABCA1

244

When we selected candidate genes such as miR-106b and ABCA1, we attempted to

245

verify their regulatory relationships experimentally. The approach to target a given

246

miRNA for gene prediction was based on potential 3′-UTR interactions with mRNA

247

using software packages such as TargetScan and DAVID (https://david.ncifcrf.gov).

248

Analysis indicated that ABCA1 has a miR-106b binding site in its 3′-UTR, suggesting

249

it may be a target gene of this miRNA. The focus on ABCA1 for functional validation

250

was due to its known known association with mammary lipid metabolism. The protein

251

encoded by ABCA1 promotes intracellular free cholesterol and phospholipid efflux by

252

hydrolyzing ATP23-25. Both mRNA and protein results indicated that ABCA1 expression

253

is up-regulated by inhibition of miR-106b, while over-expression of miR-106b leads to

254

down-regulation of ABCA1 expression (Fig. 6A, 6D). To demonstrate that miR-106b

255

directly targets ABCA1, we synthesized a 3′-UTR fragment containing the miR-106b

256

target site of ABCA1, and cloned it into a psi-CHECK2 vector to construct a 3′-UTR

257

reporter plasmid. A luciferase reporter assay indicated that miR-106b overexpression

258

decreased luciferase activity of the wild-type 3′-UTR reporter. However, there was no 8


Page 9 of 28


259

significant difference in the luciferase activity of the mutant reporter gene (Fig. 6B,

260

6C). These results suggest that miR-106b directly targets the mRNA site of ABCA1 and

261

plays a negative regulatory role.

262

Functional evaluation of the miR-106b and ABCA1 association in BMECs

263

Expression of miR-106b decreases triglyceride and cholesterol content in BMECs

264

After determining that ABCA1 is a target gene of miR-106b, we sought to uncover its

265

specific function in BMECs. Compared with the negative control, the expression of

266

miR-106b was 48 times greater in BMECs overexpressing the miRNA. In contrast,

267

compared with the negative control, the expression of miR-106b more than 99% lower

268

in BMECs in which this miRNA was silenced (Fig, S3). Milk fat is composed mainly

269

of triglycerides (TAGs). Therefore, we examined the levels of TAG and cholesterol

270

during overexpression and suppression of miR-106b, respectively. Compared with

271

negative controls, the miR-106b mimic reduced the levels of TAG to 0.23-fold (p

MicroRNA-106b Regulates Milk Fat Metabolism via ATP Binding

Recommend Documents