This is the accepted version of the journal article: Mármol-Sánchez, Emilio; Cirera, Susanna; Zingaretti, Laura M.; [et al.]. «Mod- eling microRNA-driven post-transcriptional regulation using exon-intron split analysis in pigs». Animal Genetics, Vol. 53, Issue 5 (October 2022), p. 613-626. DOI 10.1111/age.13238 This version is available at https://ddd.uab.cat/record/273541 under the terms of the license 1 Running title: Study of post-transcriptional regulation in pigs 2 3 Modeling microRNA-driven post-transcriptional regulation by 4 using exon-intron split analysis (EISA) in pigs 5 Emilio Mármol-Sánchez1*, Susanna Cirera2, Laura M. Zingaretti3, Mette Juul Jacobsen2, 6 Yuliaxis Ramayo-Caldas4, Claus B. Jørgensen2, Merete Fredholm2, Tainã Figueiredo 7 Cardoso1†, Raquel Quintanilla4, Marcel Amills1,5 8 1 9 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, 10 Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain. 2Department of Veterinary 11 and Animal Sciences, Faculty of Health and Medical Sciences, University of 12 Copenhagen, 1871 Frederiksberg C, Denmark. 3Universidad Nacional de Villa María, 13 Villa María, Córdoba, Argentina. 4Animal Breeding and Genetics Program, Institute for 14 Research and Technology in Food and Agriculture (IRTA), Torre Marimon, 08140 15 Caldes de Montbui, Barcelona, Spain. 5Departament de Ciència Animal i dels Aliments, 16 Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain. 17 18 *Emilio Mármol-Sánchez current affiliation: 1Department of Molecular Biosciences, The 2 19 Wenner-Gren Institute, Stockholm University, Stockholm, Sweden. Centre for 20 Paleogenetics, Stockholm University, Stockholm, Sweden. 21 † 22 Tainã Figueiredo Cardoso current affiliation: Embrapa Pecuária Sudeste, Empresa 23 Brasileira de Pesquisa Agropecuária (EMBRAPA), 13560-970, São Carlos, SP, Brazil. 24 1 25 Corresponding author: Emilio Mármol-Sánchez. Science for Life Laboratory, 26 Department of Molecular Biosciences, The Wenner-Gren Institute. Stockholm 27 University, Stockholm, Sweden. Email: emilio.marmol.sanchez@gmail.com 28 29 Emilio Mármol-Sánchez: emilio.marmol.sanchez@gmail.com 30 Susanna Cirera: scs@sund.ku.dk 31 Laura M. Zingaretti: m.lau.zingaretti@gmail.com 32 Mette Juul Jacobsen: pmn418@alumni.ku.dk 33 Yuliaxis Ramayo-Caldas: yuliaxis@gmail.com; yuliaxis.ramayo@irta.cat 34 Claus B. Jørgensen: clausbj@sund.ku.dk 35 Merete Fredholm: mf@sund.ku.dk 36 Tainã Figueiredo Cardoso: tainafcardoso@gmail.com 37 Raquel Quintanilla: raquel.quintanilla@irta.cat 38 Marcel Amills: Marcel.Amills@uab.cat 39 40 41 42 43 44 45 46 47 48 2 49 Abstract 50 The contribution of microRNAs (miRNAs) to mRNA post-transcriptional regulation has 51 often been explored by post hoc selection of downregulated genes and determining 52 whether they harbor binding sites for miRNAs of interest. This approach, however, does 53 not discriminate whether these mRNAs are also downregulated at the transcriptional 54 level. Here, we have characterized the transcriptional and post-transcriptional changes of 55 mRNA expression in two porcine tissues: gluteus medius muscle of fasted and fed Duroc 56 gilts and adipose tissue of lean and obese Duroc-Göttingen minipigs. Exon-intron split 57 analysis (EISA) of RNA-seq data allowed us to identify downregulated mRNAs with high 58 post-transcriptional signals in fed or obese states, and we assessed whether they harbor 59 binding sites for upregulated miRNAs in any of these two physiological states. We found 60 26 downregulated mRNAs with high post-transcriptional signals in the muscle of fed gilts 61 and 21 of these were predicted targets of upregulated miRNAs also in fed pigs. For 62 adipose tissue, 44 downregulated mRNAs in obese minipigs displayed high post- 63 transcriptional signals, and 25 of these were predicted targets of miRNAs upregulated in 64 the obese state. These results suggest that the contribution of miRNAs to mRNA 65 repression is more prominent in the skeletal muscle system. Finally, we identified several 66 genes that may play relevant roles in the energy homeostasis of the pig skeletal muscle 67 (DKK2 and PDK4) and adipose (SESN3 and ESRRG) tissues. By differentiating 68 transcriptional from post-transcriptional changes in mRNA expression, EISA provides a 69 valuable view about the regulation of gene expression, complementary to canonical 70 differential expression analyses. 71 72 Keywords: Exon-intron split analysis, microRNA, pigs, energy homeostasis. 73 3 74 Introduction 75 The post-transcriptional regulation of gene expression plays a fundamental role towards 76 shaping fine-tuned biological responses to environmental changes (Schaefke et al. 2018). 77 Such regulation can take place at multiple levels including splicing, 3’-cleavage and 78 polyadenylation, decay or translation, and its main effectors are RNA binding proteins 79 and non-coding RNAs (Schaefke et al. 2018). Of particular importance are microRNAs 80 (miRNAs), which are primarily engaged in the post-transcriptional control of gene 81 expression through inhibition of translation and/or degradation of target mRNAs (Bartel, 82 2018). 83 Multiple differential expression studies have been performed in pigs during the last 84 decade (Pérez-Montarelo et al. 2013; Óvilo et al. 2014; Pilcher et al. 2015; Horodyska et 85 al. 2018; Benítez et al. 2019). One of the main limitations of these studies is that the 86 transcriptional and post-transcriptional components of gene regulation are not 87 independently analyzed. This means that genes that are transcriptionally upregulated and 88 post-transcriptionally downregulated, or vice versa, might not be detected as significantly 89 differentially expressed. Another disadvantage of this approach is that it does not provide 90 insights about the causes of the observed downregulation of RNA transcripts. For 91 instance, studies have typically focused on specific sets of downregulated genes harboring 92 binding sites for miRNAs, in order to disentangle regulatory functions driven by miRNAs 93 (Han et al. 2017; Xie et al. 2019; Ali et al. 2021). This approach, however, cannot 94 distinguish between transcriptional or post-transcriptional repression, which is essential 95 to understand at which level of the mRNA life-cycle regulation is taking place. 96 To overcome this important limitation, Gaidatzis et al. (2015) devised the exon-intron 97 split analysis (EISA), which separates the transcriptional and post-transcriptional 98 components of gene regulation by comparing the amounts of exonic and intronic reads 4 99 from expressed mRNA transcripts. The main assumption of this method is that intronic 100 reads are predominantly derived from nascent unprocessed mRNAs or pre-mRNAs, so 101 they reflect transcriptional changes, while post-transcriptional changes can be inferred 102 from differences between the levels of the exonic and intronic fractions (Ameur et al. 103 2011; Gaidatzis et al. 2015). For instance, a gene showing little to no differences in the 104 number of sequenced intronic reads, but a strong downregulation of exonic reads after a 105 certain treatment or challenge (nutrition, infection, temperature etc.), could be subjected 106 to post-transcriptional repression (Gaidatzis et al. 2015; Cursons et al. 2018; Pillman et 107 al. 2019). Recent developments on this principle have also been applied to infer the 108 transcriptional fate of cells (La Manno et al. 2018). 109 The main goal of the present study was to investigate the contribution of miRNAs to post- 110 transcriptional regulatory responses using the EISA methodology, combined with in 111 silico prediction of miRNA-mRNA interactions and covariation analyses in porcine 112 skeletal muscle and adipose tissues. 113 114 115 Materials and methods 116 Experimental design, sampling and processing 117 Two experimental systems were used: 118 (i) Duroc pigs: Twenty-three gilts were subjected to two fasting/feeding regimes, i.e. 11 119 gilts (AL-T0) were slaughtered in fasting condition, while 12 gilts (AL-T2) were 120 slaughtered after 7 h of having access to food (Cardoso et al. 2017; Ballester et al. 2018; 121 Mármol-Sánchez et al. 2020). Immediately after slaughtering, gluteus medius (GM) 122 skeletal muscle samples were collected and snap-frozen at -80°C. 5 123 (ii) Duroc-Göttingen minipig F2 inter-cross: Ten individuals fed ad libitum with divergent 124 fatness profiles according to their body mass index (BMI, 5 lean and 5 obese) were 125 selected from the UNIK resource population (Kogelman et al. 2013; Jacobsen et al. 126 2019). Retroperitoneal adipose tissue was collected at slaughter and mature adipocytes 127 were subsequently isolated as previously reported (Jacobsen et al. 2019). UNIK minipigs 128 BMI profiles are detailed in Table S1. 129 RNA-seq and small RNA-seq expression data generated in the framework of these two 130 experimental systems have been described previously (Cardoso et al. 2017; Jacobsen et 131 al. 2019; Mármol-Sánchez et al. 2020). Briefly, sequencing reads from the RNA-seq and 132 small RNA-seq datasets were trimmed with the Cutadapt software (Martin, 2011). RNA- 133 seq reads were then mapped with the HISAT2 aligner (Kim et al. 2019) using default 134 parameters. The Bowtie v.1.2.1.1 software (Langmead et al. 2009) was used to align small 135 RNA-seq reads by considering small sequence reads specifications (bowtie -n 0 -m 20 -k 136 1 --best). The Sscrofa11.1 porcine reference assembly (Warr et al. 2020) was used for 137 mapping. 138 139 Exon/Intron quantification 140 We generated exonic and intronic-specific annotations spanning all available genes by 141 using the Sscrofa.11.1 v.103 gene annotation (Ensembl repositories: 142 http://ftp.ensembl.org/pub/release-103/gtf/sus_scrofa/). Overlapping intronic/exonic 143 regions (10 bp at both ends of introns), as well as singleton positions were removed 144 (Lawrence et al. 2013). We then used the featureCounts tool (Liao et al. 2014) to 145 independently quantify exon and intron expression levels for each mRNA encoding gene, 146 as well as expression levels of mature miRNAs. 147 6 148 Differential expression analyses 149 Canonical differential expression analyses were carried out with the edgeR tool 150 (Robinson et al. 2010) by considering only the exonic counts of mRNAs and the mature 151 miRNA expression levels measured in the two experimental systems under study. Only 152 genes showing average expression values above 1 count-per-million in at least 50% of 153 animals within each treatment group (AL-T0 and AL-T2 for GM tissue and lean and obese 154 for adipocyte isolates) were retained for further analyses. Expression filtered raw counts 155 for exonic mRNA and miRNA reads were normalized with the trimmed mean of M- 156 values normalization (TMM) method (Robinson & Oshlack 2010) and the statistical 157 significance of mean expression differences was tested with a quasi-likelihood F-test 158 (Robinson et al. 2010). Multiple hypothesis testing correction was implemented with the 159 false discovery rate method (Benjamini & Hochberg 1995). Messenger RNAs were 160 considered to be significantly differentially expressed when the absolute value of the fold- 161 change (FC) was higher than 2 (|FC| > 2) and the q-value < 0.05. Fasting Duroc gilts (AL- 162 T0) as well as UNIK obese minipigs were classified as baseline controls in differential 163 expression analyses, i.e. any given upregulation represents and overexpression of the 164 corresponding gene in fed (AL-T2) Duroc gilts or lean UNIK minipigs with respect to 165 their fasting (AL-T0) and obese counterparts, respectively. 166 Since changes in the expression of miRNAs are often subtler than those of mRNAs, the 167 thresholds to consider a miRNA as significantly differentially expressed were set to |FC| 168 > 1.5 and q-value < 0.05 (Guo et al. 2015) for fasted (AL-T0) vs fed (AL-T2) Duroc gilts. 169 Given the relatively low statistical significance of miRNA expression changes observed 170 in obese vs lean UNIK minipigs, we considered, as potential miRNA regulators, those 171 that were significantly (FC < -1.5; q-value < 0.05) and suggestively (FC > 1.5 and P- 172 value < 0.01) upregulated in lean minipigs. 7 173 174 Exon/intron split analysis (EISA). 175 We applied EISA to differentiate transcriptional from post-transcriptional changes in 176 mRNA expression in our two experimental systems (muscle and fat). Normalization was 177 performed independently for exon and intron counts as described by Gaidatzis et al. 178 (2015). Exonic and intronic gene abundances were subsequently log2 transformed, adding 179 a pseudo-count of 1 and averaged within each considered treatment group. 180 Only genes for which both exonic and intronic read counts were successfully quantified 181 were further considered. Observed differences in each ith gene were expressed as the 182 increment of exonic/intronic counts in fed (AL-T2) and obese animals with respect to 183 fasting (AL-T0) and lean animals, respectively. In this way, the increment of intronic and 184 exonic counts was calculated considering ΔInt = Int2i – Int1i and ΔEx = Ex2i – Ex1i, 185 respectively. The magnitudes of the transcriptional (Tc) and post-transcriptional (PTc) 186 changes in mRNA expression were then calculated. The Tc contribution to the observed 187 counts is explained by ΔInt (Pillman et al. 2019), while PTc can be deduced from ΔEx – 188 ΔInt. In this way, the significance of Tc scores was assessed as in canonical differential 189 expression analyses but using the intronic fraction as input. Moreover, we took advantage 190 of the generalized linear model framework from edgeR tool to assess the significance of 191 PTc scores by introducing an interaction term between the fraction type (exonic or 192 intronic) and condition type (fasting AL-T0 vs fed AL-T2 or obese vs lean). Both Tc and 193 PTc components were z-scored to make ΔEx and ΔInt estimates comparable. Tc and PTc 194 scores were considered significant when |FC| > 2 and q-value < 0.05. Multiple hypothesis 195 testing correction was implemented by using the false discovery rate approach (Benjamini 196 & Hochberg 1995). 8 197 In order to obtain a prioritized list of genes showing high post-transcriptional regulation 198 signals, the top 5% of expressed genes with the most negative PTc scores were retrieved 199 (irrespective of their statistical significance after multiple testing correction). From these, 200 we only focused on genes showing strongly reduced ΔEx values of at least 2-fold in both 201 experimental systems (i.e. ΔEx < -1 in the log2 scale). All implemented analyses have 202 been summarized in Fig. S1. A ready-to-use modular pipeline for running EISA is 203 publicly available at https://github.com/emarmolsanchez/EISAcompR. 204 205 miRNA target prediction 206 Putative interactions between the seed of expressed miRNAs (2nd-8th 5’ nts) and the 3’- 207 UTRs of expressed mRNAs were predicted on the basis of sequence identity using the 208 Sscrofa11.1 reference assembly and the seedVicious v1.1 tool (Marco 2018). The 209 annotated 3’-UTRs longer than 30 nts from porcine mRNAs were retrieved from 210 http://www.ensembl.org/biomart, while mature porcine miRNA sequences were obtained 211 from miRBase (Kozomara et al. 2019). Redundant miRNA seeds were removed, and 212 8mer, 7mer-m8 and 7mer-A1 miRNA-mRNA interactions were taken into account 213 (Bartel 2018). 214 Based on Grimson et al. (2007), in silico-predicted miRNA-mRNA interactions matching 215 any of the following criteria were removed: (i) Binding sites located in 3’-UTRs at less 216 than 15 nts close to the end of the open reading frame (and the stop codon) or less than 217 15 nts close to the end of the 3’-UTR and the beginning of the terminal poly(A) tail (E 218 criterion), (ii) binding sites located in the middle of the 3’-UTR in a range comprising 45- 219 55% of the central region of the non-coding sequence (M criterion), and (iii) binding sites 220 that lack AU-rich sequences in their immediate upstream and downstream flanking 9 221 regions comprising 30 nts each (AU criterion). A schematic representation of such criteria 222 is available at Fig. S2. 223 Covariation patterns between miRNAs and their predicted mRNA targets were assessed 224 by computing Spearman’s correlation coefficients (ρ) with the TMM normalized and log2 225 transformed expression profiles of the exonic fraction of mRNA and miRNA genes. 226 227 miRNA target enrichment analyses 228 We predicted which downregulated mRNA genes, from those with highly negative post- 229 transcriptional signals, are putatively targeted by at least one of the significantly 230 upregulated miRNAs. Subsequently, we investigated whether the sets of mRNA genes 231 identified in this way were enriched in being targets of upregulated miRNAs, compared 232 to the whole set of expressed mRNAs genes with available 3’-UTRs (control 233 background). Enrichment analyses were carried out using the Fisher’s exact test in R. 234 Significance level was set at a nominal P-value < 0.05. 235 We also tested whether downregulated mRNA genes with highly negative post- 236 transcriptional signals were significantly enriched to be targets of at least one of the top 237 5% most highly expressed miRNA genes (considering their overall average expression 238 and excluding significantly upregulated miRNAs), as well as of significantly 239 downregulated miRNAs. 240 As an additional randomized control test for enrichment analyses, we generated 100 241 random sets of 10 expressed mature miRNA genes without seed redundancy. In this way, 242 we predicted which downregulated mRNA genes, from those with highly negative post- 243 transcriptional signals, were putatively targeted by at least one miRNA from the random 244 sets defined. The distribution of odds ratios obtained after enrichment tests over each 10 245 random set of miRNAs (N = 100) were then compared with the observed odds ratios 246 obtained with enrichment analyses using the set of significantly upregulated miRNAs. 247 The P-value for the significance of the deviation of observed odds ratios against the 248 bootstrapped odds ratios distribution was defined as: 𝑟+1 249 𝑃 − 𝑣𝑎𝑙𝑢𝑒 = 1 − , where r is the number of permuted odds ratios with values 𝑘+1 250 equal or higher than the observed odds ratio, and k is the number of permutations (N = 251 100). 252 253 Gene covariation network and covariation enrichment score 254 We used edgeR to identify mRNA genes in the AL-T0 vs AL-T2 and obese vs lean 255 comparisons showing q-value < 0.05, after multiple testing correction. Then, the 256 normalized exonic and intronic estimates in the log2 scale obtained from EISA were 257 independently used to compute Spearman’s correlation coefficients (ρ). Significant 258 correlations were identified with the Partial Correlation with Information Theory (PCIT) 259 algorithm (Reverter & Chan 2008) implemented in the pcit R package (Watson-Haigh et 260 al. 2010). In this way, we calculated a covariation enrichment score (CES), as reported 261 by Tarbier et al. (2020), to assess the potential contribution of miRNAs to the observed 262 differences in covariation. This test compares the number of overall significant pairwise 263 correlations with those obtained when only considering the set of downregulated mRNA 264 genes with highly negative post-transcriptional signals and putatively targeted by 265 significantly upregulated miRNAs. Further details about the algorithm used to calculate 266 the CES values and control tests can be found in Supplementary Methods. Significant 267 differences among the set of exonic, intronic and control CES values were tested with a 268 Mann-Whitney U non-parametric test (Mann & Whitney 1947). 269 11 270 Verification of RNA-seq expression profiles by qPCR 271 By using qPCR, we have verified that the expression profiles of selected mRNAs and 272 miRNAs highlighted in the adipose tissue experiment were in accordance with those 273 obtained with RNA-seq and small RNA-seq data. Since this experiment is just a control 274 of the quality of expression estimates, further details are reported in Supplementary 275 Methods. Primers for mRNA and miRNA qPCR expression profiling are available at 276 Table S2. Raw Cq values for each assay are available at Table S3. 277 278 279 Results 280 The analysis of post-transcriptional regulation in muscle samples from fasting and 281 fed Duroc gilts 282 Differential expression and EISA 283 Total RNA and small RNA were independently sequenced in GM muscle samples from 284 fasting (AL-T0) and fed (AL-T2) Duroc gilts. About 45.2 million reads (93%) per sample 285 from protein coding and non-coding genes were successfully mapped against the 286 Sscrofa.11.1 assembly when analyzing the RNA-seq data. Besides, around 2.2 million 287 reads per sample (42%) obtained from the small RNA-seq experiment were successfully 288 mapped to 370 annotated porcine miRNA genes. 289 A total of 30,322 (based on exonic reads) and 22,769 (based on intronic reads) genes were 290 successfully quantified after splitting the reference genome assembly between exonic and 291 intronic features. Exonic counts were ~23 fold more abundant than those corresponding 292 to intronic regions. 293 By using edgeR, we detected 454 mRNA genes with q-value < 0.05 (Table S4a). Among 294 these, only genes with |FC| > 2 were retained, resulting in 52 upregulated and 80 12 295 downregulated genes (Table S4a in bold and Fig. S3a). The analysis of small RNA-seq 296 data with edgeR revealed 16 miRNAs significantly differentially expressed, of which 8 297 were upregulated in AL-T2 gilts. These 8 miRNAs, which represented 6 unique miRNA 298 seeds (ssc-miR-148a-3p, ssc-miR-7-5p, ssc-miR-30-3p, ssc-miR-151-3p, ssc-miR-374a- 299 3p and ssc-miR-421-5p; Table S5 in bold), were selected as potential post-transcriptional 300 regulators of mRNA genes. 301 On the other hand, EISA highlighted 26 mRNA genes displaying the top 5% negative 302 PTc scores with at least 2-fold ΔEx reduction (Table 1 and Fig. S3b). Eighteen out of 303 these 26 genes (69.23%) appeared as significantly downregulated (FC < -2; q-value < 304 0.05, Table 1 and Table S4b in bold) according to canonical differential expression 305 analyses. The whole list of expressed mRNA genes after canonical differential expression 306 analyses is available at Table S4c. 307 Also, we detected 133 mRNA genes with significant PTc scores (|FC| > 2; q-value < 0.05, 308 Table S6a), of which three experienced at least a 2-fold reduction of their ΔEx fraction 309 (Table S6a in bold). Two out from these three mRNA genes ranked within those with the 310 top 5% negative PTc scores (Table 1 and Table S6a). Among this set of 133 genes, only 311 seven (5.26%) were also significantly differentially expressed (Table S6a). Moreover, 312 with EISA we detected 344 genes displaying significant Tc scores (|FC| > 2; q-value < 313 0.05, Table S6b) and among these, 71 (20.63%) were also significantly differentially 314 expressed (Tables S6b). Besides, 91 out of these 344 genes (26.45%) also showed 315 significant PTc scores (Table S6b in bold), but none of them were among the mRNA 316 genes displaying the top 5% negative PTc scores with at least 2-fold ΔEx reduction (N = 317 26, Table 1). The whole lists of expressed mRNA genes after EISA and their PTc and Tc 318 scores are available at Table S6c and S6d, respectively. 319 13 320 Context-based pruning of predicted miRNA-mRNA interactions removes unreliable 321 target events 322 Before making in silico predictions of miRNA-mRNA interactions, we investigated their 323 reliability. To do so, we first evaluated the enrichment in the number of genes with 324 binding sites for at least one of the 6 non-redundant miRNAs upregulated in the GM 325 muscle of AL-T2 gilts (ssc-miR-148a-3p, ssc-miR-7-5p, ssc-miR-30-3p, ssc-miR-151-3p, 326 ssc-miR-374a-3p and ssc-miR-421-5p) over a background of all expressed genes with no 327 context-based removal of predicted binding sites (see Methods). Introducing additional 328 context-based filtering criteria to remove unreliable binding site predictions resulted in 329 an overall increased enrichment of target genes within the list of the top 1% (N = 13 330 genes, Fig. S4a) and 5% (N = 26 genes, Fig. S4b) genes with negative PTc scores and 331 displaying at least 2-fold ΔEx reduction. This enrichment was more evident when using 332 the AU criterion, as shown in Fig. S4a. 333 334 Several genes with relevant post-transcriptional signals detected with EISA are predicted 335 to be targets of upregulated miRNAs 336 Target prediction and context-based pruning of miRNA-mRNA interactions for mRNA 337 genes displaying the top 5% negative PTc scores and at least 2-fold reduction in the ΔEx 338 exonic fraction (N = 26, Table 1, Fig. 1a) made possible to detect 11 8mer, 21 7mer-m8 339 and 22 7mer-A1 miRNA binding sites for the six non-redundant seeds of miRNAs 340 significantly upregulated in AL-T2 gilts (Table S5 in bold) in 21 out of the 26 analyzed 341 mRNAs (80.77%, Table S7). Moreover, 14 out of these 21 genes (66.67%) were also 342 significantly differentially expressed (Table 1 and Table S4b in bold). 343 This set of 21 mRNA genes with putative post-transcriptional repression mediated by 344 miRNAs showed a significant enrichment in 8mer, 7mer-m8 and 7mer-A1 sites for the 6 14 345 miRNAs significantly upregulated in AL-T2 gilts, especially when combining the three 346 types of miRNA binding sites considered (Fig. 1b). The miRNAs with the highest number 347 of significant miRNA-mRNA interactions were ssc-miR-30a-3p and ssc-miR-421-5p, 348 which showed nine and eight significant interactions, followed by ssc-miR-148-3p with 349 four significant interactions (Table S7). 350 We also evaluated the enrichment of the number of mRNA genes within the list of the 351 top 5% negative PTc scores and at least 2-fold ΔEx reduction (N = 26, Table 1) to be 352 targets of at least one of the following: (i) miRNAs downregulated in AL-T2 fed gilts 353 (Table S5), (ii) top 5% most expressed miRNAs, excluding those significantly 354 upregulated (ssc-miR-1, ssc-miR-133a-3p, ssc-miR-26a, ssc-miR-10b, ssc-miR-378, ssc- 355 miR-99a-5p, ssc-miR-27b-3p, ssc-miR-30d, ssc-miR-486 and ssc-let-7f-5p), and (iii) 356 random sets (N = 100) of 10 expressed miRNAs, as a control test. We did not detect a 357 significant enrichment in any of these three additional analyses (Fig. 1b). 358 The mRNA with the highest negative and most significant PTc score was the Dickkopf 359 WNT Signaling Pathway Inhibitor 2 (DKK2), being a strong candidate to be repressed by 360 miRNAs (Table 1). Indeed, DKK2 was the only gene harboring two 8mer binding sites 361 (Table S7). Interestingly, this locus was almost significantly differentially expressed 362 (Table 1 and Table S4b). The discordance between EISA and canonical differential 363 expression results can be fully appreciated when comparing Fig. 1a (genes with high post- 364 transcriptional repression after EISA) and Fig. 1c (canonical differential expression 365 analyses), where only 18 out of the 26 mRNA genes detected with EISA appeared as 366 significantly downregulated in the edgeR-based differential expression analyses (Table 367 1). Although several of the mRNA genes shown in Table 1 were highly downregulated 368 (Table S4b), the majority were mildly to slightly downregulated or not significantly 369 differentially expressed. 15 370 371 Genes showing post-transcriptional regulatory signals predominantly covary at the 372 exonic level 373 To further elucidate whether the set of 21 mRNA genes with putative post-transcriptional 374 repression mediated by upregulated miRNAs showed covarying expression profiles, we 375 evaluated the number of significant co-expressed pairs among them and among the whole 376 set of mRNA genes with q-value < 0.05, including those from the set of 21 mRNA genes 377 with q-value > 0.05 after canonical differential expression analyses (Table 4b). 378 By calculating CES values for both exonic and intronic fractions (see Methods) of the 21 379 genes putatively targeted by the 6 significantly upregulated miRNAs in fed gilts, our 380 analyses revealed that 19 out of these 21 genes showed increased covariation in their 381 exonic fraction when compared to their intronic fraction (Table S8, Fig. 1d), and DKK2 382 was again the gene with the strongest exonic covariation fold change compared to its 383 intronic covariation (Table S8). As expected, control random sets of genes (N = 1,000) 384 displayed CES ≈ 1, indicative of no covariation (Fig. 1d). The observed CES distributions 385 of exonic and intronic sets were significantly different (P-value = 3.663E-06) after 386 running non-parametric tests (Fig. 1d), thus supporting that the majority of these 21 genes 387 might be indeed co-regulated at the post-transcriptional level by upregulated miRNAs. 388 389 Studying post-transcriptional signals in adipose tissue using the UNIK minipig 390 population 391 After pre-processing and filtering of sequenced reads from adipocyte samples, we were 392 able to retrieve ~98.1 and ~0.87 million mRNA and small RNA reads per sample, and 393 ~96.5% and ~73.4% of these reads mapped to annotated porcine mRNAs and mature 394 miRNAs, respectively. Canonical differential expression analyses revealed a total of 299 16 395 genes with q-value < 0.05, of which 52 and 95 were significantly downregulated and 396 upregulated (FC > |2|; q-value < 0.05), respectively (Table S9a in bold). Only one 397 miRNA (ssc-miR-92b-3p) was significantly upregulated in lean minipigs (Table S10), 398 while six additional miRNAs showed suggestive differential expression (P-value < 0.01), 399 of which four were upregulated (ssc-miR-148a-3p, ssc-miR-204, ssc-miR-92a and ssc- 400 miR-214-3p; Table S10 in bold). 401 After running EISA, only the sestrin 3 (SESN3) gene showed a significant PTc score, 402 having the second highest negative PTc score (Table 2 and Table S11a). Moreover, 403 SESN3 was also detected as the most significantly downregulated gene by edgeR (Table 404 S9a). A total of 44 downregulated mRNAs in lean minipigs displayed the top 5% PTc 405 scores with reduced ΔEx of at least 2-fold (Table 2 and Fig. 2a). Among them, only 12 406 (27.27%) appeared as significantly downregulated (FC < -2; q-value < 0.05, Table 2 and 407 Table S9b in bold). The whole list of expressed genes after differential expression 408 analyses are available at Table S9c. 409 Besides, 25 of these 44 (58.14%, Table S9b) mRNAs were classified as putative targets 410 of the set of miRNAs upregulated in lean minipigs (N = 4, ssc-miR-92b-3p, ssc-miR- 411 148a-3p, ssc-miR-204 and ssc-miR-214-3p; Table S10 in bold). Target prediction and 412 context-based pruning of miRNA-mRNA interactions for these 25 genes made possible 413 to detect eight 8mer, 21 7mer-m8 and 24 7mer-A1 miRNA binding sites (Table S12) for 414 upregulated miRNAs (N = 4) in lean UNIK minipigs (Table S10). Again, the SESN3 415 gene showed the highest number of predicted putative miRNA target sites in its 3´-UTR 416 (Table S12). 417 Enrichment analyses for these 25 mRNAs (Table 2) showed no significant results for 418 8mer, 7mer-m8 and 7mer-A1 miRNA binding sites (Fig. 2b). Moreover, only seven of 419 them (26.92%) appeared as significantly downregulated (FC < -2; q-value < 0.05) in the 17 420 canonical differential expression analyses with edgeR (Table S9b and Table 2, Fig. 2c). 421 The exonic fractions of 18 out of these 25 mRNA genes showed significantly increased 422 covariation (P-value = 2.703E-02) compared to the covariation observed for the intronic 423 fractions (Fig. 2d and Table S13). 424 Regarding Tc scores, a total of 195 genes showed significant transcriptional signals (|FC| 425 > 2; q-value < 0.05, Table S11b), and 48 of them were also significantly differentially 426 expressed (24.61%, Tables S9a and S11b). Moreover, three of them (ARHGAP27, CDH1 427 and LEP) were found among those with the top 5% post-transcriptional signals (Table 2 428 and Table S11b in bold). The whole lists of expressed genes after EISA and their PTc 429 and Tc scores are available at Table S11c and S11d, respectively. 430 Results obtained for qPCR verification analyses are described in Fig. S5. 431 432 433 Discussion 434 Contribution of the Tc and PTc components of gene regulation to energy 435 homeostasis in porcine muscle and adipose tissues 436 After running EISA on both muscle and adipose tissue datasets, we observed that the 437 number of genes with significant transcriptional signals (Tc) was much higher than that 438 of loci with significant post-transcriptional signals (PTc). Such difference evidences that 439 gene expression changes induced by feeding or adiposity might be mostly driven by 440 transcriptional rather than post-transcriptional modulators. It is worth noting, however, 441 that relatively few mRNAs showed post-transcriptional signals alone and were thus 442 mixed with transcriptional signals either in concordant or in opposite directions. 443 For prioritizing putative post-transcriptionally repressed mRNA genes by miRNAs, we 444 focused on those with the strongest observed downregulation based on their ΔEx values 18 445 (at least 2-fold reduction) and PTc signal (top 5% negative scores). Hence, we did not 446 consider the significance of PTc scores as a relevant criterion, as these will appear as 447 significant only when the post-transcriptional response is strong and not confounded by 448 a cooperative repressive transcriptional signal. Alternative thresholds other than the top 449 5% negative PTc scores or 2-fold for ΔEx fraction could be applied, depending on the 450 strength of the post-transcriptional signal detected. 451 Besides, it is worth noting that our RNA-seq data was generated following manufacturer’s 452 instructions for TruSeq stranded libraries. This commonly used protocol selects for 453 poly(A) mRNAs, and the majority of non-poly(A) intronic lariats generated after splicing 454 will be lost, thus producing an artificial decrease in the quantified intronic fraction 455 (intronic reads). Although the intronic yield is decreased in poly(A) RNA compared to 456 total RNA protocols (Ameur et al. 2011), the recoverable intronic fraction is still highly 457 correlated with nascent transcription (Gaidatzis et al. 2015; La Manno et al. 2018), 458 probably derived from the presence of poly(A) introns or unspliced mRNAs being 459 sequenced, albeit at low abundance. Thus, the remaining intronic reads might conform a 460 limited and indirect yet representative proxy of the transcriptional activity. Alternative 461 methods for measuring such transcriptional activity by directly monitoring transcription 462 across the genome have also been applied (Patel et al. 2020), and might be preferred over 463 the use of intronic fractions. However, the implementation of such alternatives is limited 464 so far, and the relatively simple and cost-effective usage of intronic reads present in 465 already available RNA-seq data justifies the use of the EISA approach. 466 467 Canonical differential expression analyses and EISA highlight different sets of genes 468 Few genes with significant Tc and PTc components were also classified by edgeR as 469 significantly differentially expressed. Such a discrepancy between EISA and differential 19 470 expression analyses is in agreement with the subtle regulation elicited by miRNAs, which 471 is dependent on the expression level of miRNAs and the number of binding sites within 472 a given mRNA 3’-UTR (Bartel 2018). In this way, EISA might serve as a good approach 473 to identify both strong and subtle post-transcriptional effects mediated by miRNAs that 474 canonical differential expression approaches might not be able to capture. 475 Importantly, these discrepancies were reduced when we focused on genes with top 5% 476 negative PTc scores and at least 2-fold reduction in their ΔEx values: as much as 69.23% 477 (skeletal muscle) and 27.27% (adipose tissue) of such genes were also detected as 478 differentially expressed. This increase in concordance was more pronounced in the 479 skeletal muscle (fasted vs fed gilts) experimental system. This might be due to the overall 480 stronger upregulation of miRNAs observed for this dataset when compared with that 481 generated in the adipose tissue experiment, which can be explained by intrinsic genomic 482 differences among pig breeds, the tissues analyzed and/or the metabolic challenge 483 undertaken. 484 485 Predicting the contribution of miRNAs to the post-transcriptional regulatory 486 response in porcine muscle and adipose tissues 487 Since the efficacy of miRNA-based repression of mRNA expression depends on the 488 context of the miRNA binding site within the 3’-UTR (Grimson et al. 2007), we have 489 assessed the usefulness of introducing context-based filtering criteria for removing 490 unreliable in silico-predicted binding sites for miRNAs. Using enrichment analyses, we 491 were able to link the expression of the set of mRNAs with downregulated exonic fraction 492 to the expression of upregulated miRNAs predicted to target them. 493 In the skeletal muscle system, prediction of miRNA binding sites in mRNA genes 494 displaying the top 5% negative post-transcriptional signals and at least 2-fold reduction 20 495 in their exonic fraction (N = 26) revealed that the majority of them (80.77%) harbored at 496 least one binding site for the corresponding set of significantly upregulated miRNAs (N 497 = 6). In contrast, such pattern was much less evident in the adipose tissue, which could 498 be explained by the fact that the majority of upregulated miRNAs in lean minipigs did 499 not reach significance. 500 Although we verified by qPCR the RNA-seq expression levels of selected downregulated 501 mRNAs and upregulated miRNAs in the UNIK adipose tissue experiment, further 502 experimental validation of the reported mRNA-miRNA interactions is needed. In silico 503 predictions of miRNA binding sites, as well as EISA and covariation analyses, might be 504 helpful to identify and prioritize miRNA-mRNA pairs to be experimentally validated with 505 co-transfection gene reporter assays. In this way, the yet scarce collection of validated 506 mRNA-miRNA interactions in domestic species could be expanded and improved. 507 508 Covariation patterns in the expression of downregulated mRNAs predicted to be 509 targeted by upregulated microRNAs 510 We further hypothesized that mRNA genes showing relevant post-transcriptional 511 downregulation might be repressed by the same set of significantly upregulated miRNAs, 512 which could induce shared covariation in the expression profiles of such mRNAs at the 513 exonic level. In contrast, their intronic fraction would be mainly unaffected as introns 514 would have been excised prior to any given miRNA-driven downregulation. Therefore, 515 an increased gene covariation in downregulated mRNAs with high post-transcriptional 516 signals might be detectable at the exon but not at the intron level. Indeed, our results 517 revealed an increased covariation in downregulated mRNAs with high post- 518 transcriptional signals at their exonic fraction compared with covariation patterns of their 21 519 intronic fraction, suggesting that their expression might be repressed by a common set of 520 upregulated miRNAs. 521 522 Genes displaying the strongest post-transcriptional signals in porcine skeletal 523 muscle and adipose tissue are involved in glucose and lipid metabolism 524 The mRNA genes showing the strongest post-transcriptional downregulation in fasted vs 525 fed gilts displayed a variety of relevant biological functions. The DKK2 gene showed the 526 most negative significant PTc score. This gene also displayed the strongest covariation 527 difference in its exonic fraction compared with the intronic one. This consistent post- 528 transcriptional regulatory effect might be mediated by ssc-miR-421-5p and ssc-miR-30a- 529 3p, two highly significantly upregulated miRNAs. The DKK2 protein is a member of the 530 dickkopf family, which inhibits the Wnt signaling pathway through its interaction with 531 the LDL-receptor related protein 6 (LRP6). Its repression has been associated with 532 reduced blood-glucose levels and improved glucose uptake (Li et al. 2012), as well as 533 with improved adipogenesis (Yang & Shi 2021) and the inhibition of aerobic glycolysis 534 (Deng et al. 2019). These results are consistent with the increased glucose usage and 535 triggered adipogenesis in muscle tissue after nutrient supply. Other relevant post- 536 transcriptionally downregulated mRNAs detected with EISA were: pyruvate 537 dehydrogenase kinase 4 (PDK4), interleukin 18 (IL18), nuclear receptor subfamily 4 538 group A member 3 (NR4A3), acetylcholine receptor subunit α (CHRNA1), PBX 539 homeobox 1 (PBX1), Tet methylcytosine dioxygenase 2 (TET2), BTB domain and CNC 540 homolog (BACH2), all of which are involved in the regulation of energy homeostasis 541 (Lindegaard et al. 2013; Zhang et al. 2014; Tamahara et al. 2017; Wu et al. 2018; Xu et 542 al. 2019) and lipid metabolism in muscle cells (Monteiro et al. 2011; Pearen et al. 2013) 543 in response to nutrient uptake. 22 544 On the other hand, several genes with high post-transcriptional signals were not predicted 545 to be targets of upregulated miRNAs: the circadian associated repressor of transcription 546 (CIART) and period 1 (PER1), oxysterol binding protein like 6 (OSBPL6) and nuclear 547 receptor subfamily 4 group A member 1 (NR4A1). Interestingly, all of them were 548 significantly downregulated (Cardoso et al. 2017; Mármol-Sánchez et al. 2020). An 549 explanation to this might be that although miRNAs are key post-transcriptional 550 regulators, other alternative post-transcriptional effectors, such as long non-coding 551 RNAs, circular RNAs or RNA binding proteins might be at play. Besides, indirect 552 repression via upregulated miRNAs acting over regulators of these genes, such as 553 transcription factors, could be also a major influence on their observed repression (Patel 554 et al. 2020). 555 The use of EISA on expression data from adipocytes isolated from obese vs lean UNIK 556 minipigs revealed several mRNA genes with high post-transcriptional repression, which 557 are also involved in the regulation of lipid metabolism and energy homeostasis. The gene 558 showing the highest post-transcriptional signal was the estrogen related receptor γ 559 (ESRRG), which modulates oxidative metabolism and mitochondrial function in adipose 560 tissue and inhibits adipocyte differentiation when repressed (Kubo et al. 2009). Another 561 relevant locus identified with EISA was SESN3, an activator of the mTORC2 and 562 PI3K/AKT signaling pathway that promotes lipolysis when inhibited (Tao et al. 2015). 563 This gene showed the most significant downregulation in lean minipigs, and gathered 564 multiple putative binding sites for all the four upregulated miRNAs under study. 565 Other genes showing significant post-transcriptional downregulation were: sterile α motif 566 domain containing 4A (SAMD4A), prostaglandin F2- receptor protein (PTGFR), serine 567 protease 23 (PRSS23), ring finger protein 157 (RNF157), oxysterol binding protein like 568 10 (OSBPL10), glycosylphosphatidylinositol phospholipase 1 (GPLD1), RAP1 GTPase 23 569 activating protein (RAP1GAP) and leptin (LEP), all of which are tightly linked to 570 adipocyte differentiation (Perttilä et al. 2009; Martínez et al. 2013; Chen et al. 2014; 571 Kosacka et al. 2018) or energy homeostasis (Ussar et al. 2012; Wang et al. 2018; 572 Izquierdo et al. 2019; Kuo et al. 2020). Despite the overall weak influence of putative 573 miRNA-driven downregulation on mRNAs expressed in adipocytes, we were able to 574 identify a set of genes with high post-transcriptional signals indicative of putative 575 miRNA-derived repression and tightly related to adipose tissue metabolism regulation. 576 However, non-miRNA transcriptional and post-transcriptional modulators might also 577 contribute to such repression. 578 579 580 Conclusions 581 EISA applied to study gene regulation in porcine skeletal muscle and adipose tissues 582 showed that more genes were subjected to transcriptional rather than post-transcriptional 583 regulation, suggesting that changes in mRNA expression are mostly driven by factors 584 acting at the transcriptional level. More importantly, the concordance between the sets of 585 significantly differentially expressed genes and those with significant Tc or PTc scores 586 was quite limited, but improved (mostly in the skeletal muscle experiment) when we 587 prioritized the downregulated genes with the top 5% negative post-transcriptional signals. 588 Nevertheless, many of the genes with relevant PTc signals were not among the top 589 significantly downregulated loci, thus demonstrating the usefulness of complementing 590 canonical differential expression analyses with the EISA approach. In the skeletal muscle, 591 we detected several mRNAs predicted to be co-regulated by a common set of miRNAs. 592 In contrast, in the adipose tissue such relationship was more subtle, suggesting that the 593 contribution of miRNAs to mRNA repression might be affected by tissue type, breed 24 594 and/or intrinsic experimental factors. Finally, EISA made possible to identify several 595 genes related with carbohydrate and lipid metabolism, which may play relevant roles in 596 the energy homeostasis of the skeletal muscle and adipose tissues. By differentiating the 597 transcriptional from the post-transcriptional changes in mRNA expression, EISA 598 provides a valuable view, complementary to canonical differential expression analyses, 599 about the miRNA-driven regulation of gene expression. 600 601 602 Ethics approval 603 Animal care and management procedures for Duroc gilts followed national guidelines for 604 the Good Experimental Practices and were approved by the Ethical Committee of the 605 Institut de Recerca i Tecnologia Agroalimentàries (IRTA). Animal care and management 606 procedures for UNIK minipigs were carried out according to the Danish “Animal 607 Maintenance Act” (Act 432 dated 9 June 2004). 608 609 Availability of data 610 The RNA-seq and small RNA-seq datasets from skeletal muscle tissue used in the current 611 study are available at the Sequence Read Archive (SRA) database with BioProject codes 612 PRJNA386796 and PRJNA595998, respectively. For the adipose tissue samples, RNA- 613 seq and small RNA-seq datasets are available at PRJNA563583 and PRJNA759240. 614 615 Competing interests 616 The authors declare that they have no competing interests 617 618 Funding 25 619 The present research work was funded by grants AGL2013-48742-C2-1-R and 620 AGL2013-48742-C2-2-R awarded by the Spanish Ministry of Economy and 621 Competitivity. E. Mármol-Sánchez was funded with a PhD fellowship FPU15/01733 622 awarded by the Spanish Ministry of Education and Culture (MECD). YRC is recipient of 623 a Ramon y Cajal fellowship (RYC2019-027244-I) from the Spanish Ministry of Science 624 and Innovation. 625 626 Authors’ contributions 627 The authors’ responsibilities were as follows: M.A. and R.Q. generated the total RNA 628 and small RNA sequencing data corresponding to the GM muscle. S.C., M.J.J., C.B.J. 629 and M.F produced the total RNA and small RNA sequencing corresponding to the adipose 630 tissue. M.A., R.Q., S.C., M.J.J., C.B.J., M.F., T.F.C. and E.M.-S. conducted the research. 631 S.C. performed qPCR analyses. E.M.-S. analyzed the data. L.M.Z. contributed to 632 bioinformatic analyses. Y.R.-C. contributed to critical assessment. M.A. and 633 R.Q. secured funding for the study. E.M.-S., S.C. and M.A. drafted the manuscript. All 634 authors contributed to manuscript corrections, read and approved the final manuscript. 635 636 Acknowledgements 637 The authors would like to thank the Department of Veterinary Animal Sciences in the 638 Faculty of Health and Medical Sciences of the University of Copenhagen for providing 639 sequencing data and their facilities and resources for qPCR experiments. We also 640 acknowledge Selección Batallé S.A. for providing animal material and the support of the 641 Spanish Ministry of Economy and Competitivity for the Center of Excellence Severo 642 Ochoa 2020–2023 (CEX2019-000902-S) grant awarded to the Centre for Research in 26 643 Agricultural Genomics (CRAG, Bellaterra, Spain). Thanks also to the CERCA 644 Programme of the Generalitat de Catalunya for their support. 645 646 647 References 648 Ali A., Murani E., Hadlich F., Liu X., Wimmers K. & Ponsuksili S. (2021) In utero fetal 649 weight in pigs is regulated by microRNAs and their target genes. Genes 12, 1264. 650 Ameur A., Zaghlool A., Halvardson J., Wetterbom A., Gyllensten U., Cavelier L. et al. 651 (2011) Total RNA sequencing reveals nascent transcription and widespread co- 652 transcriptional splicing in the human brain. Nature Structural & Molecular Biology 653 18, 1435–40. 654 Ballester M., Amills M., González-Rodríguez O., Cardoso T.F., Pascual M., González- 655 Prendes R. et al. (2018) Role of AMPK signaling pathway during compensatory 656 growth in pigs. BMC Genomics 9, 682. 657 Bartel D.P. (2018) Metazoan microRNAs. Cell 173, 20–51. 658 Benítez R., Trakooljul N., Núñez Y., Isabel B., Murani E., De Mercado E. et al. (2019) 659 Breed, diet, and interaction effects on adipose tissue transcriptome in Iberian and 660 duroc pigs fed different energy sources, Genes 10, 589. 661 Benjamini Y. & Hochberg Y. (1995) Controlling the false discovery rate: a practical and 662 powerful approach to multiple testing. Journal of the Royal Statistical Society 663 Series B (Methodological) 57, 289–300. 664 Cardoso T.F., Quintanilla R., Tibau J., Gil M., Mármol-Sánchez E., González-Rodríguez 665 O. et al. (2017) Nutrient supply affects the mRNA expression profile of the porcine 666 skeletal muscle. BMC Genomics 18, 603. 667 Chen Z., Holland W., Shelton J.M., Ali A., Zhan X., Won S, et al. (2014) Mutation of 27 668 mouse Samd4 causes leanness, myopathy, uncoupled mitochondrial respiration, 669 and dysregulated mTORC1 signaling. Proceedings of the National Academy of 670 Sciences (USA) 111, 7367–7372. 671 Cursons J., Pillman K.A., Scheer K.G., Gregory P.A., Foroutan M., Hediyeh-Zadeh S. et 672 al. (2018) Combinatorial targeting by microRNAs co-ordinates post- 673 transcriptional control of EMT. Cell Systems 7, 77-91.e7. 674 Deng F., Zhou R., Lin C., Yang S., Wang H., Li W. et al. (2019) Tumor-secreted 675 dickkopf2 accelerates aerobic glycolysis and promotes angiogenesis in colorectal 676 cancer. Theranostics 9, 1001-1014. 677 Friedman R.C., Farh K.K.H., Burge C.B. & Bartel D.P. (2009) Most mammalian mRNAs 678 are conserved targets of microRNAs. Genome Research 19, 92–105. 679 Gaidatzis D., Burger L., Florescu M. & Stadler M.B. (2015) Analysis of intronic and 680 exonic reads in RNA-seq data characterizes transcriptional and post-transcriptional 681 regulation. Nature Biotechnology 33, 722–729. 682 Grimson A., Farh K.K.H., Johnston W.K., Garrett-Engele P., Lim L.P. & Bartel D.P. 683 (2007) MicroRNA targeting specificity in mammals: determinants beyond seed 684 pairing. Molecular Cell 27, 91–105. 685 Guo Y., Liu J., Elfenbein S.J., Ma Y., Zhong M., Qiu C., Ding Y. et al. (2015) 686 Characterization of the mammalian miRNA turnover landscape. Nucleic Acids 687 Research 43, 2326–2341. 688 Han H., Gu S., Chu W., Sun W., Wei W., Dang X. et al. (2017) miR-17-5p regulates 689 differential expression of NCOA3 in pig intramuscular and subcutaneous adipose 690 tissue. Lipids 52, 939–949. 691 Horodyska J., Wimmers K., Reyer H., Trakooljul N., Mullen A.M., Lawlor P.G. et al. 692 (2018) RNA-seq of muscle from pigs divergent in feed efficiency and product 28 693 quality identifies differences in immune response, growth, and macronutrient and 694 connective tissue metabolism. BMC Genomics 19, 791. 695 Izquierdo A.G., Crujeiras A.B., Casanueva F.F. & Carreira M.C. (2019) Leptin, obesity, 696 and leptin resistance: where are we 25 years later? Nutrients 11, 2704. 697 Jacobsen M.J., Havgaard J.H., Anthon C., Mentzel C.M.J., Cirera S., Krogh P.M. et al. 698 (2019) Epigenetic and transcriptomic characterization of pure adipocyte fractions 699 from obese pigs identifies candidate pathways controlling metabolism. Frontiers 700 in Genetics 10, 1268. 701 Kim D., Paggi J.M., Park C., Bennett C. & Salzberg S.L. (2019) Graph-based genome 702 alignment and genotyping with HISAT2 and HISAT-genotype. Nature 703 Biotechnology 37, 907–915. 704 Kogelman L.J.A., Kadarmideen H.N., Mark T., Karlskov-Mortensen P., Bruun C.S., 705 Cirera S. et al. (2013) An F2 pig resource population as a model for genetic studies 706 of obesity and obesity-related diseases in humans: Design and genetic parameters. 707 Frontiers in Genetics 4, 29. 708 Kosacka J., Nowicki M., Paeschke S., Baum P., Blüher M. & Klö ting N. (2018) Up- 709 regulated autophagy: as a protective factor in adipose tissue of WOKW rats with 710 metabolic syndrome. Diabetology & Metabolic Syndrome 10, 13. 711 Kozomara A., Birgaoanu M. & Griffiths-Jones S. (2019) MiRBase: From microRNA 712 sequences to function. Nucleic Acids Research 47, D155–D162. 713 Kubo M., Ijichi N., Ikeda K., Horie-Inoue K., Takeda S. & Inoue S. (2009) Modulation 714 of adipogenesis-related gene expression by estrogen-related receptor γ during 715 adipocytic differentiation. Biochimica et Biophysica Acta (BBA) Gene Regulatory 716 Mechanisms 1789, 71–77. 717 Kuo C.-S., Chen J.-S., Lin L.-Y., Schmid-Schönbein G.W., Chien S., Huang P.-H. et al. 29 718 (2020) Inhibition of serine protease activity protects against high fat diet-induced 719 inflammation and insulin resistance. Scientific Reports 10, 1–11. 720 La Manno G., Soldatov R., Zeisel A., Braun E., Hochgerner H., Petukhov V. et al. (2018) 721 RNA velocity of single cells. Nature 560, 494–498. 722 Langmead B., Trapnell C., Pop M. & Salzberg S.L. (2009) Ultrafast and memory- 723 efficient alignment of short DNA sequences to the human genome. Genome 724 Biology 10, R25. 725 Lawrence M., Huber W., Pagès H., Aboyoun P., Carlson M., Gentleman R. et al. (2013) 726 Software for computing and annotating genomic ranges. PLoS Computational 727 Biology 9, e1003118. 728 Li X., Shan J., Chang W., Kim I., Bao J., Lee H.-J. et al. (2012) Chemical and genetic 729 evidence for the involvement of Wnt antagonist Dickkopf2 in regulation of glucose 730 metabolism. Proceedings of the National Academy of Sciences (USA) 109, 11402– 731 11407. 732 Liao Y., Smyth G.K. & Shi W. (2014) featureCounts: an efficient general purpose 733 program for assigning sequence reads to genomic features. Bioinformatics 30, 734 923–930. 735 Lindegaard B., Mathews V.B., Brandt C., Hojman P., Allen T.L., Estevez E. et al. (2013) 736 Interleukin-18 activates skeletal muscle AMPK and reduces weight gain and 737 insulin resistance in mice. Diabetes 62, 3064–3074. 738 Mann H.B. & Whitney D.R. (1947) On a test of whether one of two random variables is 739 stochastically larger than the other. Annals of Mathematical Statistics 18, 50–60. 740 Marco A. (2018) SeedVicious: Analysis of microRNA target and near-target sites. PLoS 741 One 13, e0195532. 742 Mármol-Sánchez E., Ramayo-Caldas Y., Quintanilla R., Cardoso T.F., González-Prendes 30 743 R., Tibau J. et al. (2020) Co-expression network analysis predicts a key role of 744 microRNAs in the adaptation of the porcine skeletal muscle to nutrient supply. 745 Journal of Animal Science & Biotechnology 11, 10. 746 Martin M. (2011) Cutadapt removes adapter sequences from high-throughput sequencing 747 reads. EMBnet.Journal 17, 10. 748 Martínez P., Gómez-López G., García F., Mercken E., Mitchell S., Flores J.M. et al. 749 (2013) RAP1 protects from obesity through its extratelomeric role regulating gene 750 expression. Cell Reports 3, 2059–2074. 751 Monteiro M.C., Sanyal M., Cleary M.L., Sengenès C., Bouloumié A., Dani C. et al. 752 (2011) PBX1: a novel stage-specific regulator of adipocyte development. Stem 753 Cells 29, 1837–1848. 754 Óvilo C., Benítez R., Fernández A., Núñez Y., Ayuso M., Fernández A.I. et al. (2014) 755 Longissimus dorsi transcriptome analysis of purebred and crossbred Iberian pigs 756 differing in muscle characteristics. BMC Genomics 15, 413. 757 Patel R.K., West J.D., Jiang Y., Fogarty E.A. & Grimson A. (2020) Robust partitioning 758 of microRNA targets from downstream regulatory changes. Nucleic Acids 759 Research 48, 9724–9746. 760 Pearen M.A., Goode J.M., Fitzsimmons R.L., Eriksson N.A., Thomas G.P., Cowin G.J. 761 et al. (2013) Transgenic muscle-specific Nor-1 expression regulates multiple 762 pathways that effect adiposity, metabolism, and endurance. Molecular 763 Endocrinology 27, 1897–1917. 764 Pérez-Montarelo D., Fernández A., Barragán C., Noguera J.L., Folch J.M., Rodríguez 765 M.C. et al. (2013) Transcriptional characterization of porcine leptin and leptin 766 receptor genes. PLoS One 8, e66398. 767 Perttilä J., Merikanto K., Naukkarinen J., Surakka I., Martin N.W., Tanhuanpää K. et al. 31 768 (2009) OSBPL10, a novel candidate gene for high triglyceride trait in dyslipidemic 769 Finnish subjects, regulates cellular lipid metabolism. Journal of Molecular 770 Medicine 87, 825–835. 771 Pilcher C.M., Jones C.K., Schroyen M., Severin A.J., Patience J.F., Tuggle C.K. et al. 772 (2015) Transcript profiles in longissimus dorsi muscle and subcutaneous adipose 773 tissue: A comparison of pigs with different postweaning growth rates. Journal of 774 Animal Science 93, 2134–2143. 775 Pillman K.A., Scheer K.G., Hackett-Jones E., Saunders K., Bert A.G., Toubia J. et al. 776 (2019) Extensive transcriptional responses are co-ordinated by microRNAs as 777 revealed by Exon–Intron Split Analysis (EISA). Nucleic Acids Research 47, 8606- 778 8619. 779 Reverter A. & Chan E.K.F. (2008) Combining partial correlation and an information 780 theory approach to the reversed engineering of gene co-expression networks. 781 Bioinformatics 24, 2491–2497. 782 Robinson M.D. & Oshlack A. (2010) A scaling normalization method for differential 783 expression analysis of RNA-seq data. Genome Biology 11, R25. 784 Robinson M.D., McCarthy D.J. & Smyth G.K. (2010) edgeR: a Bioconductor package 785 for differential expression analysis of digital gene expression data. Bioinformatics 786 26, 139–40. 787 Schaefke B., Sun W., Li Y.-S., Feng L. & Chen W. (2018) The evolution of 788 posttranscriptional regulation. Wiley Interdisciplinary Reviews 9, e1485. 789 Tamahara T., Ochiai K., Muto A., Kato Y., Sax N., Matsumoto M. et al. (2017) The 790 mTOR-Bach2 cascade controls cell cycle and class switch recombination during 791 B cell differentiation. Molecular and Cellular Biology 37, e00418-17. 792 Tao R., Xiong X., Liangpunsakul S. & Dong X.C. (2015) Sestrin 3 protein enhances 32 793 hepatic insulin sensitivity by direct activation of the mTORC2-Akt signaling. 794 Diabetes 64, 1211–1223. 795 Tarbier M., Mackowiak S.D., Frade J., Catuara-Solarz S., Biryukova I., Gelali E. et al. 796 (2020) Nuclear gene proximity and protein interactions shape transcript 797 covariations in mammalian single cells. Nature Communications 11, 5445. 798 Ussar S., Bezy O., Blüher M. & Kahn C.R. (2012) Glypican-4 enhances insulin signaling 799 via interaction with the insulin receptor and serves as a novel adipokine. Diabetes 800 61, 2289–2298. 801 Wang Y., Yan S., Xiao B., Zuo S., Zhang Q., Chen G. et al. (2018) Prostaglandin F 2α 802 facilitates hepatic glucose production through CaMKIIγ/p38/FOXO1 signaling 803 pathway in fasting and obesity. Diabetes 67, 1748–1760. 804 Warr A., Affara N., Aken B., Beiki H., Bickhart D.M., Billis K. et al. (2020) An improved 805 pig reference genome sequence to enable pig genetics and genomics research. 806 Gigascience 9, 1–14. 807 Watson-Haigh N.S., Kadarmideen H.N. & Reverter A. (2010) PCIT: an R package for 808 weighted gene co-expression networks based on partial correlation and 809 information theory approaches. Bioinformatics 26, 411–413. 810 Wu D., Hu D., Chen H., Shi G., Fetahu I.S., Wu F. et al. (2018) Glucose-regulated 811 phosphorylation of TET2 by AMPK reveals a pathway linking diabetes to cancer. 812 Nature 559, 637–641. 813 Xie S., Li X., Qian L., Cai C., Xiao G., Jiang S. et al. (2019) An integrated analysis of 814 mRNA and miRNA in skeletal muscle from myostatin-edited Meishan pigs. 815 Genome 62, 305–315. 816 Xu S., Ni H., Chen H. & Dai Q. (2019) The interaction between STAT3 and nAChRα1 817 interferes with nicotine-induced atherosclerosis via Akt/mTOR signaling cascade. 33 818 Aging 11, 8120–8138. 819 Yang J. & Shi B. yin (2021) Dickkopf (Dkk)-2 is a beige fat-enriched adipokine to 820 regulate adipogenesis. Biochemical and Biophysical Research Communications 821 548, 211–216. 822 Zhang S., Hulver M.W., McMillan R.P., Cline M.A. & Gilbert E.R. (2014) The pivotal 823 role of pyruvate dehydrogenase kinases in metabolic flexibility. Nutrition & 824 Metabolism 11, 10. 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 34 842 Table 1: mRNA genes with the top 5% post-transcriptional (PTc) scores and at least 2- 843 fold exonic fraction (ΔEx) reduction (equivalent to -1 in the log2 scale) of gluteus medius 844 skeletal muscle samples from fasting (AL-T0, N = 11) and fed (AL-T2, N = 12) Duroc 845 gilts. miRNA ID Gene log2FCa ΔExb PTcc P-value q-value DEd target ENSSSCG00000032094 DKK2 -2.010 -1.431 -4.738 1.654E-05 3.830E-03 x ENSSSCG00000015334 PDK4 -2.108 -5.250 -4.698 4.693E-03 1.330E-01 x x ENSSSCG00000015037 IL18 -1.655 -1.191 -3.682 4.787E-03 1.340E-01 x x ENSSSCG00000005385 NR4A3 -1.337 -3.082 -3.646 4.038E-02 4.098E-01 x x ENSSSCG00000003766 DNAJB4 -1.391 -1.008 -3.348 8.358E-03 1.905E-01 x ENSSSCG00000015969 CHRNA1 -1.561 -1.339 -3.341 2.606E-03 9.406E-02 x x ENSSSCG00000039419 SLCO4A1 -1.055 -2.279 -3.180 2.820E-02 3.544E-01 x x ENSSSCG00000049158 -1.107 -1.096 -3.164 3.182E-02 3.735E-01 ENSSSCG00000004347 FBXL4 -1.298 -1.126 -3.133 1.422E-03 6.520E-02 x x ENSSSCG00000004979 MYO9A -1.239 -1.003 -3.043 7.296E-03 1.731E-01 x ENSSSCG00000013351 NAV2 -1.163 -1.196 -2.863 2.605E-04 2.301E-02 x x ENSSSCG00000032741 TBC1D9 -0.913 -1.061 -2.736 1.534E-02 2.583E-01 x ENSSSCG00000031728 ABRA -1.238 -1.393 -2.704 1.295E-03 6.116E-02 x x ENSSSCG00000006331 PBX1 -0.891 -1.039 -2.480 1.135E-02 2.177E-01 x ENSSSCG00000035037 SIK1 -1.357 -1.289 -2.475 3.999E-03 1.212E-01 x x ENSSSCG00000038374 CIART -1.027 -1.321 -2.052 1.543E-02 2.587E-01 x ENSSSCG00000023806 LRRN1 -0.776 -1.013 -1.983 1.580E-01 7.074E-01 x ENSSSCG00000009157 TET2 -0.381 -1.123 -1.792 4.880E-01 9.582E-01 x ENSSSCG00000011133 PFKFB3 -0.022 -2.256 -1.785 9.712E-01 9.987E-01 x x ENSSSCG00000002283 FUT8 -0.578 -1.286 -1.784 9.887E-02 6.059E-01 x x ENSSSCG00000023133 OSBPL6 -0.432 -1.088 -1.772 3.835E-01 9.108E-01 x ENSSSCG00000017986 NDEL1 -0.767 -1.644 -1.759 1.006E-02 2.081E-01 x x ENSSSCG00000031321 NR4A1 -0.630 -1.328 -1.720 6.298E-02 5.006E-01 x ENSSSCG00000035101 KLF5 -0.519 -1.487 -1.708 2.942E-01 8.488E-01 x x ENSSSCG00000004332 BACH2 -0.714 -2.105 -1.705 9.089E-02 5.861E-01 x x ENSSSCG00000017983 PER1 -0.773 -1.073 -1.627 3.000E-02 3.662E-01 x 846 847 ENSSSCG00000049158 did not have any annotated 3´-UTR so it was excluded from further analyses. aLog2FC: 848 estimated log2 fold change for mean exonic fraction from gluteus medius skeletal muscle samples of fasted AL-T0 and 849 fed AL-T2 Duroc gilts; bΔEx: exonic fraction increment (Ex2 – Ex1) in log2 scale when comparing exon abundances in 850 AL-T0 (Ex1) vs AL-T2 (Ex2) Duroc gilts; cPTc: post-transcriptional signal (ΔEx – ΔInt) after z-score normalization. The 851 q-value has been calculated with the false discovery rate approach (Benjamini & Hochberg 1995). dDE: Significantly 852 differentially expressed (|FC| > 2; q-value < 0.05). The “x” symbols indicate significantly downregulated genes 853 according to their exonic counts, as well as those mRNA genes targeted by at least one of the significantly upregulated 854 miRNAs excluding redundant seeds (N = 6, Table S5). 35 855 Table 2: mRNA genes with the top 5% post-transcriptional (PTc) scores and at least 2- 856 fold exonic fraction (ΔEx) reduction (equivalent to -1 in the log2 scale) of adipocytes 857 from lean (N = 5) and obese (N = 5) UNIK minipigs. mRNA Gene log2FCa ΔExb PTcc P-value q-value DEd miRNA target ENSSSCG00000010814 ESRRG -0.591 -5.305 -6.425 7.364E-01 9.996E-01 x ENSSSCG00000037015 SESN3 -2.000 -1.378 -5.707 2.541E-07 2.477E-03 x x ENSSSCG00000032452 WFS1 -2.198 -2.138 -5.510 9.509E-03 9.996E-01 x ENSSSCG00000039548 PTGFR -1.634 -1.590 -4.915 8.804E-03 9.996E-01 x ENSSSCG00000013829 SYDE1 -1.670 -1.160 -4.188 5.795E-04 6.000E-01 x ENSSSCG00000002265 FAM174B -1.244 -1.726 -4.179 5.385E-02 9.996E-01 x ENSSSCG00000016233 SERPINE2 -1.735 -2.060 -3.603 5.684E-02 9.996E-01 x x ENSSSCG00000006243 PENK -0.420 -2.104 -3.573 7.628E-01 9.996E-01 ENSSSCG00000038879 RELB -1.272 -1.056 -3.512 3.659E-03 9.996E-01 x ENSSSCG00000023408 SAMD4A -1.328 -1.156 -3.509 4.486E-02 9.996E-01 x ENSSSCG00000008449 SLC3A1 -1.014 -1.154 -3.491 5.859E-02 9.996E-01 ENSSSCG00000014921 PRSS23 -1.141 -1.739 -3.360 2.719E-01 9.996E-01 x ENSSSCG00000017186 RNF157 -1.218 -2.338 -3.317 2.413E-01 9.996E-01 x x ENSSSCG00000035403 RFX2 -1.109 -1.022 -2.958 1.550E-01 9.996E-01 x ENSSSCG00000010893 -0.655 -1.352 -2.931 4.068E-01 9.996E-01 x ENSSSCG00000031819 TP53I11 -1.002 -1.711 -2.883 4.102E-01 9.996E-01 x ENSSSCG00000017137 METRNL -0.674 -1.102 -2.812 2.422E-01 9.996E-01 ENSSSCG00000032562 TMC6 -0.837 -1.152 -2.765 2.078E-01 9.996E-01 ENSSSCG00000031261 RHOQ -0.903 -1.046 -2.750 1.839E-02 9.996E-01 ENSSSCG00000001089 GPLD1 -0.872 -1.761 -2.723 4.302E-01 9.996E-01 x ENSSSCG00000034259 PMEPA1 -0.880 -1.348 -2.720 3.583E-01 9.996E-01 x ENSSSCG00000017014 PANK3 -0.614 -1.037 -2.557 2.288E-01 9.996E-01 x ENSSSCG00000003377 ACOT7 -0.790 -2.688 -2.544 3.439E-01 9.996E-01 x ENSSSCG00000010079 PPM1F -0.762 -1.035 -2.473 4.967E-02 9.996E-01 x x ENSSSCG00000040464 LEP -0.747 -2.186 -2.463 1.880E-01 9.996E-01 x x ENSSSCG00000022029 RAP1GAP -0.120 -1.109 -2.418 8.822E-01 9.996E-01 x ENSSSCG00000022099 TP53INP2 -0.628 -1.058 -2.403 3.683E-01 9.996E-01 ENSSSCG00000025652 CDH1 -0.472 -2.592 -2.372 6.533E-01 9.996E-01 x ENSSSCG00000027266 PNPLA3 -0.443 -1.386 -2.198 5.725E-01 9.996E-01 ENSSSCG00000015401 PCLO -0.674 -1.492 -2.182 4.537E-01 9.996E-01 x ENSSSCG00000020872 -1.029 -1.128 -2.090 1.340E-01 9.996E-01 ENSSSCG00000032633 FAM53A -0.749 -1.033 -2.066 4.576E-02 9.996E-01 x ENSSSCG00000015559 NCF2 -0.679 -1.221 -2.061 3.570E-01 9.996E-01 x ENSSSCG00000015766 WDR17 -0.609 -1.139 -1.998 2.093E-01 9.996E-01 ENSSSCG00000009761 NCOR2 -0.681 -1.421 -1.913 2.643E-01 9.996E-01 ENSSSCG00000016928 RAB3D -0.491 -1.142 -1.888 2.953E-01 9.996E-01 ENSSSCG00000011230 OSBPL10 -0.576 -1.594 -1.869 4.272E-01 9.996E-01 x ENSSSCG00000017298 TANC2 -0.615 -1.541 -1.846 4.896E-01 9.996E-01 ENSSSCG00000007899 -0.524 -1.036 -1.814 4.500E-01 9.996E-01 x ENSSSCG00000026421 PKD2L2 -0.463 -1.230 -1.800 5.098E-01 9.996E-01 ENSSSCG00000015332 PON1 -0.626 -1.076 -1.763 2.530E-01 9.996E-01 x x 36 ENSSSCG00000009215 ABCG2 -0.455 -1.446 -1.749 5.684E-01 9.996E-01 ENSSSCG00000017328 ARHGAP27 -0.235 -2.788 -1.699 8.113E-01 9.996E-01 x x ENSSSCG00000017199 TRIM47 -0.362 -1.057 -1.645 6.717E-01 9.996E-01 x 858 859 ENSSSCG00000016928 did not have an annotated 3´-UTR and was therefore excluded from further analyses. 860 aLog FC: 2 estimated log2 fold change for mean exonic fraction from adipocytes of lean and obese UNIK minipigs; bΔEx: 861 exonic fraction increment (Ex2 – Ex1) in log2 scale when comparing exon abundances in obese (Ex1) vs lean (Ex2) 862 UNIK minipigs; cPTc: post-transcriptional signal (ΔEx – ΔInt) after z-score normalization. The q-value has been 863 calculated with the false discovery rate (FDR) approach (Benjamini & Hochberg 1995); dDE: Significantly 864 differentially expressed (|FC| > 2; q-value < 0.05). The “x” symbols indicate significantly downregulated genes 865 according to their exonic counts, as well as those mRNA genes targeted by at least one of the significantly upregulated 866 miRNAs excluding redundant seeds (N = 4, Table S10). 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 37 887 Figures 888 889 Figure 1: (a) Scatterplot depicting genes expressed in gluteus medius skeletal muscle of 890 fasted (AL-T0, N = 11) and fed (AL-T2, N = 12) Duroc gilts according to their exonic 891 fraction (ΔEx) and post-transcriptional (PTc) scores. Genes with the top 5% negative PTc 892 scores and at least 2-fold ΔEx reduction (equivalent to -1 in the log2 scale) are highlighted 893 in purple and delimited by dashed lines. (b) Enrichment analyses comparing all expressed 894 mRNA genes and the set of mRNA genes with the top 5% negative PTc scores and at 895 least 2-fold ΔEx reduction as being putatively targeted by either significantly upregulated 896 miRNAs (FC > 1.5; q-value < 0.05), significantly downregulated miRNAs (FC < -1.5; q- 897 value < 0.05) or the top 5% most highly expressed miRNAs, excluding significantly 898 upregulated miRNAs. As indicated with the dashed line, a nominal P-value = 0.05 was 899 set as a significance threshold. (c) Scatterplot depicting genes expressed in gluteus medius 38 900 skeletal muscle of fasted (AL-T0, N = 11) and fed (AL-T2, N = 12) Duroc gilts according 901 to their exonic fraction (ΔEx) and post-transcriptional (PTc) scores. Genes significantly 902 upregulated are in green, while those being downregulated are in red (|FC| > 2; q-value < 903 0.05). (d) Covariation enrichment scores (CES) for the exonic and intronic fractions of 904 the mRNA genes with the top 5% negative PTc scores and at least 2-fold ΔEx reduction 905 predicted to harbor binding sites for upregulated miRNAs (N = 6) in the gluteus medius 906 skeletal muscle of fasted (AL-T0, N = 11) and fed (AL-T2, N = 12) Duroc gilts (Tables 1 907 and S8). The control set was established by generating 1,000 permuted lists of 21 genes 908 chosen at random and using their exonic and intronic fractions for calculating their CES 909 values. Statistical significance was assessed using a Mann-Whitney U non-parametric test 910 (Mann & Whitney 1947). The dashed line represents a CES of 1, equivalent to an 911 observed null fold change in covariation. 912 913 914 915 916 917 918 919 39 920 921 Figure 2: (a) Scatterplot depicting genes expressed in adipocytes obtained from UNIK 922 minipigs with lean (N = 5) and obese (N = 5) phenotypes (using their body mass index 923 as reference) according to their exonic fraction (ΔEx) and post-transcriptional (PTc) 924 scores. Genes with the top 5% negative PTc scores and at least 2-fold ΔEx reduction 925 (equivalent to -1 in the log2 scale) are highlighted in purple and delimited by dashed lines. 926 (b) Enrichment analyses comparing all expressed mRNA genes and the set of mRNA 927 genes with the top 5% negative PTc scores and at least 2-fold ΔEx reduction as being 928 putatively targeted by either significantly upregulated miRNAs (FC > 1.5; q-value < 929 0.05), significantly downregulated miRNAs (FC < -1.5; q-value < 0.05) or the top 5% 930 most highly expressed miRNAs, excluding significantly upregulated miRNAs. As 931 indicated with the dashed line, a nominal P-value = 0.05 was set as a significance 932 threshold. (c) Scatterplot depicting genes expressed in adipocytes obtained from UNIK 933 minipigs with lean (N = 5) and obese (N = 5) phenotypes (using their body mass index 40 934 as reference) according to their exonic fraction (ΔEx) and post-transcriptional (PTc) 935 scores. Genes significantly upregulated are in green, while those being downregulated are 936 in red (|FC| > 2; q-value < 0.05). (d) Covariation enrichment scores (CES) for the exonic 937 and intronic fractions of the mRNA genes with the top 5% negative PTc scores and at 938 least 2-fold ΔEx reduction predicted to harbor binding sites for upregulated miRNAs (N 939 = 4) in adipocytes obtained from UNIK minipigs with lean (N = 5) and obese (N = 5) 940 phenotypes (Tables 2 and S13). The control set was established by generating 1,000 941 permuted lists of 25 genes chosen at random and using their exonic and intronic fractions 942 for calculating their CES values. Statistical significance was assessed using a Mann- 943 Whitney U non-parametric test (Mann & Whitney 1947). The dashed line represents a 944 CES of 1, equivalent to an observed null fold change in covariation. 945 946 947 Supplementary Materials 948 Figure S1: Diagram depicting the routine/pipeline implemented for studying miRNA- 949 driven post-transcriptional regulatory signals applying the EISA approach and additional 950 enrichment and covariation analyses. 951 Figure S2: Diagram representing each one of the context-based filtering criteria used for 952 excluding unreliable in silico-predicted miRNA-mRNA interactions. AU: miRNA 953 binding sites with AU-rich flanking sequences (30 nts upstream and downstream). M: 954 miRNA binding sites located in the middle of the 3’-UTR sequence (45-55%). E: miRNA 955 binding sites located too close (< 15 nts) to the beginning or the end of the 3’-UTR 956 sequences. 957 Figure S3: Scatterplots depicting the exonic (ΔEx) and intronic (ΔInt) fractions of 958 expressed genes from gluteus medius skeletal muscle samples of fasting (AL-T0, N = 11) 41 959 and fed (AL-T2, N = 12) Duroc gilts. (a) Genes differentially expressed and showing 960 either significant upregulation (FC > 2; q-value < 0.05, in green) or downregulation (FC 961 < -2, q-value < 0.05, in red) in fed (AL-T2, N = 12) Duroc gilts with respect to their fasted 962 (AL-T0, N = 11) counterparts. (b) Genes with the top 5% negative post-transcriptional 963 (PTc) scores and at least 2-fold reduced exonic (ΔEx) fraction (equivalent to -1 in the log2 964 scale) are highlighted in purple. 965 Figure S4: Enrichment analyses comparing all expressed mRNA genes and the set of 966 mRNA genes with the (a) top 1% and (b) top 5% negative PTc scores and at least 2-fold 967 ΔEx reduction as being putatively targeted by significantly upregulated miRNAs (FC > 968 1.5; q-value < 0.05) from gluteus medius skeletal muscle samples of fasting (AL-T0, N = 969 11) and fed (AL-T2, N = 12) Duroc gilts. Results show the change in enrichment 970 significance (expressed as -log10 of the estimated P-value) when incorporating context- 971 based pruning of 8mer, 7mer-m8 and 7mer-A1 miRNA binding sites. R: Raw enrichment 972 analyses without any additional context-based pruning. AU: Enrichment analyses 973 removing miRNA binding sites without AU-rich flanking sequences (30 nts upstream and 974 downstream). M: Enrichment analyses removing miRNA binding sites located in the 975 middle of the 3’-UTR sequence (45-55%). E: Enrichment analyses removing miRNA 976 binding sites located too close (< 15 nts) to the beginning or the end of the 3’-UTR 977 sequences. The dashed line represents a nominal P-value of 0.05 set as the significance 978 threshold. 979 Figure S5: Quantification of selected genes and miRNAs expressed in the pig adipose 980 tissue by qPCR. (a) Barplots depicting qPCR log2 transformed relative quantities (Rq) for 981 LEP, OSBPL10, PRSS23, RNF157 and SERPINE2 mRNA transcripts measured in 982 adipocytes from the retroperitoneal fat of lean (N = 5) and obese (N = 5) UNIK minipigs. 983 (b) Barplots depicting qPCR log2 transformed relative quantities (Rq) for ssc-miR-148a- 42 984 3p, ssc-miR-214-3p and ssc-miR-92b-3p miRNA transcripts measured in isolated 985 adipocytes from the retroperitoneal fat of lean (N = 5) and obese (N = 5) UNIK minipigs. 986 All the analyzed mRNA genes showed a reduced expression in lean pigs compared with 987 their obese counterparts, and the LEP gene was the most significantly downregulated 988 gene. For miRNAs, the opposite pattern of expression was observed, being all of them 989 upregulated in lean minipigs. Moreover, ssc-miR-92b-3p showed the most significant 990 increased expression in lean minipigs, in agreement with results obtained in differential 991 expression analyses (Table S10). 992 993 Table S1: Phenotypic values of body mass index (BMI) trait and sex classification 994 recorded in 11 Duroc-Göttingen minipigs from the F2-UNIK resource population. 995 Table S2: Primers for qPCR verification of selected mRNAs and miRNAs in the F2- 996 UNIK Duroc-Göttingen minipig population. 997 Table S3: Raw Cq values obtained in qPCR analyses. 998 Table S4: Differential expression analyses of RNA-seq data using the edgeR tool and 999 comparing gluteus medius expression profiles of fasted AL-T0 (N = 11) and fed AL-T2 1000 (N = 12) Duroc gilts. (a) Differentially expressed genes (q-value < 0.05, N = 454). In bold 1001 are genes either upregulated (N = 52) or downregulated (N = 80) with |FC| > 2 and q- 1002 value < 0.05. (b) Differential expression results for genes with top 5% post-transcriptional 1003 (PTc) scores and at least 2-fold reduced exonic fraction (ΔEx) (equivalent to -1 in the log2 1004 scale). The 18 genes with FC < -2 and q-value < 0.05 are shown in bold. (c) Differential 1005 expression results showing the whole list of expressed genes with an average expression 1006 above 1 CPM in at least 50% of samples within each group (N = 9,492). 43 1007 Table S5: Differential expression analyses of microRNAs using the edgeR tool and 1008 comparing gluteus medius expression levels of fasted AL-T0 (N = 11) and fed AL-T2 (N 1009 = 12) Duroc gilts. 1010 Table S6: Post-transcriptional (PTc) and transcriptional (Tc) signals detected with EISA 1011 in genes expressed in gluteus medius skeletal muscle samples from fasted (AL-T0, N = 1012 11) and fed (AL-T2, N = 12) Duroc gilts. (a) Genes with significant PTc scores (|FC| > 2; 1013 q-value < 0.05, N = 133). In bold are genes with at least 2-fold reduced ΔEx fraction (N 1014 = 3). (b) Genes with significant Tc scores (|FC| > 2; q-value < 0.05, N = 344). In bold are 1015 genes also significant in their PTc scores (|FC| > 2; q-value < 0.05, N = 91). EISA results 1016 showing the whole list of genes with an average expression above 1 CPM in at least 50% 1017 of samples within each group (N = 9,492) and their (c) PTc and (d) Tc scores. 1018 Table S7: Binding sites in the 3’-UTRs of mRNA genes (with the top 5% negative PTc 1019 scores and at least 2-fold reduction in the exonic fraction) predicted as targets (N = 21) 1020 of non-redundant significantly upregulated miRNAs (N = 6) expressed in the gluteus 1021 medius skeletal muscle samples from fasting (AL-T0, N = 11) and fed (AL-T2, N = 12) 1022 Duroc gilts. 1023 Table S8: Covariation enrichment scores (CES) for the exonic and intronic fractions of 1024 mRNA genes (with the top 5% negative post-transcriptional signals PTc and at least 2- 1025 fold reduction in their exonic ΔEx fraction) predicted as targets of non-redundant 1026 significantly upregulated miRNAs (N = 6) expressed in gluteus medius skeletal muscle 1027 samples from fasting (AL-T0, N = 11) and fed (AL-T2, N = 12) Duroc gilts. 1028 Table S9: Differential expression analyses of RNA-seq data using the edgeR tool and 1029 comparing adipocyte expression profiles of lean (N = 5) and obese (N = 5) UNIK 1030 minipigs. (a) Differentially expressed genes (q-value < 0.05, N = 299). In bold are genes 1031 either upregulated (N = 52) or downregulated (N = 95) with |FC| > 2 and q-value < 0.05. 44 1032 (b) Differential expression results for genes with top 5% post-transcriptional (PTc) scores 1033 and at least 2-fold reduced exonic fraction (ΔEx) (equivalent to -1 in the log2 scale). The 1034 12 genes with FC < -2 and q-value < 0.05 are shown in bold. (c) Differential expression 1035 results showing the whole list of expressed genes with an average expression above 1 1036 CPM in at least 50% of samples within each group (N = 9,746). 1037 Table S10: Differential expression analyses of microRNAs using the edgeR tool and 1038 comparing adipocyte expression profiles from lean (N = 5) and obese (N = 5) UNIK 1039 minipigs. 1040 Table S11: Post-transcriptional (PTc) and transcriptional (Tc) signals detected with EISA 1041 in genes expressed in adipocytes from lean (N = 5) and obese (N = 5) UNIK minipigs. 1042 (a) Genes with significant PTc scores (|FC| > 2; q-value < 0.05, N = 1). In bold are genes 1043 with at least 2-fold reduced ΔEx fractions (N = 1). (b) Genes with significant Tc scores 1044 (|FC| > 2; q-value < 0.05, N = 195). In bold are genes also among the top 5% negative 1045 PTc scores and at least 2-fold ΔEx reduction (N = 3). EISA results showing the whole list 1046 of genes with an average expression above 1 CPM in at least 50% of samples within each 1047 group (N = 9,746) and their (c) PTc and (d) Tc scores. 1048 Table S12: Binding sites in the 3’-UTRs of mRNA genes (with the top 5% negative PTc 1049 scores and at least 2-fold reduction in the exonic fraction) predicted as targets (N = 25) 1050 of non-redundant significantly upregulated miRNAs (N = 4) expressed in adipocytes from 1051 lean (N = 5) and obese (N = 5) UNIK minipigs. 1052 Table S13: Covariation enrichment scores (CES) for the exonic and intronic fractions of 1053 mRNA genes (with the top 5% negative post-transcriptional signals PTc and at least 2- 1054 fold reduction in their exonic ΔEx fraction) predicted as targets of non-redundant 1055 significantly upregulated miRNAs (N = 4) expressed in adipocytes from lean (N = 5) and 1056 obese (N = 5) UNIK minipigs. 45