miRNAs regulate acute transcriptional changes in broiler embryos in response to modification of incubation temperature

MicroRNAs are post-transcriptional regulators that play critical roles in diverse biological processes. We hypothesize that miRNAs may be involved in regulating transcriptome responses to changes in embryonic incubation temperature in chickens affecting differentiation and proliferation processes during tissue development. Therefore, we conducted comparative transcriptome profiling of miRNAs to examine altered expression in breast and hind muscle of embryos and day 35 chickens experiencing high (38.8 °C), control (37.8 °C), or low (36.8 °C) embryonic incubation temperature during embryonic day (ED) 7–10 or ED10–13. The results revealed differential expression of miRNAs due to modification of embryonic incubation temperature in a muscle type-specific and a developmental stage-specific manner. The immediate effects of thermal change observed in embryos were substantial compared to the subtle long-term effects in chickens at day 35 post-hatch. Upregulation of miR-133 in breast muscle and downregulation of miR-199a-5p, miR-1915, and miR-638 in hind muscle post ED7–10 high-temperature treatment are functionally associated with myogenesis and body size. ED10–13 low-temperature treatment led to downregulation of let-7, miR-93, and miR-130c that are related to proliferation and differentiation. The results provide insight into the dynamics of miRNA expression at variable embryonic incubation temperatures during developmental processes and indicate a major regulatory role of miRNAs in acute responses to modified environmental conditions that affect remodelling of cells and tissues.

