In some European Mediterranean countries, suckling lamb meat is highly appreciated because of its mild flavor and soft texture (1, 2). In Spain, dairy sheep farms produce both milk and suckling lambs. Lambs produced in this way are characterized by slaughter weights of 9–14 kg and are exclusively milk-fed animals (3). Suckling lamb carcass quality has traditionally been mainly based on weight, fatness, meat and fat color, conformation and carcass composition (4). Several factors, such as breed (5, 6), sex (7, 8), slaughter body weight fatness level (9), and diet, have been reported to influence the quality of lamb carcasses.

Assaf is a semi-fat-tail breed that is a stabilized cross of the East Friesian breed and the fat-tail Awassi breed from Israel (10). Since the introduction of this breed in Spain in 1977 for milk production, the Assaf breed has largely replaced local breeds due to its high milk production potential (11, 12). Regarding lamb meat production, Spanish Assaf lamb carcasses have been shown to be shorter than carcasses of similar weight from other Spanish native breeds, such as Churra and Castellana (13, 14). In relation to fatness traits, and despite the higher fat tail accumulation of Assaf lambs, overall carcass fattening is similar between Manchega and Assaf, although carcasses of Churra and Latxa suckling lambs seem to be fatter (12, 14).

In suckling lamb carcasses, the content of internal adipose deposits, such as perirenal and pelvic fat, also known as kidney knob and channel fat (KKCF), represent a notable proportion of the total fat of lambs during the first weeks of life (15) and are vital for the neonatal animal to maintain its body temperature (16). These depots are used as classical indicators of total carcass fatness (17). Among different carcass parameters, the KKCF percentages in the carcass, either visually assessed or weighed, have been categorized as the best predictors of suckling lamb carcass quality in Spanish breeds (4, 18). The tissue composition of the leg, loin, or shoulder can also be used for carcass quality prediction; however, all abovementioned traits are difficult to measure during routine production (4). Thus, at the commercial level in Spain, it is common practice to classify carcasses in which the kidneys are completely covered by fat as high-quality suckling lamb carcasses because there is a positive correlation between perirenal fat and global carcass fatness [Commission Regulation (EC) No 2107/1999]. In addition, due to its ability to produce adipokines, perirenal fat is considered not only a deposition tissue but also an important metabolically active tissue (19, 20).

Studies based on RNA-seq technology provide a view of the extent and complexity of eukaryotic transcriptomes, which are essential for revealing the molecular constituents of cells and tissues and for understanding development and disease (21). RNA-seq also provides quantitative measures of gene expression levels and allows the discovery of novel transcribed regions (22). In livestock, the study of transcriptomes has been used to identify the expression of genes in several tissues related to productive traits, such as the mammary gland (23) or intramuscular fat (24) in dairy cattle, adipose tissue in pigs (25, 26), and fat depots in suckling lambs (27).

In relation to fatness traits associated with carcass quality, some RNA-seq studies have analyzed the transcriptome profiles associated with different fat depot deposition patterns by comparing different breeds of the same livestock species, such as cattle (28) or pigs (29). In sheep, the majority of RNA-Seq studies aimed at understanding fat metabolism have focused on the fat tail transcriptome in adult sheep (30–35). In suckling lambs, considering the interest in internal adipose deposits as indicators of carcass fatness (17) and carcass quality (4), our group has previously compared the perirenal fat transcriptome of Churra and Assaf suckling lambs, which show different carcass fatness levels, quantified based on the KKCF indicator trait (27). However, no research has yet been performed to compare suckling lambs of the same breed with different carcass KKCF contents.

In the interest of gaining a deeper understanding of the functional genomic mechanisms associated with carcass quality in suckling lambs, this study aimed to reveal the metabolic processes underlying the control of fatness levels in Assaf suckling lambs. To this end, we performed a comparative analysis between RNA-seq datasets from the perirenal fat of Spanish Assaf suckling lambs showing different KKCF carcass contents. The results presented here provide additional information to better understand the processes underlying fat deposition in sheep. Further studies must confirm whether these processes are common to other ovine breeds or whether breed-specific mechanisms may influence the fatness levels of suckling lamb carcasses.

