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

In this study, we determined the effect(s) of early shearing on Australian Merino ewes (Ovis aries) and their lambs. To test this research question, we used a suite of field and laboratory methods including GPS collars, wool cortisol, and epigenetic change between ewes and lambs identified using Illumina NovaSeq RRBS. Once shorn ewes (n = 24) were kept on their full fleece throughout the entire gestation period, whereas twice (early) shorn ewes (n = 24) had their wool shorn pre-joining. Top-knot wool sample was taken from ewes during pre-joining, day 50 (mid-gestation), and day 90 (late gestation) for laboratory analysis. Ewes were pregnancy scanned at mid-gestation to determine whether they were early or late parturition (this confirmation is provided by the pregnancy scanner based on fetus size). Top-knot wool sample was also taken from the lambs at weaning for hormone and wool quality testing. Ear tissue was taken from ewes at day 50 (mid-gestation) and from lambs at lamb marking for DNA analysis. Results showed that twice or early shorn ewes grazed 10% higher and maintained stronger body condition than once shorn ewes. Wool cortisol levels were also significantly lower in the early shorn ewes between mid- and late gestation. Lambs bred from twice shorn ewes had on average better visual wool quality parameters in terms of micron, spin finesses, and curvature. For the DNA methylation results, when comparing a group of once sheared with twice sheared ewes, we have discovered one locus (Chr20:50404014) that was significantly differentially methylated [False Discovery Rate (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. To further our understanding of the potential interaction between pregnancy status and shearing frequency of the ewes, we performed further differential methylation analysis using a combination of shearing treatment and pregnancy scanning status. The comparisons (1) late pregnancy vs. early pregnancy for ewes 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. We discovered that 36 gene loci were significantly modulated either between different shearing treatments or late vs. early pregnancy status of ewes. This result suggests that maternal pregnancy and nutritional status during gestation influence DNA methylation. We further investigated DNA methylation in lambs and identified 16 annotated gene loci that showed epigenetic modulation as a result of being born from an early or late stage pregnancy. From the genomics data, we pointed out that ewes go through epigenetic modifications during gestation, and there is a degree of intra-individual variation in the reproductive performance of ewes, which could be due to combination of intrinsic (genetic and physiological) and extrinsic (management and climatic) factors. Collectively, this research provides novel dataset combining physiological, molecular epigenetics, and digital tracking indices that advances our understanding of how Merino ewes respond to shearing frequency, and this information could guide further research on Merino sheep breeding and welfare.

In this study, we determined the effect(s) of early shearing on Australian Merino ewes (Ovis aries) and their lambs. To test this research question, we used a suite of field and laboratory methods including GPS collars, wool cortisol, and epigenetic change between ewes and lambs identified using Illumina NovaSeq RRBS. Once shorn ewes (n = 24) were kept on their full fleece throughout the entire gestation period, whereas twice (early) shorn ewes (n = 24) had their wool shorn pre-joining. Top-knot wool sample was taken from ewes during pre-joining, day 50 (mid-gestation), and day 90 (late gestation) for laboratory analysis. Ewes were pregnancy scanned at mid-gestation to determine whether they were early or late parturition (this confirmation is provided by the pregnancy scanner based on fetus size). Top-knot wool sample was also taken from the lambs at weaning for hormone and wool quality testing. Ear tissue was taken from ewes at day 50 (midgestation) and from lambs at lamb marking for DNA analysis. Results showed that twice or early shorn ewes grazed 10% higher and maintained stronger body condition than once shorn ewes. Wool cortisol levels were also significantly lower in the early shorn ewes between mid-and late gestation. Lambs bred from twice shorn ewes had on average better visual wool quality parameters in terms of micron, spin finesses, and curvature. For the DNA methylation results, when comparing a group of once sheared with twice sheared ewes, we have discovered one locus (Chr20:50404014) that was significantly differentially methylated [False Discovery Rate (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. To further our understanding of the potential interaction between pregnancy status and shearing frequency of the ewes, we performed further differential methylation analysis using a combination of shearing treatment and pregnancy scanning status. The comparisons (1) late pregnancy vs. early pregnancy for ewes 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. We discovered that 36 gene loci were significantly modulated either between different shearing treatments or late vs. early pregnancy status of ewes. This result suggests that maternal pregnancy and nutritional status during gestation influence DNA methylation. We further investigated DNA methylation in lambs and identified 16 annotated gene loci that showed epigenetic modulation as a result of being born from an early or late stage pregnancy. From the genomics data, we pointed out that ewes go through epigenetic modifications during gestation, and there is a degree of intra-individual variation in the reproductive performance of ewes, which could be due to combination of intrinsic (genetic and physiological) and extrinsic (management and climatic) factors. Collectively, this research provides novel dataset combining physiological, molecular epigenetics, and digital tracking indices that advances our understanding of how Merino ewes respond to shearing frequency, and this information could guide further research on Merino sheep breeding and welfare.

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 postparturition (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).

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)(22)(23)].

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.

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 TM 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 2 PO 4 H 2 O, 61 mM NaHPO 4 , 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 TM ) 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 2 O 2 , 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 2 SO 4 and Milli-Q water) added to all wells in the Nunc Maxi-Sorp TM plates. To determine hormone concentration in each sample plates were read at 450 nm (reference: 630 nm) on an Elx800 (BioTeck TM ) 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 R 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 R 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: 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. Means and standard errors are rounded to two decimal places where necessary. * is used to denote significant difference (p < 0.05) in mean condition score of ewes between once and twice shearing events.

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.

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 (

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).

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.

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).
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.   (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.
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).

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 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. 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.
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 underand 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.