Embryonic incubation temperature is a key factor for optimal physiological and developmental processes that may a have long-term influence on adult chickens. Incubation temperature profoundly influences physiological responses via alteration of biochemical reaction rates and protein structures as well as catalytic enzyme functions 1 . Within a limited range, it is critical for broilers to optimize body temperature during pre-and post-hatch processes 2 .
Manipulation of incubation temperature during specific stages of development can result in immediate transcriptomic changes in embryos, although changing temperature beyond critical thresholds can be lethal. Previous studies demonstrated that high temperatures during embryonic day (ED) 7-10 positively associate with improvement of slaughter and breast muscle weights in male broilers, but do not influence meat quality 3 . Our previous experiments showed acute and long-term transcriptomic changes with temperature manipulation during muscle fibre formation. Thermal incubation treatments influence several biological functions and pathways depending Siloam Springs, USA) and equally randomly assigned 1,001 hatching eggs to 6 experimental groups. All eggs were incubated in commercial incubators with 12 automatic turns per day (HEKA, Euro-Lux, Riet berg, Germany) at 37.8 °C with 55% relative humidity (RH) until three days prior to hatch. At ED7 all eggs were candled and unfertilized eggs and dead embryos (in total 108 eggs) were removed. During the last days of incubation RH was adjusted to 65% and turning was omitted until hatching (control conditions; continuously for control groups). The experimental thermal profiles comprised of transient adjusting temperature to 38.8 °C (high temperature) or 36.8 °C (low temperature) at early development (ED7-10) or late muscle development (ED10-13), respectively. Accordingly, the experimental design delivered 6 experimental groups in total (38.8 °C, 65% RH, ED7-10 (H10); 2) 37.8 °C, 55% RH, ED7-10 (C10); 3) 36.8 °C, 55% RH, ED7-10 (L10); 4) 38.8 °C, 65% RH, ED10-13 (H13); 5) 37.8 °C, 55% RH, ED10-13 (C13); 6) 36.8 °C, 55% RH, ED10-13 (L13) with 2 repeated batches each. The hatchlings were reared in barn system and fed a standard diet ad libitum until day 35 (D35; slaughter) (Supplementary Table S1). Broilers were slaughtered at experimental poultry abattoir of the Department of Animal Sciences (Goettingen, Germany). The processes were proceed by electric stunning (0.12 A, 5 to 10 sec), bleeding, scalded (58 to 60•C, 3 min) and defeathered. Carcasses were manually eviscerated then dissected and weighted. For each group, breast muscle (M. pectoralis) and hind muscle (M. gastrocnemius) tissue samples were collected at the respective embryonic stages (ED10 or ED13) and at D35 post-hatch.
Tissue samples were immediately dissected, snap frozen in liquid nitrogen, and stored at −80 °C until use. Zootechnical and biochemical traits were examined as previously described 13 . All animals were sexed, and 6-8 (sex-balanced) animals per experimental group at ED10, ED13, or D35 were used for expression analyses. Study design and sample collection procedures were approved by the Institutional Animal Care and Use Committee (IACUC) of the Department of Animal Sciences of the University of Goettingen, Germany and the Leibniz Institute for Farm Animal Biology and conducted according to the guidelines of the German Law of Animal Protection and the "EU Directive 2010/63/EU for animal experiments".
Small RNA isolation. Total RNA was isolated from individual samples (n =144; 6 biological replicates × 6 treatment groups × 2 muscle tissues at ED10 or ED13; 8 biological replicates × 6 treatment groups × 2 muscle tissues at D35) using Tri-Reagent (Sigma-Aldrich, Taufkirchen, Germany), and the small RNA fraction was retained using miReasy and RNeasy MinElute Cleanup kits (Qiagen, Hilden, Germany) with an on-column DNase treatment according to the manufacturer's protocol. RNA integrity was assessed by capillary electrophoresis using the Agilent Small RNA kit and an Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany). Total RNA concentration was determined using a NanoDrop ND-1000 spectrophotometer (PEQLAB, Erlangen, Germany). Additionally, absence of trace DNA contamination was verified by PCR amplification of Data processing and statistical analysis. Raw probe signal intensity was pre-processed and normalized using Perfect Matched and Detection Above Background features of Affymetrix Expression Consol software. Data were submitted to the MIAME-compliant database Gene Expression Omnibus (accession number: GSE83703-GSE83704) accessible via the National Center for Biotechnology Information (www.ncbi.nlm.nih.gov/ geo). Differential expression of miRNAs was computed by analysis of variance (JMP Genomics, SAS-Institute, Cary, NC, USA). Independent calculation was performed for each tissue. Fixed effects of temperature, treatment period, and their interactions were modelled in statistical tests. For analysis in D35 samples, slaughter weight was used as covariance in the statistical model. Differentially expressed miRNAs were identified by comparing treatment groups and corresponding controls for ED7-10 or ED10-13. Significance threshold was set at nominal p ≤ 0.05. To account for multiple testing, FDR values were calculated based on the p-value distribution using an approach designed for large data sets 14 . Functional miRNAs and potential target genes. To identify functional miRNAs and potential target genes, we integrated miRNA expression data from the present study and mRNA expression data of samples from our previous publications 4,5 using correlation analysis and target prediction (Fig. 1). Firstly, correlation analyses between signal intensities of differentially expressed miRNAs and all mRNA probe-sets were calculated. All significant negative correlations between miRNA-mRNA pairs were retained (p ≤ 0.05; 5% FDR). Secondly, all differentially expressed miRNAs were scanned for potential target genes against all available chicken mRNA sequences in the NCBI database using Target Scan software 15 . Predicted targets were further filtered using RNA Hybrid software 16 with an energy threshold cut-off of ≤ -25 kcal/mole. Accordingly, a "functional" miRNA was defined as a differentially expressed miRNA that negatively correlated with mRNA transcriptional level and was predicted to bind at the respective target genes. Pathway analysis. Differentially expressed miRNAs and their potential targets were mined for biological functions and gene regulatory networks using Ingenuity Pathway Analysis (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis). Statistical significance was determined based on Fisher's exact test (p ≤ 0.05). Significant biological pathways for each list of miRNAs and target genes were aggregated into new major categories based on shared GO subterms to simplify results, while maintaining a comprehensive view of biological processes.
Eight major biological functional groups were defined: (1) cell maintenance, proliferation, differentiation, and replacement; (2) organismal, organ, and tissue development; (3) nutrient metabolism; (4) genetic information and nucleic acid processing; (5) molecular transport; (6) cell signalling and interaction; (7) small molecule biochemistry; and (8) response to stimuli. Significant pathways were further considered "activated" or "deactivated" based on positive or negative Ingenuity z-scores that predict the direction of regulation of a pathway. Selected genes were used to generate regulatory networks based on best z-scores, p-values, and biological functions related to tissue development and myogenesis.
Quantitative real time PCR (qPCR) validation. Six miRNAs (gga-miR-199, gga-miR-30d, gga-miR-460a, gga-miR-100, gga-miR-222, gga-miR-133a) differentially expressed in embryo breast were validated by qPCR of the same individual samples that were used for miRNA-chips (n = 32). The cDNA was synthesized from 250 ng isolated miRNAs using Mir-X ™ miRNA First-Strand Synthesis and SYBR ® qRT-PCR Kit (Clontech Laboratories, Mountain View, USA) according to manufacturer's protocols. All measurements were performed in duplicates. The thermal parameters were 50 °C for 2 min, 95 °C for 2 min, followed by 45 cycles of 95 °C for 15 sec and 60 °C for 1 min. The universal qPCR primer was provided in the kit and the miRNA-specific forward primers were designed for the miRNAs (Supplementary Table S2). Internal standard of Ce_miR-39_1 miScript Primer Assay (QIAGEN, Germany) was used for normalize the miRNA expression value. The designed primer sequence information is accessible in Additional file. Correlation coefficient analysis between the miChip and qPCR data was calculated using SAS 9.3 (SAS Institute).