Some carcass parameters, such as kidney knob and channel fatness scores or loin and leg tissue compositions, have been proven to be useful predictors of the composition and quality of suckling lamb carcasses (4). However, the molecular mechanisms of fat metabolism are complex, and characterizing the basis of fat deposition in livestock species may be relevant for maximizing the advantages that can be achieved through genetic selection. In general, previous studies have investigated the transcriptomic profiles of fat tissue between different sheep breeds (27, 33, 35, 49). However, there are no available studies comparing suckling lambs of the same breed with differences in fat deposition.

Research on different mammalian species has revealed dynamic features of fat tissue. Depending on the cellular lineage from which they originate, there are at least three different categories of adipocytes: brown, white, and brite or beige adipocytes (50, 51). In both large and small mammals, brown adipose tissue (BAT) depots are very important in the early stage after birth (52, 53). In fetal sheep, beginning in mid-gestation, multilocular adipocytes, characteristic of BAT tissue, are clearly visible (50). Then, close to term, the fat deposits increase in size, and a mixture of brown and white adipocytes is observed (54). After this time point, coinciding with the time of birth, brown adipocytes gradually disappear until all fat deposits consist of white adipose tissue (WAT). Beige or brite adipocytes, found in beige adipose tissue (BeAT), are very similar to brown adipocytes; however, they come from different cell lines, sharing common cell precursors with white adipocytes (55, 56). Focusing on the BAT to WAT transition, it is not clear how long the period of gradual BAT disappearance lasts. Basse et al. (57) reported that this transition occurs before the first 15 days of life. In contrast, based on comparing the perirenal fat of Churra and Spanish Assaf suckling lambs, Suárez-Vega et al. (27) reported that BAT was still present in some animals at approximately 20 days of age, mainly in Assaf lambs, which were younger at slaughter. These observations agree with the findings of Clarke et al. (58), who reported a large increase in lipid deposition from 7- to 30-day-old ewe lambs, associated with a decline in the amount and activity of UCP to basal values in the ewe lambs. These authors suggested that the precise time scale of this transition process could be regulated in part by the lamb’s body temperature, which determines whether adipose tissue is required for heat production (i.e., BAT) or as an endogenous energy source (i.e., WAT). In this context, considering the expression levels of some of the most representative markers of BAT, BeAT and WAT (according to various authors) in the KKCF suckling lamb samples analyzed here shown in Figure 2, we can see that although WAT markers show the highest expression levels, the expression of some BAT gene markers, such as UCP1 and CIDEA, can still be identified. Hence, to properly interpret the results of further analyses, it is important to consider that the perirenal fat samples of suckling lambs analyzed in this work may still be within the BAT to WAT transition period.

As a first step in the interpretation of the transcriptome analysis reported here, we considered it important to determine which genes showed the highest proportions of accumulated gene expression in the fat transcriptome of suckling lambs. The identified CRGs were associated with energy production in the form of ATP in the mitochondrion based on the enrichment analysis performed herein (see the complete list in Supplementary Table S4). This organelle also appeared in some terms in the list (inner mitochondrial membrane protein complex, respiratory chain complex, mitochondrial respirasome, etc.). These results hint at the high metabolic activity occurring in the adipose tissue analyzed in all analyzed samples and agree with the central role that mitochondria play in the metabolism of adipose tissue (59). Moreover, one of the BAT-specific functions is transforming mitochondrial energy into heat in adaptive thermogenesis (59). Our previous study already highlighted this result by comparing the perirenal fat transcriptome of Spanish Assaf and Churra breed suckling lambs (27). From the total of 670 genes identified in that study as CRGs for these two breeds (expression levels >150 FPKM in both breeds), a total of 318 genes (~ 47.5%) are included in the list of Spanish Assaf suckling CRGs reported in this study. The fact that almost half of the genes expressed in the perirenal fat tissue in both studies are common suggests that there is consistency among the fat transcriptomes of suckling lambs, regardless of their breed.

Among the CRGs, focusing on the top 10 most highly expressed genes identified in the two groups of KKCF tissue samples analyzed in this work (Table 1), it was notable that the most highly expressed gene in both the High- and Low-KKCF groups was that encoding fatty-acid-binding protein 4 (FABP4; Table 1; Figure 2). FABP4 is one of the most abundant proteins in mature adipocytes and is considered a genetic marker of WAT (60). However, whether its localization is limited to mature adipocytes is not clear. The expression of FABP4 has been reported to increase with the progression of adipogenesis in fat tissue beneath the lower dermis in mice (61). FABP4 physically interacts with PPARG (62), which controls the expression of genes responsible for the differentiation of preadipocytes into mature white adipocytes and is considered a key regulator of adipogenesis (63). The expression level of PPARG, which is also considered a WAT marker, in the samples analyzed here was also >150 FPKM (Figure 2). Because FABP4 negatively regulates PPARG at the posttranscriptional level and represses preadipocyte differentiation (64), lower expression of FABP4 has been associated with obesity in different studies (65, 66). In our samples, the expression levels of these two WAT markers, FABP4 and PPARG, were lower in the Low-KKCF group than in the High-KKCF group (see Figure 2), although none of them were identified as DEGs.

