All animal procedures included in the current study were approved by the Animal Breeding Ethics Committee at the Desert Research Center (DRC) in Egypt (reference number: AB/NO2022) and the Research Ethics Review Committee, Faculty of Agriculture, University of Alexandria, Egypt (Reference: Alex. Agri. 082305307). All methods and protocols in this study comply with the ARRIVE 2.00 guidelines. The experiment did not include animal euthanasia, and all animals were released to the goat herd after the end of the experiment.

The experiment was conducted at the Agriculture Research Station, Animal Production Research Institute, Agricultural Research Center (ARC), Sakha, Kafr El Sheikh, Egypt. A total of 20 lactating goats were selected to represent both breeds, ensuring sufficient statistical power for comparative analysis. The animals comprised Shami goats (SH) with an average body weight of 40.6 ± 1.39 kg (n = 10) and Zaraibi goats (ZA) with an average body weight of 25.3 ± 1.57 kg (n = 10). The experiment lasted 45 days, with a 15-day adaptation period followed by 30 days of data collection. The adaptation period aimed to adapt the animals to individual housing conditions and feeding systems. All the animals used in this study were the offspring from the Sakha Agriculture Research Station goat herd in Kafr El Sheikh, Egypt. The body weight of the animals was recorded at the beginning of the experiment, and the animals were housed individually and had free access to drinking water. All animals in both groups received the same diet throughout the experiment: fresh Egyptian clover (Trifolium alexandrinum) ad libitum and commercial concentrates (2.5% of live weight per head/day). The quantities of offered and refused clover and concentrates feed mixture were estimated daily for each goat throughout the experiment. Subsequently, dry matter intake (DMI) was calculated as the difference between the offered and refused feed. After 20 days, approximately 200 mL of rumen fluid sample was collected from each doe via a stomach tube before morning feeding. The rumen fluid samples were filtered through a two-layer cheesecloth. The pH of rumen samples was measured using a digital pH meter (WPA CD70, ADWA, Szeged, Hungary). Then, the samples were separated into two portions, which were frozen at −20°C for DNA extraction and to analyze rumen ammonia (NH₃-N) and volatile fatty acids (VFA). Milk yield was measured individually on days 20, 35, and 45 using twice-daily hand milking, morning and evening. Approximately 100 mL of representative milk samples were collected to conduct milk chemical composition. At the end of the experiment, all animals were released back to the goat herd without euthanasia.

The gross FE was calculated using the equation FE = milk production (kg/d) divided by DMI (kg/d). Fat-corrected milk for 3.5% fat content was calculated using the formula: milk yield (MY) 3.5% = (0.432 + 0.1625 × % milk fat) × milk yield, kg/d (Sklan et al., 1992). Subsequently, adjusted FE was calculated using the formula: FE = 3.5% fat-corrected milk yield (kg)/dry matter intake (kg). Milk net energy (Milk NE) was calculated using the equation described in National Research Council (NRC) (2001): Milk NE (Mcal of NE_L/d) = Milk Production x (0.0929 * Fat % + 0.0563 * Protein +0.0395 * Lactose %). Finally, FE for lactation was calculated using the formula: Mcal/kg = Milk NE (Mcal /d) / DMI /d.

Diets: Fresh Egyptian clover, concentrate feed mixture, and refused Egyptian clover were analyzed (replicate for every sample”) according to the method of AOAC (1997) to measure dry matter (DM, method 930.15), crude protein (CP, method 954.01), crude fiber (CF, method 962.09), ether extract (EE, method 920.39), and ash contents (method 942.05; Table 1). Nitrogen-free extract (NFE) was estimated by the difference from the sum of the protein, fat, ash, and crude fiber content.

Chemical composition of milk: The percentage of milk protein, fat, lactose, and total solid, in addition to somatic cell count (SCC, cells/ml), were analyzed using a MilkoScan (130 A/SN. Foss Electric, Hilleroed, Denmark) using three aliquots (6 mL/aliquot) for each sample. Before the testing, milk samples were heated to 40°C and homogenized by vortexing for 20 s.

