Interplay Between Stress and Reproduction: Novel Epigenetic Markers in Response to Shearing Patterns in Australian Merino Sheep (Ovis aries)

Introduction

In the past decade, there has been an increase in scientific reporting into the effects of adversity in early life on the participant's DNA profile within both human and animal studies (1, 2). This emerging research area continues to be driven by scientists globally to better understand a wide range of effects caused by a variety of intrinsic and extrinsic influences on the DNA profile of ongoing generations within observed genotypes. From within the cell of the early stage developing embryo, transcriptional and epigenetic changes to the cell are occurring via remodeling and reprogramming within the cell nucleus. What is still unclear is if/how external factors such as management practices could alter DNA during critical life history phases such as reproduction. This science is what is known as epigenetics (3). Generally, epigenetics represents the genome-wide study of the distribution of methylated and unmethylated nucleoside residues within the genome (4). Daxinger and Whitelaw (5) concluded that epigenetics refers to effects on phenotype (or on patterns of gene expression) that are passed from one generation to the next by molecules in the germ cells and that cannot be explained by Mendelian genetics (or by changes to the primary DNA sequence). In production animal research, recently, the need of epigenetic evaluation of maternal-fetal nexus especially in relation to environmental factors (e.g., climate), nutrition, and post-natal development and growth of progeny has been highlighted (6). Recently, researchers have employed DNA methylation tools to evaluate stress responses in mammals (7).

Early biomedical studies into animals' epigenetics have been focused on mice due to their ability to reproduce quickly and for researchers to gain fast results with multiple offspring from the same female. Because of the nature of sheep growth and time to reach puberty, a predominant single offspring, and the lack of funding for epigenetic research into the sheep, there have been very limited research studies in this field [see (8, 9)]. However, because of the foresight by early researchers, there have been substantial advances in current genomic technologies to allow for development of genome analysis and sequencing in livestock (10). To our knowledge, there is no peer-reviewed scientific literature of the epigenetic effect on the female sheep' DNA caused by a basic management intervention such as shearing. Our earlier research showed that sheep embryos from the same sire and dam that were placed into surrogate ewes and raised under the same environmental and management regimes had significantly varied phenotype (wool quality) (11, 12). Maternal stress at the early stages of pregnancy appears to result in greater implication of the epigenetic changes than that trigged post-parturition (13). The placenta acts as a connection between the mother and the developing fetus and stress, activates the maternal hypothalamo-pituitary adrenal (HPA) axis, and triggers stress hormone or glucocorticoid (GC) synthesis that reaches the fetus by transplacental passage (14). Prenatal stress and prenatal exposure to GCs have been shown to have long-term effects on the expression of genes associated with HPA function, neurologic function, and phenotype (15). In sheep, latest research has demonstrated that pharmacologically elevated cortisol can also result in negative effects on the animal's body condition and phenotype (wool fiber diameter) (16).

Shearing is a well-validated management intervention and known to generate acute stress responses in sheep (17). Benefits of early or even mid-pregnancy shearing are well known across various breeds of sheep such as increase in shelter seeking and cleaner udder areas with lesser wool, resulting in lambs born heavier and with improved survival chances (18). However, there have been no previous studies conducted in sheep on the influence of shearing pattern (once a year vs. twice a year) on maternal stress physiology, grazing behavior, and epigenetic effects on lamb phenotype. Our interest was to non-invasively monitor the ewe's HPA axis activity between shearing treatments during gestation; therefore, we employed a well-established wool cortisol assay (19, 20). We hypothesized that shearing frequency will influence the grazing activity of pregnant ewes together with underlying molecular and physiological changes and changes to lamb phenotype (wool quality).

Materials and Methods

Ethics Approval

This research was approved by the Western Sydney University ACEC Protocol (A12610) and ratified by The University of Queensland ACEC Approval Protocol (SAFS/544/19).

Experimental Design

The study was conducted on a sheep property in Cattai NSW 2756, Australia. In January 2019, a total of 100 Merino ewes (mixed age) participated in natural joining. A pre-joining shearing was performed on each ewe to provide baseline samples. Top-knot wool sample provides a convenient method of wool collection in sheep for research purposes [see (19)]. Once joined successfully, 48 Merino ewes were used for this experiment (the experimental ewes were run together with the rest of the dry ewes). Ewes were bought into the pens by the researcher and visually assessed and conditioned scored by the researcher. Upon confirmation that the ewe was sound (not ill, or lame), it proceeded to inclusion in the experiment, if the ewe was unsound, it was not included in the research and removed from the experiment at this point.

The experimental design is shown in Figure 1. There were two treatment groups (once shorn or twice shorn) and 24 ewes per treatment group. All ewes were run as a flock, and shearing frequency was the only main factor which was different between the two groups. Ewes (and lambs) were fitted with light weight battery powered collar tags attached with tri-axial accelerometer to discriminate between grazing, standing, and walking activity [see (21–23)].

FIGURE 1

Figure 1. Illustration summarizing the field sampling design presented with the sampling timeline, groups, and sample types.

On Farm Assessment

Pregnancy Scanning of Ewes

On day 50 (mid-gestation) and day 90 (late gestation), ewes were scanned by experienced operator using a “walk-through” system that required minimal setup and consisted of the operator's crate or crush which ewes enter/exit as they were scanned. The ultrasound technicians measured the lengths of the fetus(es), thus providing the approximate time of gestation. For sampling, the ewes were classified as being in early or late pregnancy on the first scanning event on day 50. Pregnancy is characterized by the presence of a fetus(es) with a heartbeat. Ultrasound technicians can also apply an approximation of length of pregnancy from conception to time of scanning—this is known as early and late from time of conception to scanning.