Results
Differentially expressed miRNAs. Although the Affymetrix GeneChip miRNA 3.0 array contains multiple mature miRNAs from diverse species that may resemble a chicken miRNA family, we treated each mature miRNA probe-set as an entity (feature) in statistical tests and then aggregated significant probe-sets into unique mature miRNAs. The number of differential probe-sets and mature miRNAs (unique miRNAs) for each tissue type, temperature treatment condition, and sampling stage are summarized in Table 1.
Numerous differentially expressed mature miRNAs were detected in all comparisons. In analyses comprising breast muscle (embryonic stage) and hind muscle (embryonic stage, adult stage) the threshold level for nominal p-values corresponded to FDR values of 0.18; for analyses comprising the breast muscle at adult stage, nominal p-values corresponded to a FDR value of 0.71 (Table 1). Quantitative real time PCR of exemplarily selected miR-NAs in breast muscle revealed significant correlations with microarray data (Supplementary Figure S1). The lists of miRNAs differentially expressed across comparisons of temperature treatments, tissue types, and sampling stages, with a comparable number of upregulated and downregulated miRNAs, offer a rich source for subsequent filtering based on the biological context and function-related criteria. The results suggest that volatile miRNA changes in terms of direction of regulation and timing regulate transcriptional alterations due to modification of embryonic incubation temperature.

Functional miRNAs and potential target genes.
Integrating miRNA expression data with previous gene expression profiles of matched samples (www.ncbi.nlm.nih.gov/geo; accession number: GSE76670) using correlation analysis revealed "functional" miRNAs with negatively correlated miRNA-mRNA relationships. From our previous study, we found that increasing temperature from 37.8 °C to 38.8 °C during ED7-10 (H10) and decreasing temperature from 37.8 °C to 36.8 °C during ED10-13 (L13) resulted in considerable immediate transcriptomic changes (based on the abundance of differentially expressed genes) in embryos, whereas decreasing temperature during ED7-10 (L10) as well as ED10-13 (L13) showed large long-term effects in D35 chickens 4,5 . Therefore, we focused on these treatment conditions that showed most prominent changes of the mRNA expression.
Numbers of miRNA-mRNA pairs, covering potential miRNA-targeted genes negatively correlated with miRNAs and also predicted as miRNA binding sites, and unique miRNAs are presented in Table 2. Further information on criteria and pipelines are shown in Fig. 1. Overall, embryos had higher numbers of potential miRNA-mRNA relationships, ranging from 104 miRNA-mRNA pairs for L13 in hind muscle to 941 pairs for H10 in breast muscle, compared to D35 chickens (2 pairs for L13 in hind muscle and 19 pairs for L10 in breast muscle).
For H10, 421 unique genes were negatively correlated with 40 miRNAs and 200 unique genes were predicted as target candidates for 38 miRNAs in breast and hind muscles, respectively. Number of miRNA-mRNA pairs in breast muscle was higher than in hind muscle at in-ovo stages. A similar number of miRNAs and potential targets were also identified in L13 in both muscle types. Only up to 10 miRNA-mRNA pairs were found among treatment conditions and muscle types at D35 (Table 2). Differential miRNAs and potential target genes were further used to generate a hierarchical clustering based on expression level to demonstrate an overall dominantly negative correlation between miRNA and mRNA expression levels (Fig. 2). Altogether, these results suggest that miRNAs may play an essential regulatory role on the immediate transcriptome response to modification of incubation temperature.