Fatty acids in milk: The total lipid content of milk was extracted using the Folch et al. (1957) method. Briefly, 2 mL of milk sample were transferred into 15 mL screw cap tubes. Then, 6 mL of a mixture containing chloroform and methanol in a ratio of 2:1 were added to the tubes, followed by vortexing for 3 min. Next, 2 mL of deionized water was added to the tubes and vortexed for 3 min. The tubes were then centrifuged for 30 min at 5000xg. The lower phase containing the extracted lipids was transferred to clean tubes and allowed to dry at room temperature. To methylate the fatty acids, sodium methoxide (2 M) was used following the method suggested by Kramer et al. (1997). The extracted lipids dissolved in 1 mL of hexane were mixed with 200 μL of sodium methoxide and vortexed. The mixture was kept for 10 min at room temperature, and then the clear top layer was transferred to gas chromatography (GC) vials. The fatty acid methyl esters (FAME) were analyzed using GC with a mass spectrometer detector (Trac 1,300, Thermo Fisher Scientific, Waltham, United States) and a TG-5MS Zebron capillary column. The FAME compounds were identified using AMDIS software¹ based on their retention times matching the NIST library database.

Rumen fermentation: For volatile fatty acid (VFA) and ammonia analysis, 1 mL of rumen fluid samples were transferred into 1.5 mL Eppendorf tubes. The samples were then acidified by 200 μL of meta-phosphoric acid 25% (w/v) and stored at −20°C for later analysis. Upon thawing, samples were centrifuged at 30,000 × g (15,000 rpm, JA-17 rotor) for 20 min, and then the supernatant was used for VFA and ammonia determination. A 750 μL of supernatant was transferred to GC vials for VFA analysis by injecting 1 μL into gas chromatography (TRACE 1300, Thermo Fisher Scientific, Waltham, United States) using a capillary column (TR-FFAP 30 m × 0.53 mmI D × 0.5 μm). A standard with known concentrations of VFA was used for calibration. The other 250 μL of supernatant was used for ammonia determination colorimetrically using an ammonia assay kit (Biodiagnostic Company, Dokki, Giza, Egypt). The predicted methane was determined using the following equation: Methane yield = 316/propionate +4.4, according to Williams et al. (2019).

DNA extraction and PCR amplification: Total microbial DNA was extracted from 500 μL of rumen sample. Briefly, the sample was centrifuged at 13,000 rpm. According to the manufacturer’s instructions, DNA was extracted from the precipitated solid material using a QIAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany). The quantity and quality of DNA were assessed using agarose gel electrophoresis and a Nanodrop spectrophotometer 2000 (Thermo Scientific, Massachusetts, United States). The V4 region of the 16S rRNA gene was amplified with the primers 515F and 926R using the following PCR conditions: 94°C for 3 min; 35 cycles of 94°C for 45 s, 50°C for 60 s, and 72°C for 90 s; and a final extension at 72°C for 10 min. The rumen archaeal community was studied using primers Ar915aF and Ar1386R, and the PCR amplification was conducted under the following conditions: 95°C for 5 min; 30 cycles of 95°C for 20 s, 55°C for 15 s, 72°C for 5 min, and a final extension at 72°C for 10 min. PCR amplicons were purified and sequenced using the Illumina MiSeq system.

The generated paired-end sequence reads were analyzed using the DADA2 pipeline in the R platform (Callahan et al., 2016). The fastq files of sequence reads were demultiplexed, and their quality was evaluated. Then, the sequences were filtered, trimmed, and dereplicated, followed by merging read 1 and read 2 together to get denoised sequences. Generate denoised amplicon sequence variants (ASVs), and chimeric ASVs were then removed. Taxonomic assignment of ASVs was performed using the “assign taxonomy” and “addSpecies” functions, as well as the SILVA reference database (version 138). Alpha diversity metrics were measured, including observed ASVs, Chao1, Shannon diversity, and inverse Simpson diversity indices.

Additionally, the beta diversity of the bacterial and archaeal communities was determined as principal coordinate analysis (PCoA) using Bray–Curtis dissimilarity, and figures were created using the phyloseq and ggplot R packages. The differences in alpha diversity indices, relative abundances of bacterial and archaeal phyla and dominant bacterial and archaeal genera, feed intake, rumen fermentation parameters, milk yield and composition, and FE were examined using an unpaired t-test at p < 0.05. Function prediction of microbial communities associated with Shami (SH) and Zaraibi (ZA) goat breeds was conducted based on the 16S rRNA data using PICRUSt2 (Douglas et al., 2020) based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The false discovery rate (FDR) method was used for multiple-comparison correction, and p-values below 0.05 were considered significant.