Considering the remaining top 10 highly expressed genes in the two groups of samples (Table 1), it seems worth highlighting that many of the genes common to the two lists are related to ribosomes (RPLP1, RPLP2, RPL37, and RPS2), most of which are constituents of this organelle. The synthesis of ribosomes and translation are essential for the growth, differentiation, proliferation and development of animal cells (67). A study using human adipose-derived stromal cells (hASCs) reported an overall reduction in protein synthesis after adipogenic induction and a downregulation of the expression and translational efficiency of ribosomal proteins (68). Based on this, the high expression levels of ribosome-related genes identified in both groups of samples under study (High- and Low-KKCF) may indicate that extensive cellular remodeling is taking place in the studied adipose tissue because of the early stage of growth of the suckling lambs analyzed, independent of the proportion of carcass KKCF.

The expression levels of two other important fat tissue gene markers, the UCP1 and CIDEA genes, are also represented in Figure 2. The UCP1 gene, which is a classical gene marker of BAT, showed very heterogeneous expression levels among samples, and even within groups of samples, there were important differences (209.06 ± 245.06 FPKM). The expression of this gene reached the 150 FPKM threshold only in the Low-KKCF group and not in the High-KKCF group. The identification of gene expression of this BAT marker in some of our suckling lamb samples contrasted with the results reported by Basse et al. (57), who did not observe UCP1 expression in four-day-old lambs. However, our results agree with those of Suárez-Vega et al. (27) and Yuan et al. (69), who reported expression levels in the perirenal fat adipose tissue of sheep at 17–35 days of age and sheep of 2–12 months of age, respectively. The previously indicated higher expression levels observed for UCP1 in the Low-KKCF group compared with the High-KKCF group (Figure 2) can be explained by the fact that the Low-KKCF group included animals that were slaughtered at a younger age and could suggest that more thermogenesis was taking place in the adipocytes of the Low-KKCF animals, as UCP1 represents the key point of this pathway. Indeed, this gene was one of the Low-KKCF DEGs identified in the DE analyses and will be discussed later. On the other hand, the CIDEA (cell death-inducing DNA fragmentation factor α-like effector A) gene, which was suggested by Garcia et al. (70) for use as a BeAT marker, as shown in Figure 2, showed consistently high expression levels among samples (226.89 ± 86.73 FPKM) and was included in the list of CRGs reported in this work. This gene is a marker of thermogenesis in adipose tissue and plays an important role in the regulation of this process by inhibiting UCP1 activity (71, 72). In addition, the CIDEA protein colocalizes around lipid droplets (73). Human CIDEA inhibits lipolysis, causing a significant decrease in glycerol release (74). Because CIDEA expression occurs specifically in BAT and BeAT, this is an interesting marker gene for performing further studies on its role in the development of adipose tissue.

The experimental design of this study gave us the opportunity to identify a set of 80 genes whose expression levels may influence fat deposition levels in perirenal and visceral depots in Spanish Assaf suckling lambs. We acknowledge that the number of samples analyzed in the present study was limited, although it was higher than the minimum number of biological replicates recommended by Conesa et al. (75) for RNA-Seq studies and similar to the number of samples analyzed in other recent studies evaluating the fat transcriptome in sheep (30, 35, 76). However, an additional limitation to our experimental design in this specific study was the significant differences in age in the two groups of animals initially considered based on their different levels of perirenal fat accumulation. Although these differences supported the general observation that the lambs slaughtered at an older age show more fat accumulation than younger lambs, individual deviations from this general observation are of special interest when studying the “fatness” phenotype. Hence, the presence of the two animals in the study that did not conform to the fatness-age correlation (one animal was Low-KKCF and Old, and one other animal was High-KKCF and Young) led us to include the “age” category in our model. Although the contribution of the inclusion of this second category in our DE model may appear limited due its nonuniform distribution in the two considered fatness groups (with only 2 animals that did not show the fatness-age correlation), we consider these animals to have provided valuable information that helped us to identify a set of DEGs and enriched terms that were specifically associated with fatness levels in the KKCF depots of suckling lambs, independent of the time needed to reach the slaughter weight defined by market criteria.