Pathway analysis.
To functionally link miRNAs to the physiological effects of modifications of embryonic incubation temperatures, differentially expressed miRNAs and potential target genes were analysed using Ingenuity and its Knowledge Base software. To simplify and comprehend the resulting complex biological pathways and key genes, we aggregated significant terms and pathways into "major functional categories" based on shared terms and keywords of those pathways. Pathways in significant functional categories for each treatment condition at embryonic stages are given in Supplementary Tables S3 and S4.
In breast muscle of embryos, H10 showed activation of pathways related to cell maintenance and proliferation, organismal and tissue development, and nutrient metabolism, while L13 showed stimulated cell maintenance and proliferation, organismal and tissue development, genetic information and nucleic acid processing, cell signalling, and interaction and response to stimuli. Hind muscle of embryos showed deactivation of cell maintenance and proliferation, while L13 affected cell maintenance, proliferation, and differentiation pathways. Detailed information of pathway analyses can be found in Supplementary Tables S3 and S4. Functional miRNAs obtained from  Table 2. Functional miRNAs and potential target genes for selected treatment conditions. H10 and L13 that were highly associated with various target genes and that therefore were related to at least 4 out of 8 major categories are provided in Tables 3 and 4.
Regulatory networks predicting potential physiological effects. Representative miRNA-mRNA regulatory networks were modelled for H10 (Fig. 3) and L13 (Fig. 4) treatment in breast and hind muscles. The networks integrate miRNAs, potential target genes, and Ingenuity biofunctions and display an enrichment of miRNA-mRNA pairs related to major functional categories of cell maintenance, proliferation, and differentiation as well as tissue and organ development (Supplementary Tables S3 and S4).  Representative miRNA-mRNA regulatory networks demonstrate complex connectivity and relationships between the two molecular features. Sets of miRNAs target several genes that assemble into biological pathways and hence regulate these pathways. For example, miRNAs derived from H10 revealed activation of cytoskeletal organization and inhibition of cytoskeletal formation, demonstrating the overall fine tuning and balancing impact of miRNA posttranscriptional regulation (Fig. 3A,C). We also observed activation of pathways involved in white blood cell quantity, vasculogenesis, and thermogenesis as well as stimulation of body size. Pathways involved in reduced organismal death and perinatal death are shown in Fig. 3B,D.
Breast muscle that experienced reduced incubation temperature during ED10-13 showed stimulation of pathways related to proliferation, activity, formation, differentiation, and homeostasis of white blood cells. Regulated pathways included reduced growth of neurites (Fig. 4A). In hind muscle, we identified pathways involved in activation of apoptosis and cell survival (Fig. 4C). In breast muscle, we identified pathways related to white blood cell quality and development of body trunk. We also identified deactivation of organismal death and bone size (Fig. 4B) and, in hind muscle, inhibition of growth of connective tissue (Fig. 4D). Additional information for all pathways derived from miRNA-mRNA relationships indicated from integrated data analysis is available in Supplementary Figs S2-S5.