For animal genotyping, blood samples were collected using vacutainer tubes containing EDTA from each animal’s jugular vein. Samples were directly transferred in an icebox to the Molecular Genetics laboratory at DRC for DNA extraction. Genomic DNA was extracted using a Puregen Core Genomic DNA extraction from blood (Qiagen^®, Germany) according to the manufacturer’s protocol. The quantity and quality of extracted DNA were assessed using a Nanodrop spectrophotometer 2000 (Thermo Scientific, Massachusetts, United States). High-quality DNA samples (≥ 50 ng/μL) were used for genotyping at the Research Institute for Farm Animal Biology, Dummerstorf, Germany, using the Illumina®Inc. Goat_IGGC_65K_v2 Infinium HD array (Illumina, San Diego, CA, United States). SNP locations reported in this paper are based on the latest goat genome version of Capra hircus available from the National Center for Biotechnology Information (NCBI) database (ARS1, NCBI). The genotyping BeadChip contained 59,727 SNPs, evenly distributed throughout the entire genome. Genotype calling was performed using GenomeStudio software (Illumina Inc., San Diego, CA, United States) according to the manufacturer’s protocols. The quality control of genotyped SNPs was performed using PLINK v1.9 software (Chang et al., 2015), with the following filtering criteria: (i) SNPs showing significant deviation from Hardy–Weinberg Equilibrium (HWE) were excluded (p < 10⁻⁶); (ii) SNPs with a minor allele frequency (MAF) ≤ 0.01 were removed; (iii) markers with a genotype call rate < 99% and individuals with a call rate < 90% were filtered out. Additionally, SNPs mapped to unknown chromosomal positions or duplicate positions on the same chromosome were excluded to ensure data integrity for subsequent analyses.

Genetic variance explained by markers. The single-step GBLUP implemented in the BLUPF90 family (Legarra et al., 2014) was used to estimate the SNP effects from genomic estimated breeding values (GEBVs) of genotyped animals using the postGSf90 software of the BLUPF90 package (Aguilar et al., 2010). SNP effects were calculated as: û = DZ’ [ZDZ’]⁻¹ â_g, where û is the vector of the SNP effect; D is the diagonal matrix for weighting factors of the SNP effect; Z is the matrix of genotypes, and â_g is the vector of breeding values predicted for genotyped animals. The variance explained by each SNP was calculated as σ² = û² 2p(1 − p), where û is the SNP effect described above, and p is the allele frequency of the SNP (Zhang et al., 2010). The percentage of genetic variance explained by a window segment of 5 adjacent SNPs was calculated as: (Var (a_i)) / (σ²_a) x100%, where a_i is the genetic value of the i-th region that consists of 5 adjacent SNPs, and σ²_a is the total genetic variance (Wang et al., 2014). Manhattan plots of minus log₁₀ of SNP p-values versus chromosomal location were drawn using the qqman package in R (Turner, 2018).

For genome windows above the threshold (i.e., explained >0.1% of the total genetic variance), functional annotation of genes was obtained from BioMart at the Ensembl Genome Browser² (Kinsella et al., 2011). Gene functions and protein domains were identified by the UniProt OMIA (Online Mendelian Inheritance in Animals) and the GeneCards databases (Safran et al., 2021). The genes that overlapped with the identified genomic interval of the candidate genome were considered candidate genes for FE in goats and were enriched using ShinyGO v. 0.77 (Ge et al., 2020) software. The analysis was based on gene ontology (GO; Ashburner et al., 2000) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Kanehisa et al., 2023) against the goat gene set ontologies. The program’s default parameters were applied, and the results were adjusted to a false discovery rate (FDR < 0.05).

Diversity of bacterial community: The total high-quality non-chimeric reads in the SH and ZA goats were 570,677 and 666,870, respectively, and the difference in the average reads number between the goat groups was not significant (p > 0.05; Table 5). Alpha diversity was calculated as observed ASVs, Chao1, Shannon, and inverse Simpson. (Table 5) showed that the ZA goat group had numerically lower alpha diversity metrics without significant differences (p > 0.05). Principal coordinate analysis (PCoA) based on the Bray-Curtis distance was conducted to show the similarity between microbial communities (Figure 1). The results revealed that bacterial communities were clustered separately based on the goat breeds.

Structure of bacterial community: At the phylum level, the bacterial community in the rumen of goats was affiliated with eight bacterial phyla, including Bacteroidota, Cyanobacteria, Desulfobacterota, Fibrobacterota, Firmicutes, Planctomycetota, Spirochaetota, and Verrucomicrobiota, which were shared between all rumen samples (Table 5). In addition, one bacterial phylum, Synergistota, was observed only in the SH goats. The goat breed affected the relative abundance of bacterial groups (Table 5). Bacteroidota and Firmicutes were the most abundant bacterial phyla, representing more than 98% of the bacterial community. The Bacteroidota phylum dominated the bacterial community and accounted for 85.06 and 87.54% in the SH and ZA goats, respectively. This phylum was classified mainly into the families Prevotellaceae, Rikenellaceae, Bacteroidales RF16 group, Bacteroidales p-251-o5, Bacteroidales F082, Muribaculaceae, and Bacteroidales BS11 gut group (Table 6). Except for the family Prevotellaceae, other families showed higher relative abundance in the SH goats. Family Prevotellaceae represented the most abundant family and accounted for 60.38% in the SH goats and 66.57% in the ZA goats with a significant difference (p < 0.05); this family was classified mainly into two genera, Prevotella and Prevotellaceae NK3B31 group, that were higher in the ZA goats. Family Rikenellaceae was dominated by the Rikenellaceae RC9 gut group, which was higher in the SH goats (Table 6).