In relation to the DEGs identified in the two groups defined based on fatness accumulation in KKCF depots in this study, the results of the enrichment analyses (graphically summarized in Figures 3, 4) highlight relevant metabolic and developmental processes that clearly differentiate the two groups. Considering the LowKKCF-DEGs, the significant enriched terms were globally summarized as being related to fatty acid metabolic processes and other metabolic pathways associated with obtaining ADP and ATP (“purine nucleoside and ribonucleoside biphosphate metabolic process” and “ribonucleoside and nucleoside biphosphate and phosphate metabolic process”), whereas for the High-KKCF group, the enriched terms were globally grouped with terms related to vessel endothelial development (regulation of blood vessel cell endothelial cell proliferation, fibroblast migration and regulation) and cellular response mechanisms toward acid chemicals (e.g., “response to amino acid,” “cellular response to amino acids” among others). Below, we present a discussion of these results, attempting to explain their biological meaning in relation to the analysis performed.

Focusing first on the Low-KKCF-related results, the genes supporting most of the enriched terms identified in this group were directly related to lipid metabolism and different metabolic processes related to nucleoside and ribonucleoside (especially purine) bisphosphates (Supplementary Table S6). Some of these genes were ACSS2, ACLY, ELOVL6, DGAT2, FASN and GPAM. Among the lipid metabolism-related terms, we found a large variety of processes, such as “fatty acid biosynthetic process” and “cholesterol biosynthetic process.” These significantly enriched pathways suggest that the synthesis of many types of lipids occurs in the adipose tissue of Low-KKCF animals. For instance, the FASN gene encodes fatty acid synthase, which is a central enzyme in de novo lipogenesis. FASN is an enzyme that catalyzes a rate-limiting step in de novo fatty acid synthesis in ruminants, as its main function is to catalyze the synthesis of long-chain saturated fatty acids (77). In addition, the increased expression in the Low-KKCF group of genes such as DGAT2 and GPAM, which encode key enzymes for the synthesis of triacylglycerols (78, 79), also suggests that in the perirenal depots of these animals, the biosynthesis of these compounds also takes place. The mentioned synthesis processes are linked with ELOVL6, a gene that supports a remarkably large number of significantly enriched terms identified in the Low-KKCF group (see Supplementary Table S6). This gene is considered a key gene involved in regulating the balance of the fatty acid composition, as it is involved in the elongation of long-chain saturated and monounsaturated fatty acids (80–82). The identification of ELOVL6 as a LowKKCF-DEG highlights the potential importance of remodeling processes, in addition to fatty acid synthesis, in the perirenal adipose tissue of these animals compared with the High-KKCF group.

Considering these results and a detailed review of gene expression and transcriptomic studies focusing on fat deposition in livestock species, such as those comparing animals with higher and lower fatness levels, we can see that the expression levels of all the mentioned genes are generally associated with the groups of animals showing higher fat deposition levels. This has been shown in relation to different fat depots in pigs (29, 83), calves (84), and sheep (35) and suggests that these genes are general genes that are components of the universal genetic basis of fat deposition. Hence, we acknowledge that the fact that all the mentioned genes related to biosynthesis and lipid remodeling showed higher expression levels in the animals with lower KKCF values in our study appears to be a discrepancy. However, as previously stated, for the interpretation of our results, it is important to consider the critical growth stage of the suckling lambs analyzed in the present study, which were slaughtered at a very early age (17–27 days for the Low-KKCF group), resulting in a completely different situation than that of the animals analyzed in the previously mentioned studies on livestock fat deposition, in which fat tissue development had already been completed. Indeed, as previously noted, the lambs analyzed here may still be within the BAT to WAT transition period (based on Figure 2).