Discussion
Accumulating evidence suggests that modification of embryonic incubation temperature can result in phenotype variations of D35 chickens, such as adaptation to environmental conditions like heat stress. We have previously reported that changing incubation temperature during embryonic myogenesis influences weight gain and meat quality of broilers 13 . High temperature (38.8 °C) between ED 7 to 10 or ED 10 to 13 was found to be associated with higher body weights of D35 broilers compared to broilers from complementary normal (37.8 °C) or lower (36.8 °C) temperature groups 13 . Late increase of body weight at adult stage was shown to be due to increased size of breast muscle rather than hind muscle 3 .In fact, whereas breast muscle showed a 5.7% increased weight after high temperature compared to the control (351 g vs. 331 g) leg muscle showed only 3.9% increase (381 g vs. 367 g). This is in line with the fact that the treatments coincide with secondary fibre development and thus may have stronger impact on the breast muscle dominated by type II fibres originated from secondary fibres than on the hind muscle mainly consisting of type I fibres originated from earlier developing primary fibres. Moreover, lowered incubation temperature had in immediate effect on the in-ovo body weight that was reduced at the respective time points ED10 and ED13; higher incubation temperature slightly but non-significantly increased the embryo weight 17 . Also the D35 animals used in this study had higher body weight after transiently incubated at higher temperature between ED7-10 or ED10-13, whereas lower incubation temperature showed no effect on body weight later in life.
Further, we have previously demonstrated that both increasing and decreasing incubation temperature (1 °C from the control, 37.8 °C) immediately affects transcriptome profiles of embryonic muscle and associates with transcriptional changes of muscle of D35 chickens, indicating potential long-term transcriptomic effects of treatment during embryonic incubation temperature to adult stage 4 . This study now establishes posttranscriptional regulation by miRNAs in the above phenomenon.
Indeed, we found many differentially expressed miRNAs after thermal incubation treatments at the embryonic stage, compared to only a few differential miRNAs in D35 chickens. These results suggest that miRNAs are primarily, immediately affected and contribute to the regulation of gene expression at in-ovo stages in acute response to environmental circumstances during embryonic development, when thermoregulatory systems are not yet fully functional. In contrast there are only subtle changes of miRNAs in D35 chickens associated with embryonic thermal treatment. Mechanisms other than miRNA regulatory networks seem to be more important for long-term transcriptional changes in response to in-ovo thermal treatment such as epigenetic modifications which might have more persistent effects on mRNA expression.
For functional analysis, we focused on those thermal treatment conditions ((H10) and (L13) in embryos (L10) and (L13) in D35 chickens) that were previously found to show the most prominent mRNA-transcript changes 4,5 . Pairs of miRNAs and mRNAs were identified by in-silico target prediction complemented by experimental evidence based on the correlation analysis of miRNA and mRNA expression data obtained from identical samples. Considerable numbers of miRNAs and potential mRNA targets that were both regulated in the context of the experimental conditions were detected for H10 (high temperature during ED7-10) and L13 (low temperature during ED10-13) for in-ovo stages. Biological functions of the identified differential miRNAs at H10 and L13 in the regulatory contexts with their potential target genes were considered using Ingenuity analysis software and Knowledge Base. While most knowledge about the functional role of miRNAs comes from studies of pathological conditions, in particular cancer, our results provide evidence for miRNAs' role in ontogenetic proliferation and differentiation processes. In fact, H10 treatment consistently shifts expression of miRNAs related to these cellular developmental processes at in-ovo stages. Moreover, H10 treatment promotes pathways related to organismal survival and carbohydrate metabolism. In particular, thermogenesis, which is involved in thermoregulation, is important in breast muscle but not hind muscle. At the level of biofunctions, L13 treatment affects pathways related to cellular and organismal development via processes of proliferation, differentiation, and death. This is similar to H10 conditions; however, largely different miRNAs and target genes are shifted, indicating that alternative pathways are addressed to keep conditions close to homeostasis.
In breast muscle, two differential miRNAs, miR-138 (Supplementary Tables S3 and S4) and miR-3017a (Supplementary Tables S3), relate to genes that were enriched for major category 2 and predicted for activating thermogenesis at H10. Among those target genes was CSPG4, which is known to be involved in vasculogenesis. Interestingly, 12 miRNAs, including miR-133, miR-199, and miR-212, were associated with genes in major category 3 that are related to metabolism and synthesis of carbohydrates: upregulated ADRB2, CX3CL1, and PPP1R3B; and downregulated CHPF and CHST3. While PPP1R3B was found in liver and skeletal muscle tissues, it is also involved in regulating glycogen synthesis by forming a glycogen-targeting subunit for phosphatase PP1 18 .
In hind muscle at H10 conditions, differentially expressed miRNAs, including miR-199a-5p, miR-1915, and miR-638, were related to 14 genes in major category 2, which is associated with accumulated body size and reduced perinatal death, including CUL4B, ITSN1, MLL5, and MYH11. Interestingly there is evidence that miR-1915 is involved in apoptosis via interaction with the antiapoptotic protein Bcl-2 and that miR-638 relates to vascularization 19,20 . In fact, pathways analysis here indicates corresponding links to perinatal death and cellular development (Fig. 3). Further, the results show that hydrolysis of carbohydrates mediated by genes such as MTM1, NT5E, PLCB1, and MGLL, which are related to nutrition metabolism in major category 3, are targeted by miRNAs shifted at H10 in hind muscle (miR-199a-5p, miR-212, and miR-222).
Interestingly, both tissues of the H10 group showed stimulation of major category 1, organization of cytoskeleton and cytoplasm, which is linked to several miRNA-targeted genes, including CDKN1B, MTUS1, PIK3R1, PLXNB2, and SYK. Especially PIK3R1, which is represented in multiple functional categories, was predicted to be a potential target of miR-739 in breast muscle and miR-2861, miR-3960, and miR-4592 in hind muscle. In addition, several target genes were associated with improved survival in major category 2 by deactivating organismal death and perinatal death in both breast and hind muscle.
Among the differentially expressed miRNAs miR-133 is prominent as it is known as a muscle-specific miRNA also called "myomiR," affecting muscle proliferation, myotube formation, and differentiation 21,22 . Chen et al. (2006) showed that upregulation of miR-133 was associated with myoblast proliferation but reduced cell differentiation 23 . In chickens, miR-133a was shown as stimulatory factor in late-stage development in response to myogenin 24 . In this study miR-133 was upregulated in breast muscle but not in hind muscle. This is in line with the differential expression of miR-133 found in porcine M. longissimus dorsi and M. psoas major, i.e. white and red muscles with different proportions of white (type II) and red (type I) muscle fibres like breast and hind muscle 20 . Moreover, differential expression of miR-133 in the two muscles may also be due to different degree of proliferation and differentiation processes in the muscles reflecting the subsequent development of white fibres (FTG) in breast muscle and red fibres (STO) in hind muscle.
MiR-212 revealed pleiotropic properties in various tissues so far mostly in relation to pathological states like cancer and inflammation and mostly related to neuronal cells 25 . We previously found miR-212 among co-expressed miRNAs that affect muscle properties related to meat quality in pigs 26 . Here miR-212 showed increased abundance in breast and hind muscle at H10 where it was assigned functional roles in several pathways related to cellular and organismal growth (Fig. 3).
Overall, H10 treatment effects tended to promote high body and muscle weight compared to low temperature treatment, which is in line with previous observations 17 . Accordingly, shifts of miRNA expression mainly affect pathways related to cell survival, organization of cytoskeleton and cytoplasma, angiogenesis and vascularization and also more specific pathways of myogenesis. L13: Low temperature during ED10-13. Compared to H10 treatment, L13 had less differentially expressed genes. Low temperature treatment during ED10-13 (i.e., formation of secondary fibres during myogenesis) activated several miRNAs that in turn regulate major categories 1 and 2. In breast muscle, major category 1 showed potential activation of 23 biological functions, such as cellular activity of formation and engulfment. Formation of cells was influenced by upregulation of RNF2 and TCF3 and downregulation of RUNX2. These genes were targeted by miR-130c, miR-263a-star, and miR-312-5p. Moreover, other differentially expressed genes, including JAM3, PATZ1, PICK1, SOAT1, and SRF, were also associated with these miRNAs.
It is interesting that differential miRNAs from L13 treatment related to the predicted tendency of reduced bone size. Target genes NOTCH2, HIVEP3, AMER1, and RUNX2 were regulated by downregulation of let-7, miR-93, and miR-130c. L13 treatment is associated with low body weight of embryos compared to the high temperature 17 . Muscle size is a major force inducing adaptation of bone size. It is speculated that reduced bone size could be well correlated and of consequence of low body weight.
Furthermore in breast muscle at L13 eight differentially expressed genes and 16 differential miRNAs, including let-7, miR-92, and miR-93, belonged to major category 3, with biofunctions of uptake of carbohydrate and D-glucose.
Let-7 has been a major topic of discussion for functional roles of miRNAs. miRNAs of the let-7 family play major roles in regulating proliferation during ontogenesis and particularly during myogenesis 22 . Let-7 family members are associated with aging in humans and downregulation of cell cycle control, such as cellular proliferation and differentiation pathways 27 . Moreover, during myogenesis, let-7 can suppress Dicer and HMGA2, which have roles in adipogenesis and mesenchymal differentiation 28 . Interestingly, long non-coding RNA H19, possessing multiple let-7 binding sites, is proposed to prohibit let-7 from binding to other targets 29 . Here let-7 showed different abundance due to L13 treatment in breast muscle and was assigned to various pathways related to cellular, organismal and metabolic processes. Let-7 was not found differentially expressed in hind muscle at L13. miR-93 can potentially downregulate AKT3, which reduces proliferation and facilitates differentiation of myoblasts in skeletal muscle development 30 .
We reported that L13 treatment is associated with low body weight of embryos compared to high temperature treatment 17 . Bone and muscle are highly associated during developmental stages 31 , so reducing bone size could also correlate with low body weight in-ovo. miRNAs shifted at H10 and L13. Due to their functional annotation to several biofunctions and functional networks and significant response to both treatments and/or in both muscles a number of miRNAs, including members of the miR-199 family, miR-222, miR-460, and miR-5109, play important roles in fine tuning gene expression at the background of various mRNA expression profiles given in breast and hind muscle at day 10 and day 13 of in-ovo development.
Members of the miR-199 family are involved in multiple roles, including stem cell differentiation and embryo development 32 . As reviewed by Gu and Chan (2012), miR-199 family is differentially expressed in various tissues and pathological states as well as developmental processed, providing an example of the versatile function and regulation of miRNAs 32 . Accordingly, in our study miR-199 family members showed different abundances in breast and hind muscle at H10 and in breast muscle at L13 and were assigned to developmental processes at the cellular and organismal levels but also to metabolic processes as shown in Figs 3 and 4 and supplemental materials.
Due to its impact on cell proliferation and apoptosis in cancer, functions in physiological and pathological processes in the cardiovascular system and responsiveness to environmental exposure miR-222 has been proposed as a biomarker for cancer, cardiovascular diseases and environmental stressors [33][34][35] . Here we found regulation of miR-222 breast and hind muscle at H10 and L13 conditions, thus in samples where proliferative process are important, i.e. at early in-ovo development, and that were exposed to environmental stressor, albeit physical stress not chemical stressor as in studies reviewed by Vrijems et al. (2015). Interestingly, miR-222 mediates atrophy in denervated fast muscles in rates, i.e. processes that involve rearrangement of fast and slow muscle fibres as during myogenesis that is addressed in our study 36 .
For miR-5109 and miR-460 that were both differentially expressed in breast muscle at H10 and L13 knowledge of their function is scarce. MiR-460 showed lowered abundance in breast muscle at H10 and L13. Interestingly, miR-460 was detected as a growth related miRNA in tilapia 37 . Our study of functional mi-RNAs with their linked mRNAs revealed that miR-460 is involved in many structural and metabolic cellular processes but not directly linked body growth. MiR-5109 showed increased abundance at H10 but decrease at L13. There is no further information about miR-5109 available. The divergent regulation at H10 and L13 might be due to the direction of temperature change or due to different developmental stages that both are associated with specific mRNA expression profiles.
For miR-199, miR-222, miR-460 we found the same direction of regulation in breast and hind muscle and/or at H10 or L13 conditions. This reflects the complexity of miRNA expression and its role for regulation of mRNA transcript abundance; where due to the numerous potential target genes of each single miRNA changes of their frequency contribute to the maintenance of homeostasis at the background of different mRNA profiles. MiRNAs can be assigned to functional pathways via their link to target genes. Therefore, in order to get an comprehensive knowledge, global studies of expression profiles at various biology-based systems are required 32 . Our study contributes to this collection of knowledge and particularly provides in-vivo data from a non-pathological organism at different developmental stages.

Summary.
In the present study, we have demonstrated that modification of embryonic incubation temperature immediately affects miRNA expression profiles of breast and hind muscles of chicken embryos and is associated with altered expression of miRNAs in D35 chickens. An integration analysis of miRNA data and previous matched-sample mRNA data revealed functional miRNAs and enabled assembly of miRNA-mRNA regulatory networks related to biological pathways and potential physiological effects, though the functionality and targeting of each single miRNA and corresponding mRNAs remains to be validate by adequate in-vitro studies Differentially expressed miRNAs and targeted mRNAs showed treatment condition specificity in terms of timing of treatment (ED7-10 or ED10-13), tissue type, and stage of development. The large repertoire of miRNA-mRNA pairs that are shifted in various experimental groups but that finally fine-tune similar biofunctions reflects a considerable functional biodiversity and resilience.
This study reveals substantial immediate alterations of miRNAs due to experimental environmental conditions, whereas long-term miRNA responses were minor. This indicates a major regulatory role of miRNAs in acute responses to modified environmental conditions. The study adds to the cumulating knowledge of the function of miRNAs and particularly provides data of healthy organisms at different ontogenetic stages.