Phylum Firmicutes accounted for 14.77% of the SH goats and 11.03% of the ZA goats, with a significant difference (p < 0.05; Table 5). This phylum was classified mainly into four families, Ruminococcaceae, Oscillospiraceae, Christensenellaceae, and Lachnospiraceae, which were higher in the SH goats (p < 0.05; Table 6). Family Oscillospiraceae was affiliated mainly with the genus NK4A214 group, family Christensenellaceae was affiliated mainly with the genus Christensenellaceae R-7 group, and family Lachnospiraceae was affiliated mainly with the Lachnospiraceae NK3A20 group that was higher in the ZA goats (p < 0.05) in addition to the genus Butyrivibrio that was higher in the SH goats (Table 6). Phylum Fibrobacterota and Spirochaetota were more abundant in the ZA goats compared to the SH goats (p < 0.05; Table 5).

Diversity of archaeal community: After data processing and the chimeras were removed, a total of 40,313 remained in group SH, and 61,061 reads remained in the ZA goats without significant difference (p > 0.05; Table 7). Alpha diversity, including observed ASVs, Chao1, Shannon, and Invisimpsone, was similar between the goat breeds, and the ZA goats showed numerically higher observed ASVs, Chao1, and Shannon, and the SH goats exhibited higher Invisimpsone (Table 7; p > 0.05). Principal coordinate analysis (PCoA) based on the Bray-Curtis of archaeal community across the goat breeds demonstrated that the samples were clustered based on the animal breeds (Figure 2).

Structure of archaeal community: The taxonomic analysis of the archaeal community revealed that the ASVs were assigned to two phyla, Thermoplasmatota, which was further classified to the family Methanomethylophilaceae, and phylum Euryarchaeota was affiliated with the family Methanobacteriaceae (Table 7). Goat breed affected the archaeal community, and the family Methanomethylophilaceae had two genera, Candidatus Methanomethylophilus, that accounted for 18.45% in group Shami goats (ASH) and 64.60% in Zaraibi goats (AZA) group with a significant difference (p < 0.05), and unclassified Methanomethylophilaceae that was higher in the SH goats compared to AZA. The family Methanobacteriaceae was further categorized into two genera. Methanobrevibacter constituted 61.44% of the archaeal community in the SH goats and 23.04% in the AZA group, showing a significant difference (p < 0.05). Additionally, the genus Methanosphaera was more prevalent in the SH goats compared to AZA (Table 7).

Function prediction of microbial community: The metabolic pathways of rumen microbial communities in Shami (SH) and Zaraibi (ZA) goats were performed using PICRUSt2. Principal components analysis (PCA) of PICRUSt2 function prediction was generated using the relative abundance of metabolic pathways (Figure 3), which revealed that samples were clustered based on animal breeds. The samples of the SH goats revealed significantly higher (q < 0.05) relative abundances of metabolic pathways related to heme biosynthesis-II, L-arginine degradation (ARGORNPROST-PWY), P562-PWY, polyamine biosynthesis I (POLYAMSYN-PWY), ARG + POLYAMINE-SYN, L-histidine degradation I (HISDEG-PWY), biotin biosynthesis-PWY, the glyoxylate cycle (GLYOXYLATE-BYPASS), pyruvate fermentation to acetone (PWY-6588), the glyoxylate cycle and TCA cycle (GLYCOLYSIS-TCA-GLYOX-BYPASS), anaerobic gondoate biosynthesis (PWY-7663), coenzyme B biosynthesis (TCA-GLYOX-BYPASS), and the tricarboxylic acid (TCA) cycle IV (2-oxoglutarate decarboxylase; P105-PWY; Figure 4). Additionally, the samples of group ZA revealed significantly higher (q < 0.05) relative abundances of metabolic pathways related to ascorbate and aldarate metabolism, carbohydrate metabolism, nitrogen metabolism, sulfur metabolism, lactose and galactose degradation (LACTOSECAT-PWY), amino acid metabolism (BRANCHED-CHAIN-AA-SYN-PWY), L-ornithine biosynthesis I (GLUTORN-PWY), glycolysis III (ANAGLYCOLYSIS-PWY), the sucrose degradation pathway (PWY-5384), nucleoside and nucleotide degradation (PWY-5532), GDP-d-glycero-alpha-d-manno-heptose biosynthesis (PWY-6478), and the tricarboxylic acid cycle citric acid cycle (PWY-7254; Figure 4).