In this regard, it is interesting that the role of ELOVL6 in the regulation of fatty acid length is important for BAT and BeAT function (85). In mouse adipose tissue, the loss of the ELOVL6 gene causes reduced expression of components of the mitochondrial electron transport chain, which leads to a reduction in the thermogenic capacity (86). Indeed, this gene, together with UCP1, SCD, LEP and PPARGC1B, contributed to the enrichment of terms such as “positive regulation of cold-induced thermogenesis,” “cold-induced thermogenesis,” and “adaptive thermogenesis” among the LowKKCF-DEGs (Supplementary Table S6). The identification of these enriched terms may help to explain the previously mentioned discrepancy of higher lipid biosynthesis in the Low-KKCF group compared with the High-KKCF group because at least some of the fatty molecules generated in these biosynthesis processes could be used for such thermogenic activity. The increased expression of the UCP1 gene in the Low-KKCF group also contributed to the enrichment of the “mitochondrial inner membrane” term, which is directly related to UCP1 because this protein is inserted in this membrane. Additionally, UCP1 was involved in the enrichment of terms such as “fat cell differentiation,” “brown fat cell differentiation,” “response to dietary excess,” “response to nutrient levels,” and “response to fatty acid levels,” among others (Supplementary Table S6). The enrichment of these terms was supported by other genes, including LEP or DGAT2 in many cases.

The LEP gene encodes leptin, an adipokine produced by adipose tissue as a feedback mechanism that leads to appetite reduction and has been shown to inhibit glucose conversion to lipids and fatty acid incorporation by porcine adipocytes (87). In livestock studies focused on adult animals higher expression levels of this gene are normally associated with animals with higher fat deposition levels (29, 88). However, in young animals, leptin has been shown to affect growth, and periods of rapid growth require high levels of LEP (89–91). Hence, different studies have suggested that leptin acts as a skeletal growth factor and an inducer of longitudinal growth (92, 93). All of these findings and inferences are in accord with the identification of increased expression of LEP in Low-KKCF animals, which in general reached the commercial slaughter weight of 11–12 kg at a younger age. Additionally, our results agree with those of Xing et al. (94), who compared the Landrace breed of pigs, characterized by high growth rates and lean carcass percentages, with the Songliao breed, with higher fat deposition. In this case, the Landrace breed showed the highest levels of LEP expression.

Finally, in relation to the previously mentioned increase in UCP1 expression observed in the Low-KKCF group, we should take into account that while free fatty acids activate UCP1, purine nucleotides, such as ATP, strongly inhibit the function of UCP1 (95). In this context, we also found enriched terms for the LowKKCF-DEGs associated with purine nucleotide/ribonucleotide metabolic processes (Figure 3A), supported by genes such as ACSS and ACLY, which are both related to acetyl-CoA metabolism (96). Hence, although our results suggested that the Low-KKCF animals were at an earlier stage of adipose tissue development than the High-KKCF animals, where BAT-related processes such as thermogenesis were still observed, the potential inhibition of the UCP1 function identified in the Low-KKCF group suggested that they may have been in the last stages of the BAT to WAT transition. This was also supported by the overall higher expression of WAT gene markers than BAT gene markers shown for both groups in Figure 2.

Focusing now on the enrichment results obtained for the HighKKCF-DEGs (Figure 4; Supplementary Table S7), we found several terms related to blood vessel endothelial cell proliferation connected to terms associated with the cellular response to acid chemicals through a central gene, VEGFA. This gene is a member of the PDGF/VEGF growth factor family and encodes vascular endothelial growth factor A (VEGFA), which induces the proliferation and migration of vascular endothelial cells and is essential for both physiological and pathological angiogenesis (97). VEFG-A is also a crucial molecule controlling vasculogenesis, vascular permeability and tissue remodeling (98). Angiogenesis, the process of forming new blood vessels, is crucial to adipogenesis (99), and VEGF-A is responsible for most of the proangiogenic activity in fat tissue (100, 101). Since the vasculature is essential for adipose tissue expansion, as it supplies oxygen and nutrients and is also responsible for the transport of fatty acids to other tissues, angiogenesis is considered to be a rate-limiting step in fat tissue expansion (102). VEGF-A levels are regulated by both external factors, such as exercise (103) and dietary protein levels (104), and internal factors, such as hypoxia and various growth factors (105, 106). Adipocytes store lipids in lipid droplets when nutrients are abundant, and as they increase in size, oxygen becomes less available, which leads to a hypoxic status. In this context, hypoxia has been described as the most important mechanism for VEGF modulation (107). This would explain the connection of VEFG with the enriched terms “response to acid chemical” and “cellular response to acid chemical.” All of these findings and inferences agree with the observed increase in VEGF-A gene expression levels in the High-KKCF group compared with the Low-KKCF animals, as their larger fat depots could have reached the condition of moderate hypoxia with the subsequent activation of the angiogenesis process in response (108). In relation to angiogenesis, VEGF-A is connected to terms related to fibroblast migration through the THBS1 gene. Thrombospondin 1 (THBS1) is a large adhesive extracellular matrix glycoprotein (109) that is expressed predominantly in visceral adipose tissues, and its expression has been found to be elevated in obese humans (110). According to Kong et al. (111), THBS1 overexpression promotes adipose tissue development and weight gain in mice. This biological role agrees with our results, as THBS1 was one of the HighKKCF-DEGs identified in our analysis. Interestingly, neither VEGF-A nor THBS1 have been reported as DEGs in general studies comparing adult animals with high and low fatness levels. This suggests that the vasculogenesis process associated with these genes may be a critical stage of fat tissue development in young animals, such as the suckling lambs analyzed here.