Collection of Wool Fiber and DNA Samples

Ewes were bought into the pens and visually assessed and conditioned scored by the researcher. As part of the normal shearing regime on farm, a sample of wool was collected from the fleece on the top knot (closest to the skin) as this is an area that can be accessed readily on the animal and recently validated for wool cortisol evaluation by the researchers (19). To minimize undue stress on the sheep, we only used ear tissue for the DNA analysis, and this choice of sampling method (in comparison to blood collection) also goes well with the tagging of sheep. An ear notch tissue (<1 g) was collected from the ewes and lambs using the Allflex Tissue Sampling Unit (TSU) (source: https://www.allflex.global/au/product/tsu-applicator/). Ewes were sampled for ear tissue at pregnancy scanning (day 50), and lambs were ear-notched at weaning. Weaning of lambs occurred at 12 weeks of age and involved castration of males, ear marking, ear tagging, and vaccination.

Laboratory Methods

Hormone Analysis

The wool cortisol concentration in each sample was determined by colormetric analysis using polyclonal anticortisol antiserum (R4866—supplier UC Davis, CA, USA) diluted in ELISA coating buffer [Carbonate-Bicarbonate Buffer capsule (Sigma C-3041) and 100 ml of Milli-Q water, pH 9.6], working dilution 1:15,000. This was followed by reactivity with Horseradish Peroxidase (HRP)–conjugated cortisol label (CJM, UC Davis) diluted 1:80,000, and cortisol standards diluted serially (1.56–400 pg per well). Nunc Maxi-Sorp™ plates (96 wells) were coated with 50 μl of cortisol antibody solution and incubated for a minimum of 12 h at 4°C. Standards, including zeros and nsbs (non-specific binding wells), were prepared serially (two-fold) using 200 μl of standard working stock and 200 μl of assay buffer (39 mM NaH₂PO₄H₂O, 61 mM NaHPO₄, and 15 mM NaCl).

For all assays, 50 μl of standard and (1:10) diluted 90% ethanol extracted wool samples were added to each well, followed by 50 μl of the cortisol HRP. Each plate was loaded in under 10 min. Plates were covered with acetate plate sealer and incubated at room temperature for 2 h. After incubation, plates were washed four times using an automated plate washer (Elx50, BioTek™) with phosphate-buffered saline solution (0.05% Tween 20) and then blotted on paper towel to remove any excess wash solution. Substrate buffer was prepared by combining 1 μl of 30% H₂O₂, 75 μl of 1% tetramethylbenzidine (TMB), and 7.425 μl of 0.1 M acetate citrate acid buffer, pH 6.0 per plate. The TMB substrate was added to each well that contained a standard sample at 50 μl to generate color change. The plates were covered with an acetate plate sealer and left to incubate at room temperature for 15 min. The reaction was stopped with 50 μl of Stop solution (0.5 M H₂SO₄ and Milli-Q water) added to all wells in the Nunc Maxi-Sorp™ plates. To determine hormone concentration in each sample plates were read at 450 nm (reference: 630 nm) on an Elx800 (BioTeck™) microplate reader. Cortisol concentrations were presented as nanogram per gram.

Wool Laserscan Assessment

The wool sampling method is exactly as per the protocols available here (http://www.wooltesters.com.au/faqs).

Parentage Testing

DNA-based parentage determination was done by the Neogen^® laboratory, Gatton, Queensland (source: https://genomics.neogen.com/en/parentage-testing-products). A total of 18 lambs were matched to ewes from the once shorn treatment, and 13 lambs were matched to ewes from the twice shorn treatment.

Molecular Epigenetic Analysis and Bioinformatics

Sheep DNA analysis and quality control were performed using Illumina NovaSeq RRBS (Reduced Representative Bisulfite Sequencing) data of a barcoded 100-bp single end run. RRBS library was produced following the Nuggen's Ovation^® RRBS Methyl-Seq System in the Australian Genome Research Facility. The primary bioinformatics analysis involved quality control, trimming adapter, contamination and low-quality fragments, customized Nuggen's adapter trimming, read mapping, customized Nuggen's post-alignment processing, DMR (differential methylation region) analysis, and annotating differentially methylated genes. Briefly, reads were assessed with FastQC v.0.11.8 and trimmed with Trim Galore v0.5.0. Additional trimming was performed using Nugen's diversity trimming script with default values. Mapping was carried out with Bismark v0.21.0 to methylation-converted Oar3.1 genome. Alignments were performed with Bowtie2 v2.3.4 aligner with default parameters allowing 0 mismatch in a 20-bp seed.

DNA Methylation

EdgeR is a package used to detect and quantify differential methylation of digital RRBS data, that is, counts of reads support methylated and non-methylated cytosines for each locus of a given organism. The library sizes were corrected by the average of the total read counts for the methylated and unmethylated libraries in edgeR. An organizational log linear model was then used to quantify the differential expression between the groups. P-values were adjusted for multiple hypothesis testing using Benjamini–Hochberg method. Nearest transcriptional start site (TSS) was annotated to valid methylation loci using nearest TSS function in edgeR. DMR and all analyzed loci table contained the following fields: Venn diagrams show the numbers of co-upregulated and co-downregulated loci between different comparisons including only loci that have a FDR < 0.05 and logFC > 0.

Raw sequencing reads was firstly processed with TrimGalore to remove adapter/primer, low quality fragments. Trimmed reads then went through a second trimming to remove Nugen's RRBS specific primer. MultiQC report has been provided to show the statistics of clean reads. Note that the duplicate ratio here looks high because reads digested from the same locus of different molecules have the same sequence in RRBS. These reads were barcoded in library preparation to distinguish from true PCR replicates.

Clean reads were then mapped to reference genome (Oar3.1 build) and spike in Lambda DNA. PCR duplicates were removed from the alignments in the following de-duplication step, and identical reads from different molecules with unique barcodes were retained. MultiQC report has been provided to summarize the mapping results.

Non-conversion Rate

Non-methylated Lambda DNA was used as a spike in to estimate the non-conversion rate of bisulfite-treated genomic DNA. The non-conversion rate was computed by dividing the number of reads support methylated cytosine by the total reads for that cytosine on the non-methylated Lambda DNA, as follows:

$\begin{array}{l}\text{Rnc}(\%)=\text{Nmc}/(\text{Nmc}+\text{Nc})\times 100\end{array}$

where Rnc is non-conversion rate, Nmc is the number of reads supporting methylated cytosine, and Nc is the number of reads supporting non-methylated cytosine.

Differentially Methylated Regions

In this workflow, one of the most popular Bioconductor packages edgeR pipeline is used for assessing DMRs in RRBS data. It is based on the negative binomial (NB) distribution, and it models the variation between biological replicates through the NB dispersion parameter. The analysis was restricted to CpG sites that have enough coverage for the methylation level to be measurable in a meaningful way at that site. We require a CpG site to have a total count (both methylated and unmethylated) of at least 10 in every sample before it is considered in the analysis. The number of CpGs kept in the analysis is ensured large enough for our purposes after filtering. edgeR linear models are used to fit the total read count (methylated plus unmethylated) at each genomic locus. Differential methylation is assessed by likelihood ratio tests; so, here, we consider FDR value of < 0.05 as significant DMR.

Statistical Analysis of Other Field and Lab Parameters (Hormone, Body Condition, and Grazing Activity)

Statistical analysis was done to test the hypothesis that (1) parameters (hormone, body condition, or grazing activity) will be significantly varied between the once or twice shorn pregnant merino ewes. First, data were checked for homogeneity of variances and were log-transformed before analysis. Statistical analysis was done using a one-way ANOVA using sample, date of sampling and sheep number as the factors, and log-transformed data as the dependent variables. Level of significance for all statistical analysis was p < 0.05.

Results

Ewe Body Condition Scores

Ewes from both treatments (shearing groups) showed variability in their body condition score. However, the twice shorn ewes maintained their body condition during mid- and late pregnancy on average higher than the once only shorn ewes (Table 1, p < 0.05). Average range of condition scores (2.5–3.5) is typical of healthy ewes.

TABLE 1

Table 1. Median, mean, and standard error values for ewe condition score across once and twice shorn treatments.

Wool Cortisol Profiles

Mean Wool Cortisol Metabolites (WCM) of Merino ewes initially decreased between pre-joining and mid-pregnancy, however increased between April (mid gestation) and May (late gestation) across both treatments. However, once shorn ewes recorded the highest average WCM results (Figure 2, p < 0.05).

FIGURE 2

Figure 2. Mean Wool Cortisol Metabolites (WCM) of Merino ewes shorn once and twice. Error bars represent standard error of individual data set. N = 24 per group.

Grazing Activity

Collar data were analyzed between post-shearing (April 14, 2019) and May 24, 2019 (late gestation) to determine where shearing once (n = 9 sheep with collar data) or twice (n = 11 sheep with collar data) led to a significance change in grazing activity between the treatments (p < 0.05, Figure 3). Statistical analysis showed a significant difference in grazing frequency between once and twice shorn ewes (p < 0.05). Twice shorn sheep were grazing proportionally more than once shorn sheep by >10% during mid-late gestation.

FIGURE 3

Figure 3. Box-whisker plot showing the percentage activity of ewes in the once (1) or twice (2) shearing treatments. G, grazing; S, standing; and W, walking. Level of significant difference between the groups was p < 0.05 for each activity type.

Epigenetic DNA Methylation

Genewise NB generalized linear models (glmLRT) were used to identify differences in methylation between different shearing treatments and pregnancy durations. In total, 22 samples with one shearing treatment and 23 with two shearing treatments were collected in this study. Moreover, seven and 27 samples were Collected from sheep with a short and long pregnancy, respectively.

One locus (Chr20:50404014) was significantly associated with different shearing treatments (once or twice shorn ewes) (FDR = 0.005). This locus is upstream of a protein-coding gene (ENSOARG00000002778.1), which shows similarities to the forkhead box C1 (FOXC1) mRNA using BLAST searches. No significant differences were measured for pregnancy duration (early or late) using glmLRT (FDR > 0.134). A multidimensional scaling (MDS) plot did not separate the ewe samples according to pregnancy duration or shearing treatment (Figure 4).

FIGURE 4

Figure 4. Multidimensional scaling (MDS) plot based on all methylation sites. Each dot represents a ewe sample and is colored by the pregnancy duration (early vs. late) and shearing treatment (once vs. twice). Only the 50 most significant loci (p < 0.05) were used for the MDS analysis.

Differential methylation analysis was repeated using a combination of the shearing treatment and pregnancy status. The comparisons (1) late pregnancy vs. early pregnancy for sheep with one shearing treatment and (2) late pregnancy vs. early pregnancy for sheep with two shearing treatments were carried out to identify Associations between loci and pregnancy duration for sheep with either one or two shearing events.

As shown in Figure 5A, the Venn diagram illustrates the loci that were up/down in both or in each comparison (1) and (2). Eight loci are significantly upregulated, whereas six loci were downregulated between late vs. early pregnancy in sheep that were sheared once or twice. The annotation of these loci is provided as a bar graph in Figure 6A.

FIGURE 5

Figure 5. Venn diagram illustrating the upregulation and downregulation of loci between different comparisons. The Venn diagrams illustrate the core genes for (A) “Once sheared and late pregnancy vs. once sheared and early pregnancy” as well as in “Twice sheared and late pregnancy vs. twice sheared and early pregnancy” (B) Early pregnancy and once sheared vs. early pregnancy and twice sheared as well as late pregnancy and once sheared vs. late pregnancy and twice sheared sheep. Only features with a FDR < 0.05 and logFC > 0 were included.

FIGURE 6

Figure 6. Bar graph plot showing the annotations and significantly different overlapping methylation sites (in reference to the Venn diagrams, Figures 5A,B). (A) A total of 14 (8/6) loci were upregulated or downregulated for comparisons of pregnancy scanning for once or twice sheared sheep. All loci have not a gene symbol (NA, not available). (B) In total, 22 (12/10) loci were upregulated or downregulated between shearing treatment including APOBECC3Z1 gene and U6.

The comparisons (3) once vs. twice sheared for sheep with early pregnancy and (4) once vs. twice sheared for sheep with late pregnancy were calculated to identify associations between loci and shearing treatment for sheep with either early or late pregnancy. The core loci, i.e., loci that are either significantly upregulated or downregulated between different shearing treatments for sheep with either pregnancy duration are visualized in a Venn diagram (Figure 5B) and bar graphs (Figure 6B). In total, 12 loci were upregulated and 10 loci were downregulated.

Lamb Phenotype and Molecular Data

The results in Table 1 show wool characteristics of the lambs that were matched to their ewes using the DNA parentage test. The lambs from the shorn twice cohort of ewes had a visually finer wool (Average micron column), higher average comfort factor of 0.4% (Ave COMFF column), and spinning fineness difference between shearing frequency groups was 0.9 microns (Ave Spin Fineness Column).

A MDS plot did not separate the ewe samples according to pregnancy duration or shearing treatment (Figures 7A,B).

FIGURE 7

Figure 7. (A,B) Multidimensional scaling (MDS) plot based on all methylation sites. Each dot represents a lamb sample and is colored by the pregnancy duration (early vs. late) and shearing treatment (once vs. twice). Only the 50 most significant loci were used for the MDS analysis.

Discussion

This research aimed to study the molecular epigenetics, stress physiology, and behavior of pregnant Merino ewes using a combination of field and lab-based methods to test the original hypothesis that shearing frequency (twice shorn vs. single shorn) could influence the grazing activity of pregnant ewes with underlying changes to molecular epigenetic profiles of the ewes, wool cortisol, and body condition, with flow-on positive benefits of early pregnancy shearing on lamb phenotype (wool quality).

On the basis of the results from the experiment, it was determined that shearing frequency did result in significant molecular epigenetic change in the pregnant ewes. We discovered one locus (Chr20:50404014), which was significantly associated with different shearing treatments (once or twice shorn ewes) (FDR = 0.005). This locus is upstream of a protein-coding gene (ENSOARG00000002778.1), which shows similarities to the FOXC1 mRNA using BLAST searches. There was lower methylation level in once vs. twice shorn ewes, which indicated that the locus is downregulated. FOXC1 gene provides instructions for making a protein that attaches (binds) to specific regions of DNA and regulates the activity of other genes. On the basis of this action, the FOXC1 protein is called a transcription factor. The FOXC1 protein plays a critical role in early development of the eye and normal development of heart, kidneys, and brain. It also plays critical role in adult management of oxidative stress associated with the visual cortex. The FOXC1 gene has been identified in Chinese native sheep breeds in a recent genomic profiling study, as a novel series of vision-related genes (24). The FOXC1 gene belongs to forkhead family of transcription factors characterized by a distinct DNA-binding forkhead domain. These transcription factors are involved in regulating embryonic and ocular development, and FOXC1 mutations are associated with various glaucoma phenotypes, including primary congenital glaucoma and Axenfeld–Rieger syndrome, which characterized by specific ocular anomalies (25). To our knowledge, our research is the first to demonstrate the epigenetic modulation of FOXC1 gene during pregnancy in Merino sheep in a natural grazing study, and this should be explored further in relation to the consequences on embryonic and ocular development of Merino lambs.

We found that no significant differences were measured for pregnancy duration (early or late) using glmLRT (FDR > 0.134). One explanation may be loss of power. By separating the data into early and late, we reduced the number of samples for the comparisons. A MDS plot did not separate the ewe samples according to pregnancy duration or shearing treatment. Differential methylation analysis was repeated using a combination of the shearing treatment and pregnancy status to identify associations between loci and pregnancy or shearing events. We discovered that eight loci were significantly upregulated, whereas six loci were downregulated between late vs. early pregnancy in sheep that were sheared once or twice (all loci are novel and not previously annotated with BLAST search).2, in total, 12 loci were upregulated and 10 loci were downregulated between different shearing treatments for sheep with either pregnancy duration.

The annotated epigenetically modulated genes for Merino ewes as determined using gene cards (see Supplementary Table 1) included the following genes. APOBEC3 gene (this is Apolipoprotein B MRNA Editing Enzyme Catalytic Subunit). U1, U4, and U6—RNA gene affiliated with the snRNA class. Among its related pathways are mRNA splicing—major pathway and RNA transport. U6 is a type of small nuclear RNA (snRNA) and is highly conserved among species. U6 snRNA located at the heart of the spliceosome participates in the processing of mRNA precursors. U6 is very stable because of the combination of small nuclear ribonucleoprotein complexes, a 5′ cap, a 3′U-rich tail, and the capacity for self-and/or U4 hybridization. The half-life value is ~ 24 h. U6 is one of the most widely used internal reference genes for miRNA. U6 has been used as an internal reference gene in renal tissue, cell lines, and peripheral blood mononuclear cells in patients with kidney disease [see (26)]. SNORD16 (Small Nucleolar RNA, C/D Box 16) is an RNA Gene and is affiliated with the snoRNA class. 5S_rRNA is an integral component of the large ribosomal subunit in all known organisms with the exception only of mitochondrial ribosomes of fungi and animals. It is thought to enhance protein synthesis by stabilization of a ribosome structure. CRYBB3 is a protein-coding gene. Diseases associated with CRYBB3 include Cataract 22, Multiple Types, and Cataract 24. Gene Ontology annotations related to this gene include structural constituent of eye lens. An important paralog of this gene is CRYBB2. THRB—The protein encoded by this gene is a nuclear hormone receptor for triiodothyronine. It is one of the several receptors for thyroid hormone and has been shown to mediate the biological activities of thyroid hormone. Knockout studies in mice suggest that the different receptors, while having certain extent of redundancy, may mediate different functions of thyroid hormone. Mutations in this gene are known to be a cause of generalized thyroid hormone resistance, a syndrome characterized by goiter and high levels of circulating thyroid hormone (T3–T4), with normal or slightly elevated thyroid-stimulating hormone. NFKBIA—This gene encodes a member of the NF-κB inhibitor family, which contain multiple ankrin repeat domains. The encoded protein interacts with REL dimers to inhibit NF-κB/REL complexes that are involved in inflammatory responses.

The annotated epigenetically modulated genes for Merino lambs as determined using gene cards (see Supplementary Table 2) included the following genes: 5S_rRNA; IGFBP7; oar-mir-29b-1, U6, B4GALT7, SNORA70, UCP2, CDKN1A, SNORD14, MRCL3, GATD3A, NOG, NPAS4, 7SK, ST3GAL4, BLOC1S4, HADHA, VEGFA, ADIG, RPS25, PDIA5, and PRDX2 (see Table 2). Some of these interesting candidate genes are discussed.

TABLE 2

Table 2. Average raw wool micron, CC%, comfort factor average spinning fineness or average curvature of lambs grouped to ewes (group column).

Insulin-like growth factor–binding protein 7 (IGFBP7) is an interesting gene relevant to pigmentation in merino sheep and located within a 3-Mb window of ovine chromosome 6 (OAR6), in a region that also contains the KIT gene to which piebald was first associated (27). IGFBP7 could be studied further to boost molecular markers for the elimination of contamination of white wool with contaminated fibers. This is possible with availability of IGFBP7 transcripts now available for sheep (28). NOG gene has earlier been identified as one of the candidate inhibitor genes of secondary follicle differentiation in developing Merino fetus in utero (29). 5S rRNA gene is involved in ribosome functioning (30). SNORA70 gene is important in relation to feed intake, body composition, and live weight, and earlier study based on the local adaptation of Creole cattle genome diversity in the tropics reported an association between SNORA70 populations adapted either to cold or heat conditions, indicating its importance for metabolic adaptation in both thermal conditions (31, 32). Our results for epigenetic modulation of SNORA70 gene also confirms the differential expression work done earlier by (32), whereby maternal under- and overnutrition during gestation in Merino ewes affected small RNA species in the offspring lamb including SNORD113 and SNORA70. SNORD14 is highly significant in relation to heat stress (33). These are small, nuclear, noncoding RNAs that are responsible for guiding RNA for post-transcriptional modifications. Further research is warranted to understand how access to adequate maternal nutrition during gestation impacts of fetal development. MicroRNA such as oar-mir-29b-1 is a non-coding RNA that has been shown to play a major role in neurological brain development in rodents and human models through its influence on the expression of DNA methyltransferase Dnmt3a, which is responsible for the catalyzation of CH methylation. Methyl-CH marks have been hypothesized to regulate neuronal diversity by fine-tuning gene transcription across the genome (34). B4GALT7 (beta-1,4-galactosyltransferase 7) has a key role in collagen network formation. B4GALT7 is involved in glycosylation, has a role in connective-tissue disorders, and is related to disturbed fibril organization and proteoglycan synthesis. B4GALT7 is highly expressed in the growth plate, especially in the proliferative zone; mutations have been shown to cause dwarfism in livestock, e.g., horses (35). Uncoupling protein (UCP) 2 is a widely expressed mitochondrial protein whose precise function is still unclear but has been linked to mitochondria-derived reactive oxygen species production. Thus, the chronic absence of UCP2 has the potential to promote persistent reactive oxygen species accumulation and an oxidative stress response as well as impaired glucose stimulated insulin production (36). Myosin regulatory subunit (MRLC3) plays an important role in regulation of both smooth muscle and nonmuscle cell contractile activity via its phosphorylation. It is implicated in cytokinesis, receptor capping, and cell locomotion. MRLC2 expression has been studied in skeletal muscle contraction of goats (37). Neuronal PAS domain–containing protein 4 (Npas4), a member of the PAS family characterized by a conserved basic-helix-loop-helix motif and PAS domain, acts as an inducible immediate early gene activated with minutes of stimulation to regulate the formation of inhibitory synapses. Importantly, the physicochemical properties reveal the neuro-modulatory role of Npas4 in crucial pathways involved in neuronal survival and neural signaling (38). In a context relevant to future research in livestock welfare, studies in rodent model have shown changes in Npas4 expression levels in response to stress such as fear, sleep deprivation, anxiogenic environments, and maternal separation, and the absence of Npas4 has been associated with higher expression levels of cytokines IL-6 and TNF-α and modulation of brain circuits important for cognitive and emotional sense (38). ADIG is also known to be a lipid-producing protein is a novel transcription factor for regulating and controlling the differentiation and proliferation of adipocytes. The sequence of the ADIG gene and protein of human, mouse, rat, pig, and cattle is highly conserved, and it has similar regulatory and biological functions in the mammalian body. ADIG is a newly focused transcription factor in the field of animal husbandry, which regulates the expression of PPARc and participates in the regulation and control of the production of fat cells. It is closely associated with glycolipid metabolism and sugar metabolism, and targeted regulation of ADIG could be applied to better understand the mechanism of IMF deposition in livestock (39). Peroxiredoxin-2, an antioxidant protein has been recently studied in relation to meat quality parameters such as tenderness (40). These candidate genes discovered in this research could establish the basis for stable reference biomarkers of lamb phenotype, which will be a significant boost for welfare programs.

Collectively, on average, ewes that were shorn early in pregnancy demonstrated higher grazing activity, better body condition, and lower stress levels than once shorn ewes. Furthermore, first generation lambs matched to twice shorn ewes expressed visually finer wool with better comfort scores than those crossbreeding lambs that were matched to once shorn ewes. Early shearing of the ewes resulted in improved grazing activity, which could have supported the elevated nutritional plane during mid- to late gestation of the ewes, future research could benefit from the combination of wool cortisol hormone profiling and body condition scores to better evaluate the production characteristics of Merino ewes) [see (21)]. Although we had not included twin bearing ewes from our trial, it is also important to note that twin bearing ewes give birth to lambs with broader wool in comparison to lambs born from single bearing ewes. This is because of the effects of maternal nutritional partitioning on the fetal growth and development of primary and secondary wool follicles (41, 42). Further research should also inspect sire evaluation to determine the potential influence on ram genetics on fiber diameter.

Overall, the research outcomes contribute significant new knowledge to the Australian sheep production industry, and it will be a valuable tool for sheep health and welfare assessments in future.

Significance Statement

This study provides novel data on behavioral, physiological, and molecular epigenetic signatures in Merino ewes during gestation after going through once or twice shearing intervention. We have discovered DNA methylation profiles from ewes and lambs of genes that are directly associated with whole-animal physiology, development, and growth. The baseline data can provide a useful resource for future research in welfare that will benefit sheep and wool production.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Files, further inquiries can be directed to the corresponding author/s.

Ethics Statement

The animal study was reviewed and approved by Western Sydney University ACEC Protocol (A12610) which was later ratified by the University of Queensland ACEC Approval Protocol (SAFS/544/19).

Author Contributions

EN and GS conceptualized this research. GS was the field research manager. EN conducted the laboratory hormone analysis and supervised the project team. AT provided valuable input during manuscript writing. Molecular analysis and bioinformatics were conducted by the Australian Genome Research Facility. All authors contributed to the article and approved the submitted version.

Funding

Funding for this research was obtained from the Australian Wool Innovation (AWI). DF and RS were awarded the Australian Wool Education Trust (AWET) scholarship.

Conflict of Interest

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.

Publisher's Note

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2022.830450/full#supplementary-material

Supplementary Table 1. Annotated epigenetically modulated genes for Merino ewes as determined using gene cards.

Supplementary Table 2. Annotated epigenetically modulated genes for Merino lambs as determined using gene cards.

References

1. Bock J, Wainstock T, Braun K, Segal M. Stress in utero: prenatal programming of brain plasticity and cognition. Biol Psychiatry. (2015) 78:315–26. doi: 10.1016/j.biopsych.2015.02.036

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Lutz PE, Turecki GDNA. DNA methylation and childhood maltreatment: from animal models to human studies. Neuroscience. (2014) 264:142–56. doi: 10.1016/j.neuroscience.2013.07.069

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Suzuki MM, Bird A. DNA methylation landscapes: provocative insights from epigenomics. Nat Rev Genet. (2008) 9:465–76. doi: 10.1038/nrg2341

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Thompson RP, Nilsson E, Skinner MK. Environmental epigenetics and epigenetic inheritance in domestic farm animals. Anim Reprod Sci. (2020) 220:106316. doi: 10.1016/j.anireprosci.2020.106316

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Daxinger L, Whitelaw E. Understanding transgenerational epigenetic inheritance via the gametes in mammals. Nat Rev Genet. (2012) 13:153–62. doi: 10.1038/nrg3188

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Sinclair KD, Rutherford KM, D Wallace, JM Brameld, JM Stöger, R Alberio, et al. (2016). Epigenetics and developmental programming of welfare and production traits in farm animals. Reprod Fertil Dev. 28:1443–78. doi: 10.1071/RD16102

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Crossman CA, Barrett-Lennard LG, Frasier TR. An example of DNA methylation as a means to quantify stress in wildlife using killer whales. Sci Rep. (2021) 11:1–8. doi: 10.1038/s41598-021-96255-1

PubMed Abstract | CrossRef Full Text | Google Scholar

8. González-Recio O, Toro MA, Bach A. Past, present and future of epigenetics applied to livestock breeding. Front Genet. (2015) 6:305. doi: 10.3389/fgene.2015.00305

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Goddard ME, Whitelaw E. The use of epigenetic phenomena for the improvement of sheep and cattle. Front Genet. (2014) 5:247. doi: 10.3389/fgene.2014.00247

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Ross PJ, Sampaio RV. Epigenetic remodeling in preimplantation embryos: cows are not big mice. Anim Reprod. (2018) 15:204–14. doi: 10.21451/1984-3143-AR2018-0068

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Narayan E, Sawyer G, Parisella S. Faecal glucocorticoid metabolites and body temperature in Australian merino ewes (Ovis aries) during summer artificial insemination (AI) program. PloS ONE. (2018) 13:e0191961. doi: 10.1371/journal.pone.0191961

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Sawyer G, Narayan EJ. A review on the influence of climate change on sheep reproduction. In: Narayan E, editor. Comparative Endocrinology of Animals. London: IntechOpen (2019). doi: 10.5772/intechopen.86799

CrossRef Full Text | Google Scholar

13. Oberlander TF, Weinberg J, Papsdorf M, Grunau R, Misri S, Devlin AM. Prenatal exposure to maternal depression, neonatal methylation of human glucocorticoid receptor gene (NR3C1) and infant cortisol stress responses. Epigenetics. (2008) 3:97–106. doi: 10.4161/epi.3.2.6034

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Glover V. Prenatal stress and its effects on the fetus and the child: possible underlying biological mechanisms. In: Perinatal Programming of Neurodevelopment. New York, NY: Springer (2015). p. 269–83

PubMed Abstract | Google Scholar

15. Ord J, Fazeli A, Watt PJ. Long-term effects of the periconception period on embryo epigenetic profile and phenotype: the role of stress and how this effect is mediated. Periconception in Physiol. Med. (2017) 117–35. doi: 10.1007/978-3-319-62414-3_7

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Weaver SJ, Hynd PI, Ralph CR, Edwards JH, Burnard CL, Narayan E, et al. Chronic elevation of plasma cortisol causes differential expression of predominating glucocorticoid in plasma, saliva, fecal, and wool matrices in sheep. Domest Anim Endocrinol. (2021) 74:106503. doi: 10.1016/j.domaniend.2020.106503

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hargreaves AL, Hutson GD. Changes in heart rate, plasma cortisol and haematocrit of sheep during a shearing procedure. Appl Anim Behav Sci. (1990) 26:91–101. doi: 10.1016/0168-1591(90)90090-Z

CrossRef Full Text | Google Scholar

18. Corner RA, Kenyon PR, Stafford JK, West DM, Oliver MH. The effect of mid-pregnancy shearing or yarding stress on ewe post-natal behaviour and the birth weight and post-natal behaviour of their lambs. Livest Sci. (2006) 102:121–9. doi: 10.1016/j.livprodsci.2005.12.001

CrossRef Full Text | Google Scholar

19. Sawyer G, Fox DR, Narayan E. Pre-and post-partum variation in wool cortisol and wool micron in Australian Merino ewe sheep (Ovis aries). PeerJ. (2021) 9:e11288. doi: 10.7717/peerj.11288

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Pehlivan E, Kaliber M, Konca Y, Dellal G. Effect of shearing on some physiological and hormonal parameters in akkaraman sheep. Asian Australas J Anim Sci. (2020) 33:848. doi: 10.5713/ajas.19.0400

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Alvarenga F, Borges I, Palkovi IˇcL, Rodina J, Oddy V, Dobos R. Using a three-axis accelerometer to identify and classify sheep behaviour at pasture. Appl Anim Behav. Sci. (2016) 181:91–9. doi: 10.1016/j.applanim.2016.05.026

CrossRef Full Text | Google Scholar

22. Mclennan K, Skillings E, Rebelo C, Corke M, Moreira M, Morton J, et al. Technical note: validation of an automatic recording system to assess behavioural activity level in sheep (Ovis aries). Small Rumin Res. (2015) 127:92–6. doi: 10.1016/j.smallrumres.2015.04.002

CrossRef Full Text | Google Scholar

23. Umstatter C, Waterhouse A, Holland J. An automated sensor-based method of simple behavioural classification of sheep in extensive systems. Comput Electron Agric. (2008) 64:19–26. doi: 10.1016/j.compag.2008.05.004

CrossRef Full Text | Google Scholar

24. Wang W., Zhang X, Zhou X, Zhang Y, La Y, Zhang Y, et al. Deep genome resequencing reveals artificial and natural selection for visual deterioration, plateau adaptability and high prolificacy in Chinese domestic sheep. Front Genet. (2019) 10:300. doi: 10.3389/fgene.2019.00300

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Nishimura DY, Swiderski RE, Alward WL, Searby CC, Patil SR, Bennet SR, et al. The forkhead transcription factor gene FKHL7 is responsible for glaucoma phenotypes which map to 6p25. Nat Genet. (1998) 19:140–7. doi: 10.1038/493

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Duan ZY, Cai GY, Li JJ, Bu R, Wang N, Yin P, et al. U6 can be used as a housekeeping gene for urinary sediment miRNA studies of IgA nephropathy. Sci Rep. (2018) 8:1–7. doi: 10.1038/s41598-018-29297-7

PubMed Abstract | CrossRef Full Text | Google Scholar

27. García-Gámez E, Reverter A, Whan V, McWilliam SM, Arranz JJ, International Sheep Genomics Consortium, et al. Using regulatory and epistatic networks to extend the findings of a genome scan: identifying the gene drivers of pigmentation in merino sheep. PLoS ONE. (2011) 6:e21158. doi: 10.1371/journal.pone.0021158

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Mingliang Zhou Pinggui Yang Dengjun Wu Xiangyu Zhang. Cloning and Sequence Analysis of IGFBP-7 Gene in Sheep. Hubei Agri Sci. (2015) 4:38–42.

Google Scholar

29. Menzies M, Stockwell S, Brownlee A, Cam G, Ingham A. Gene expression profiles of BMP4, FGF10 and cognate inhibitors, in the skin of foetal Merino sheep, at the time of secondary follicle branching. Exp Dermatol. (2009) 18:877–9. doi: 10.1111/j.1600-0625.2008.00837.x

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Kouvela EC, Gerbanas GV, Xaplanteri MA, Petropoulos AD, Dinos GP, Kalpaxis DL. Changes in the conformation of 5S rRNA cause alterations in principal functions of the ribosomal nanomachine. Nucleic Acids Res. (2007) 35:5108–19. doi: 10.1093/nar/gkm546

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Freitas PH, Wang Y, Yan P, Oliveira HR, Schenkel F S, Zhang Y, et al. Genetic diversity and signatures of selection for thermal stress in cattle and other two bos species adapted to divergent climatic conditions. Front Genet. (2021) 12:102. doi: 10.3389/fgene.2021.604823

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Hoffman ML, Peck KN, Wegrzyn JL, Reed SA, Zinn SA, Govoni KE. Poor maternal nutrition during gestation alters the expression of genes involved in muscle development and metabolism in lambs. J Anim Sci. (2016) 94:3093–99. doi: 10.2527/jas.2016-0570

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Dikmen S, Cole JB, Null DJ, Hansen PJ. Genome-wide association mapping for identification of quantitative trait loci for rectal temperature during heat stress in Holstein cattle. PloS ONE. (2013) 8:e69202. doi: 10.1371/journal.pone.0069202

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Swahari V, Nakamura A, Hollville E, Stroud H, Simon JM, Ptacek TS, et al. MicroRNA-29 is an essential regulator of brain maturation through regulation of CH methylation. Cell Rep. (2021) 35:108946. doi: 10.1016/j.celrep.2021.108946

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Boegheim IJ, Leegwater PA, van Lith HA. Back W. Current insights into the molecular genetic basis of dwarfism in livestock. Vet. (2017) 224:64–75. doi: 10.1016/j.tvjl.2017.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Pi J, Bai Y, Daniel KW, Liu D, Lyght O, Edelstein D, et al. Persistent oxidative stress due to absence of uncoupling protein 2 associated with impaired pancreatic β-cell function. Endocrinology. (2009) 150:3040–8. doi: 10.1210/en.2008-1642

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Xu HG, Xu GY, Wan L, Ma J. Molecular cloning and functional analysis of MRLC2 in tianfu, boer, and chengdu ma goats. Genet Mol Res. (2013) 12:3510–20. doi: 10.4238/2013.March.15.7

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Fu J, Guo O, Zhen Z, Zhen J. Essential functions of the transcription factor Npas4 in neural circuit development, plasticity, and diseases. Front. Neurosci. (2020) 14:1262. doi: 10.3389/fnins.2020.603373

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Liu S, Huang J, Wang X, Ma Y. Transcription factors regulate adipocyte differentiation in beef cattle. Anim Genet. (2020) 51:351–7. doi: 10.1111/age.12931

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Johnson LG, Steadham EM, Huff-Lonergan EJ, Lonergan SM, Johnson LG, Huff-Lonergan E. Partial purification of peroxiredoxin-2 from porcine skeletal muscle. Muscle Biol. (2021) 5:1–9. doi: 10.22175/mmb.12408

CrossRef Full Text | Google Scholar

41. Thompson AN, MB Ferguson, DJ Gordon, GA Kearney, CM Oldham, BL Paganoni. “Improving the nutrition of Merino ewes during pregnancy increases the fleece weight and reduces the fibre diameter of their progeny's wool during their lifetime and these effects can be predicted from the ewe's liveweight profile.” Anim Prod Sci. (2011) 51: 794–804. doi: 10.1071/AN10161

CrossRef Full Text | Google Scholar

42. Tao X, Xu Z. MicroRNA transcriptome in swine small intestine during weaning stress. PloS ONE. (2013) 8:e79343. doi: 10.1371/journal.pone.0079343

PubMed Abstract | CrossRef Full Text | Google Scholar