Genetic variance explained by markers: A significant genome window was considered to contribute to the genetic variance of FE traits if it accounted for more than 0.1%. A total of 26 genome windows, each explaining >0.1% of the total genetic variance explained by all windows, with a total of 2.32% of genetic variance (Figure 5). These windows were located on chromosomes 1 (6, 110, and 132 Mb), 2 (9 and 116 Mb), 3 (24 Mb), 7 (88 and 94 Mb), 8 (60 Mb), 10 (76 Mb), 11 (94 and 95 Mb), 12 (47 Mb), 13 (53 Mb), 15 (31 Mb), 18 (21, 23, 25, and 36 Mb), 21 (25 and 41 Mb), 22 (49 Mb), 24 (33 Mb), 25 (1 and 22 Mb), and 26 (3 Mb; Table 8; Supplementary file 1).

Gene enrichment analysis: The gene enrichment analysis identified 26 significant (q < 0.05) biological pathways using a list of the identified candidate genes in the current study (Supplementary Excel file S1). The identified candidate genes were significantly enriched in biological pathways related to the regulation of the transport process of many cellular components such as ions, lipids, organic substances, amino acids, and the general transporter activities mechanism (Figure 6).

The rumen bacterial community in goats was dominated by the phyla Bacteroidota, Firmicutes, Spirochaetota, and Fibrobacterota (Table 5), aligning with findings from previous studies on goats (Giger-Reverdin et al., 2020; Luo et al., 2022). The Firmicutes phylum was higher in the SH goats with lower FE, which was also reported in low-efficiency beef cattle by Brooke et al. (2019).

Host genetics significantly influence the rumen microbiome, with approximately 35% of microbial taxa being heritable (Li et al., 2019). This implies that the difference in the microbial communities between the two goat breeds in the current study is expected. In this context, our PCoA analysis (Figure 1) showed that the bacterial communities in SH and ZA breeds were distinct, which agrees with the findings of Paz et al. (2016) in Holstein and Jersey breeds. Similar results were obtained by Mani et al. (2021) in the rumen of Damara and Meatmaster sheep breeds in South Africa. In the current study, ZA showed numerically low alpha bacterial diversity compared to the SH goats (p > 0.05; Table 5), which is similar to other findings in dairy cows (Shabat et al., 2016).

The variation in the composition of rumen microbiota results in changes in rumen fermentation efficiency, wherever the ZA goats had lower ammonia than the SH goats (24.51 vs. 37.10 mmol/L, p < 0.05). Furthermore, the ZA goats showed higher total VFA than the SH goats (82.76 vs. 72.83 mM, p < 0.05) as well as higher acetic and propionic (Table 4; p < 0.05). VFAs represent the main energy source of the host animal (Paz et al., 2016). Nitrogen and carbohydrate metabolism are the main microbial functions in the rumen (Xue et al., 2022). Higher VFA production in ZA was associated with higher fiber-degrading bacteria and higher carbohydrate metabolism pathways (Table 4; Figure 4), which agrees with Tapio et al. (2017) and Xue et al. (2022), indicating that Zaraibi goats have higher efficiency in carbohydrate metabolism. This finding is supported by higher relative abundances of Fibrobacters, Prevotellaceae, and Lachnospiraceae in the rumen of the ZA goats. These bacterial groups have essential roles in complex carbohydrates and protein metabolism (Thoetkiattikul et al., 2013; Paz et al., 2016; Neumann et al., 2018). Xue et al. (2022) reported that pathways related to carbohydrate metabolism could be used as prediction markers to differentiate efficient and inefficient rumen microbiomes or animals, which helps the future selection of high-efficiency animals. These findings were supported by Wang et al. (2023), who reported higher carbohydrate metabolism pathways in goats with higher growth performance. On the other hand, the ZA goats showed higher relative abundances of branched-chain amino acids (BCAA) pathways. The increase in the relative abundance of branched-chain amino acids (BCAA) pathways could refer to a higher concentration of BCAA in the rumen. Zhang H. L. et al. (2013) reported that BCAA stimulates fiber-degrading microorganisms and VFA production.