Other enriched terms found for the HighKKCF-DEGs related to acid chemical responses were “cellular response to amino acid stimulus,” “response to leucine,” “response to leucine starvation,” “cellular response to leucine” and “cellular response to leucine starvation” (Supplementary Table S7). The most relevant genes supporting these terms, SENS1 and SESN2, encode sestrins. Sestrins are a family of stress-inducible proteins that regulate metabolism based on sensing nutrient levels and redox status in cells that participate in the negative regulation of mammalian target of rapamycin complex 1 (mTORC1). Several terms related to mTOR were also enriched in the analysis of the HighKKCF-DEGs (“TOR signaling,” “negative regulation of TORC1 signaling,” “regulation of TORC1 signaling,” “negative regulation of TOR signaling”; Supplementary Table S7). mTOR is a complex that also plays important roles in adipocyte development and functional maintenance (112), for which a positive role in adipogenesis has been classically described (113). In the High-KKCF animals, increased levels of sestrin expression, leading to negative regulation of mTOR signaling (as suggested by the identified enriched terms previously mentioned), would be associated with the inhibition of adipogenesis. This would indicate that the developmental stage of the perirenal fat of these animals had moved beyond the preadipocyte differentiation phase toward an expansion phase, for which vasculogenesis is essential. In turn, the lower expression of sesrins in the Low-KKCF group would suggest that adipogenesis is activated in this group.

Overall, the transcriptomic analysis of perirenal fat in suckling lambs presented herein suggests that a higher content of KKCF in the carcass of suckling lambs is associated with a more mature stage of perirenal adipose tissue, associated with increased vasculogenesis and expansion of the tissue, where thermogenesis and adipogenesis are inhibited, compared with the Low-KKCF group. We have identified genes such as VEGF-A, THBS, SESN1 and SESN2 as being associated with this more advanced developmental stage of perirenal fat depots with higher KKCF values, which are associated with better carcass quality. While we acknowledge that the limited number of samples analyzed means that these results should be confirmed by future studies involving large population sizes, the findings presented here represent a valuable first attempt to understand the molecular mechanisms underlying fat deposition in suckling lambs. In addition, the genes indicated above are proposed as candidate genes, and their polymorphisms should be evaluated for potential direct associations with carcass fatness traits in suckling lambs.

The data presented in the study are deposited in the BioStudies repository, accession number E-MTAB-11839, and can be found at: https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-11839.

Ethical review and approval was not required for the animal study because the lambs included in this experiment were not subjected to any experiment, and their management was carried out following the usual management practices on farms raising suckling lambs with artificial lactation. The management, transport, and slaughter of the animals at a local slaughterhouse were in accordance with Spanish and EU legislation [Spanish Laws 32/2007, 6/2013, and RD 37/2014; Council Regulation (EC) 199/2009].

MA-G: formal analysis and writing—original draft. AS-V: conceptualization, validation, and writing—review and editing. PF, RP, and HM: investigation. JM: methodology and investigation. J-JA: conceptualization and investigation. BG-G: conceptualization, methodology, and writing—review and editing. All authors contributed to the article and approved the submitted version.

The research described here has been funded by the project EpiMilksheep (RTI2018-093535-B-100) funded by the Spanish Ministry of Science and Innovation. MA-G is funded by a predoctoral fellowship from the Junta de Castilla and León Government and the European Social Fund.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.