Furthermore, the increment in acetic acid production was accompanied by an increment in the PWY-7254 pathway (tricarboxylic acid cycle citric acid cycle), which is involved in the synthesis of precursors of acetic acid (De Sales-Millán et al., 2024). Group ZA showed lower rumen ammonia and higher relative abundances of pathways related to nitrogen metabolism, which could indicate the efficient use of nitrogen in microbial protein synthesis and the decline in nitrogen loss (Firkins et al., 2007). Higher synthesis of microbial protein increases the protein supply to the host animals Zhang H. L. et al. (2013). Furthermore, the ZA goats showed higher sulfur metabolism, which is linked to lower methane production (Wu et al., 2021), as sulfate-reducing bacteria are able to compete with methanogens for H₂ in the rumen, which inhibits the methanogenesis (Zhao and Zhao, 2022).

Genus Prevotella is a major player in ruminal metabolism and predominated the rumen microbiome in several ruminant species (Betancur-Murillo et al., 2022). A higher relative abundance of genus Prevotella was associated with the highly efficient breed (i.e., ZA; Tables 2, 6). This finding is supported by Brooke et al. (2019), who reported that the Prevotella is a microbial marker identifying cattle with higher FE. Genus Prevotella includes members that utilize a wide range of substrates, such as hemicellulose and protein (Matsui et al., 2000). Myer et al. (2015) found a higher relative abundance of Prevotella in the gut of steers with high daily gain and intake. De Vadder et al. (2016) indicated that the genus Prevotella produces succinate, which is the precursor to propionate. The propionate synthesis uses hydrogen molecules, reducing its availability for methane production by methanogenic archaea (Betancur-Murillo et al., 2022). Thus, Prevotella has the potential to be used as an anti-methanogenic agent.

The members of the family Lachnospiraceae have fibrolytic and cellulolytic activities and were associated with higher milk production in Holstein cows (Thoetkiattikul et al., 2013; Paz et al., 2016). In this study, the candidate genus Lachnospiraceae NK3A20 was higher in the ZA goats; this genus was reported to produce butyric acid that promotes the development of the rumen (Huang et al., 2021). Boeckaert et al. (2009) indicated that the Lachnospiraceae NK3A20 group plays an important role in the efficiency of energy utilization. Genus Succiniclasticum is one of the dominant bacteria in the rumen of heifers (Liu et al., 2017), and it converts the succinate to propionate (van Gylswyk, 1995). Shabat et al. (2016) explained that cattle with higher propionate-producing rumen bacteria, such as Succiniclasticum and Prevotella, could achieve efficient feed utilization.

The rumen of the SH goat breed had a significantly higher proportion of the genus Christensenellaceae R-7 group, Mogibacterium. Christensenellaceae R-7 group was associated with FE, rumen digestion and absorption of nutrients, fiber digestion, and protein metabolism, and this genus produces acetic and butyric acids (Bach et al., 2019; Huang et al., 2021). Genus Mogibacterium was higher in the low FE group (i.e., SH). A previous study reported that Mogibacterium was associated with higher methane-emitting cattle (Wallace et al., 2015). Furthermore, Mogibacterium cannot degrade carbohydrates for energy, and it had a higher proportion in lower weight gain in steers (McLoughlin et al., 2020).

Furthermore, the rumen of the SH goats showed a numerically higher proportion than that of the family Ruminococcaceae. This family was linked to lower FE in beef cattle (Brooke et al., 2019), which supports the current findings. The function prediction of the rumen microbiome of the SH goats is enriched with hem and biotin synthesis pathways and L-arginine degradation. Hem is an important component of hemoglobin and a cofactor for many enzymes; therefore, it has an important role in oxygen transportation and various physiological processes (Yang et al., 2023). Arginine plays a key role in urea cycle regulation and protein synthesis, as ammonia is the main product of arginine degradation (Chacher et al., 2012). Biotin vitamin acts as a cofactor responsible for carbon dioxide transfer in carboxylases and was enriched in the rumen of cows with high milk yield and protein (Xue et al., 2022).

Rumen archaea: Methane emission from farm animals is one of the main determiners of FE. The composition of the rumen methanogens was associated with the animal breed, methane emission, and FE (Bharanidharan et al., 2018; McLoughlin et al., 2023). Archaea represent 1–2% of the microbial community and are considered the sole producer of methane in the rumen (Rabee et al., 2022; McLoughlin et al., 2023). Methane represents a 2–12% loss in FE to the host (Johnson and Johnson, 1995). PCoA analysis demonstrated (Figure 2) that the diversity of the archaeal community was influenced by the breed. The previous studies indicated that the abundance of methanogenic archaea in the rumen was associated with methane emissions and FE (Wallace et al., 2015; Zhou et al., 2010).

Methanobrevibacter accounted for 61.44% of the archaeal community in the SH goats, significantly more than in the ZA goats (Table 7; p < 0.05), which correlates with the higher predicted methane emissions in the SH goats. Meanwhile, ZA goats demonstrated a higher relative abundance of Candidatus Methanomethylophilus (64.60%) than the SH goats (18.45%). Candidatus Methanomethylophilus is an H₂-dependent methylotrophic methanogen that derives its energy from the metabolism of methanol and methylamine (Borrel et al., 2012). This genus was associated with improved FE, higher daily gain, and lower methane emissions in different sheep breeds and Charolais steers (McLoughlin et al., 2023), which supports the current findings. Methanobrevibacter is the main methane producer in the rumen (Tapio et al., 2017). It uses hydrogen molecules, besides other substrates, such as acetate and formate, to produce methane (Jeyanathan et al., 2011). Bharanidharan et al. (2018) reported that Hanwoo steers (low efficiency) showed higher methane emissions and a higher abundance of Methanobrevibacter compared to Holstein steers (high efficiency), which supports our findings. The same conclusion was reported by Zhou et al. (2009) on beef steers. Previous studies reported that Methanobrevibacter uses acetate in methane production, leading to higher CH₄ energy loss in low-FE animals (Nkrumah et al., 2006; Hegarty et al., 2007). Propionic acid production suppresses methane production by consuming the hydrogen from the rumen ecosystem. Furthermore, sulfur metabolism could reduce the availability of H₂ for methane production (Zhao and Zhao, 2022). These findings explain the relationship between higher efficiency in the ZA goats and lower Methanobrevibacter and higher Prevotella as the main propionic producer in the rumen (Tapio et al., 2017; Bharanidharan et al., 2018). Thus, microbial communities in goats varied according to breed. Higher FE in the ZA goats could be attributed to higher relative abundances of lignocellulolytic bacteria that affect feed utilization positively. Moreover, a lower relative abundance of major methane-producing archaea might decline the methane energy emission (Difford et al., 2018; Gonzalez-Recio et al., 2018).

Variance explained by markers. FE is a quantitative trait that may be affected by a few genes with large or modest effects. Otherwise, it is affected by multiple genes with small effects (Zhang et al., 2020). In our study, we identified 26 genome windows. Each explained >0.1% of the total genetic variance explained by all windows, with a total of 1.89, 1.98, 2.39, and 1.78% of genetic variance for gross FE (G), adjusted FE (AFE), milk net energy (MNE), and FE for lactation (FEL), respectively. This suggests that the trait is likely controlled by multiple SNPs with small effects (Zhang et al., 2020). Additionally, most of the identified genes seem to have pleiotropic effects and should be considered in selective breeding programs (Sallam, 2023). Alternatively, there is a high genetic correlation between the studied traits that were undertaken to represent FE (e.g., G, AFE, MNE, and FEL). Thus, they shared common genome windows with the same positional candidate genes.

Our GWAS analysis identified multiple genomic regions that may be significantly associated with FE in Egyptian goats. The 132 megabases (Mb) on Capra hircus autosome (CHI) 1 was the most interesting genomic region that contributed to FE with approximately 0.26, 0.28, 0.27, and 0.3% of the total genetic variance explained for GFE, AFE, MNE, and FEL, respectively. This genomic region harbored the propionyl-CoA carboxylase subunit beta (PCCB) gene, which was involved in pathways related to the metabolism of water-soluble vitamins. Additionally, it is one of the two subunits of the biotin-dependent propionyl-CoA carboxylase (PCC), a mitochondrial enzyme involved in the catabolism of odd-chain fatty acids, branched-chain amino acids isoleucine, threonine, methionine, valine, and other metabolites (Wongkittichote et al., 2017).

An interesting genomic region was also located at 95 Mb on CHI11, which contributed to the genetic variance of G, AFE, MNE, and FEL with approximately 0.2, 0.16, 0.17, and 0.16%, respectively. This genomic region harbored the Golgin A1 (GOLGA1) gene, which is involved in vesicular trafficking at the Golgi apparatus level and endosome-to-Golgi trafficking (Wong and Munro, 2014), and the GTPase activating protein and VPS9 domains 1 (GAPVD1) gene, which enables GTPase activating protein binding activity and is involved in the regulation of protein transport. Additionally, the RABGAP1 (RAB GTPase activating protein 1) and DENND1A (DENN domain-containing 1A) genes are involved in pathways related to vesicle-mediated transport and trafficking regulation (Yoshimura et al., 2010; Samarelli et al., 2020). Thus, they were putative candidate genes for FE.

Two potential candidate genes were identified on CHI2 (116 Mb): the Metaxin 2 (MTX2) gene, which is involved in the transport of proteins into the mitochondrion (Del Castillo et al., 2021), and the Alkylglycerone Phosphate Synthase (AGPS) gene, which encodes a protein that catalyzes the second step of ether lipid biosynthesis (Chornyi et al., 2023). Similarly, the TMEM241 (transmembrane protein 241) was identified on CHI24 (33 Mb) and was involved in carbohydrate transport and transmembrane transport (Rodríguez et al., 2016). On CHI21 (40 Mb), the SCFD1 (Sec1 Family Domain Containing 1) gene was observed. It is involved in the regulation of protein transport (Hou et al., 2017), in addition to the Adaptor Related Protein Complex 4 Subunit Sigma 1 (AP4S1) gene, which was included in pathways related to transporter activity and obsolete protein transporter activity (Dell'Angelica et al., 1999). On CHI12 (35 Mb), the COMM Domain Containing 6 (COMMD6) gene, which is included in pathways related to protein metabolism, and TBC1D4 (TBC1 Domain Family Member 4), which is considered an insulin transport.

Similarly, the Enoyl-CoA Delta Isomerase 1 (ECI1) gene, involved in fatty acid metabolism pathways, was identified on CHI25, which spans 1 Mb. Important genes associated with fat and obesity metabolism were located on CHI18 (23 Mb), including the Iroquois Homeobox 3 (IRX3) gene, which acts as a regulator of energy metabolism and is linked to pathways involving the FTO obesity variant mechanism. Additionally, the Alpha-Ketoglutarate Dependent Dioxygenase (FTO) gene, which is known as the fat mass and obesity-associated gene, was also identified (Frayling et al., 2007).

On CHI22, which encompasses 49 Mb, the CACNA2D2 (Calcium Voltage-Gated Channel Auxiliary Subunit Alpha2delta 2) gene was found. This gene is involved in the regulation of calcium channels that control the entry of calcium ions into the cell (Bartuzi et al., 2013). Also identified on this chromosome is the TMEM115 gene, which a role in retrograde transport of proteins from the Golgi to the endoplasmic reticulum (Ong et al., 2014).

CHI8 spans 60 Mb and includes the Glucosylceramidase beta 2 (GBA2) gene, which plays a role in carbohydrate transport and metabolism (Woeste and Wachten, 2017). On CHI15, which spans 31 Mb, the STARD10 (StAR-related lipid transfer domain containing 10) gene is mapped and is predicted to be involved in lipid transport (Floris et al., 2019). The Acyl-CoA dehydrogenase family member 8 (ACAD8) gene, encoding a member of the acyl-CoA dehydrogenase family, catalyzes the dehydrogenation of acyl-CoA derivatives in fatty acids and amino acid metabolism (Zhuang et al., 2022). Additionally, the Galactosidase Beta 1 Like 3 (GLB1L3) gene, also located within this region, is predicted to be involved in carbohydrate metabolic processes (Nicoli et al., 2021).

The 36 megabases of CHI18 harbored several putative candidate genes: (1) the NRN1L (Neuritin 1-like) gene, which is involved in protein metabolism pathways (Zhang W. J. et al., 2013), and (2) the Lecithin-Cholesterol acyltransferase (LCAT) gene, which encodes an extracellular cholesterol esterifying enzyme essential for cholesterol transport (Kris-Etherton et al., 2023).

Lastly, multiple members of the SLC (solute carrier family) gene family transporters (e.g., SLC12A4 12, SLC7A6, SLC17A9, SLCO4A1, SLC38A3, and SLC35G2) were identified in most of the significant genome windows. SLC transporters participate in many important physiological functions, such as nutrient supply, and are highly expressed in different organs, including the kidney, brain, liver, gut, and heart (Zhang et al., 2019). Moreover, our gene enrichment analysis showed that the candidate genes identified in the current study were enriched in biological pathways related to the regulation of the transport process of many cellular components such as ions, lipids, organic substances, amino acids, and the general transporter activities mechanism. These findings support our hypothesis of the significant contribution of the animal genotype in FE in Egyptian goats. Thus, the findings of this study indicated that highly efficient animals have a higher representation of animal gene markers and microbial markers related to the metabolism of carbohydrates, protein, and lipids. Therefore, this information could be a useful aid in selective breeding programs that aim at improving FE in goats.