The study was conducted at the Universidade Federal de Viçosa (Viçosa, MG, Brazil), located in the Southwest region of Minas Gerais State (20,7604400S, 42,8,608,200 W), Brazil. The experiment was conducted during the summer season, from December 2016 to April 2017. The experiment was conducted following the standard procedures for Animal Care and Handling stated in the Universidade Federal de Viçosa guidelines process number 24/2018. Data on weather conditions were obtained from the Local Weather Station located 1.39 km from the experimental area. The experiment was composed of three periods of 21 days after an adaptation period to the diet and management of 45 days in a grazing system of Guinea grass. Before the trial, all animals received a preventive application of an anti-parasitic (Eprinex Pour-On, Boehringer Ingelhein, São Paulo, São Paulo, Brazil) at 1 mL/10 kg body weight (BW). Sixteen heifers were used in this trial, and they were divided into four groups (n = 4) composed of four Holstein and four Holstein × Gyr (½ Holstein × ½ Gyr) heifers with two different BW assigned in a randomized block design with a repeated measure scheme. The groups averaged 258.6 ± 24.80 kg and 157.1 ± 24.99 kg BW and grazed two separate sets of 16 paddocks each. Blocks were considered a combination of a set of paddocks and animals’ weights and treatments were randomized within each block (set of paddocks 1 had the heavier animals and set of paddocks 2 had the lighter animals).

Animals were managed in a rotational grazing system comprising 32 Guinea grass (Panicum maximum Jacq. cv. Mombaça) paddocks with an 800 m² area, fertilized with 200 kg of N/ha/yr. and 150 kg of K₂O/ha/yr., for 1 d grazing period. Each set of paddocks had a feeding area with a feed bunk, drinking fountains, and 4 m²/animal of shade freely accessed by the animals. Each heifer group received a concentrate supplement (0.5% BW on DM basis) at 1200 h and mineral salt ad libitum. The concentrate composition is presented in Supplementary Table S1, and was manufactured in an animal feed factory at the Universidade Federal de Viçosa.

On the last day of each period, ruminal samples were taken without fasting at 0800 h using the stomach tube technique (Lodge-Ivey et al., 2009; Henderson et al., 2013; Kumar et al., 2015). The initial volume of rumen fluid (~200 mL) was discarded to avoid saliva contamination. A sample of approximately 50 mL was stored at −80°C for DNA extraction.

Genomic DNA was extracted from ruminal samples as described by Stevenson and Weimer (2009), using the phenol/chloroform method with bead beating for mechanical disruption of microbial cells. Samples were processed and sequenced separately for each experimental period. The V4 hypervariable region of the 16S rRNA gene was paired-end sequenced on the Illumina MiSeq platform following the manufacturer’s guidelines (Illumina, Inc., San Diego, California, United States) at the University of Wisconsin-Madison (United States).

After sequencing, raw reads were demultiplexed and adapters were removed in Illumina Miseq software. Afterward, the sequences were filtered, and chimeras were removed, as well as reads shorter than 200 bp or longer than 500 bp. In addition, ambiguous sequences and homopolymers with more than 8 bp were removed. The paired-end sequences were joined, filtered, and cleaned using Mothur v.1.47.0 following MiSeq SOP (Schloss et al., 2009). V4 sequences were aligned using the SILVA v.138 16S rRNA gene reference database (Quast et al., 2012). Lastly, sequences were grouped into operational taxonomic units (OTUs) using uncorrected pairwise distances clustered with the furthest neighbor method based on a similarity cut-off of 97%. OTUs were normalized across samples, using the sample with the lowest number of sequences. For the bioinformatic analyses, sequences of the same treatment collected from three experimental periods were analyzed together and results are representative of the three sampling periods. The sequences obtained for all samples in the present study were submitted to Sequence Read Archive on the National Center for Biotechnology Information¹ under the accession number PRJNA 956552.

One muscle tissue sample was biopsied at d 80 by sampling 1 g of Longissimus dorsi from each animal using local subcutaneous lidocaine (10 mL at 5%). Samples were washed with a sterile saline solution (0.9%) and stored overnight (4°C) in RNA Later (Qiagen, Hilden, North Rhine, Westphalia, Germany) and subsequently at −80°C. DNA was extracted using the method described by Sambrook and Russell (2001).

Primers for 2 exons of gene IL-1β (GenBank AY851162.1; Supplementary Table S2) were designed using Primer-BLAST. PCR was performed with GoTaq® Master Mix (Promega Corporation, Madison, WI, United States) in a Veriti™ thermocycler (Applied Biosystems, Thermo Fisher Scientific, California, United States). The temperature cycles were as follows: initial denaturation at 95°C for 2 min, 40 denaturation cycles at 95°C for 1 min, annealing at 62°C and 63°C for 1 min for primers one and two, respectively, and extension at 73°C for 1.3 min followed by a final extension at 73°C for 5 min. PCR products were purified from agarose gel using the Wizard SV Gel and PCR clean-up System kit (Promega). Amplification products were sequenced using an AB3500 Genetic Analyzer (Applied Biosystems/Thermo Fisher Scientific, California, United States). Amplicons were forward and reverse sequenced in-house on a Sanger 3,500 Genetic Analyzer platform following the manufacturer’s guidelines (Applied Biosystem Thermo Fisher Scientific, San Diego, California, United States).

Amplicon gene sequencing of the IL-1β gene data was processed using the BioEdit Sequence Alignment Editor version 7.0.5.3, and Muscle software version 3.8.31 ² was used to compare the sequences obtained with the reference gene IL-1β sequence (GenBank AY851162.1).

Heifers were weighed on three consecutive days at the beginning and the end of the experiment. Before weighing, heifers were fasted for 12 h with free access to water. After weighing, animals were placed in paddocks until the next fasting period. At the end of each experimental period, tick counts were performed on the left side of the animal’s body according to the methodology described by Wharton and Utech (1970). After tick collection, all animals received preventive treatment with a pour-on anti-parasitic (Eprinex® Pour-On, Boehringer Ingelhein, São Paulo, São Paulo, Brazil) at a dose of 1 mL/10 kg BW.

Blood samples were collected at d 80, by jugular vein puncture using vacuum tubes (BD Vacutainer® SST® II Advance®, São Paulo, Brazil). Blood serum was separated by centrifugation (2,700 × g for 20 min at 4°C) and frozen at −20°C. Urea (K056), glucose (K082), total protein (K031), albumin (K040), and triglycerides (K117) were quantified using an automatic biochemical analyzer (Mindray, BS200E, Shenzhen, China) and Bioclin kits. IGF-1 was analyzed by chemiluminescence using the UniCell Dxl Access Immunoassay System (Beckman Coulter Inc., Brea, United States). Total T3 and total T4 were analyzed using Beckman kits (ref. number 33,830 and 33,800, respectively, Beckman Coulter®, Brea, United States) in the Access 2 Immunoassay System (Beckman Coulter Inc., Brea, United States).

Statistical analysis was performed using the GLIMMIX procedure of SAS (SAS University Edition) using a randomized block design model with repeated measures:

where: μ = general mean, T_i = fixed effect of treatment (breed) i, P_k = fixed effect of period k, T × P_ik = fixed effect of the interaction between treatment i and period k, B_l = random effect of block l (combined effect of weight and paddock), and ɛ_ijklm = random errors with a mean 0 and variance σ 2 , which is the variance between measurements within animals. Fifteen covariance structures were tested for each response variable and the one that provided the lower AIC was used. Data on tick counts and weight exhibited a non-normal distribution; all available distributions in the GLIMMIX procedure were tested using the same model described above and SHIFTED T (3) distribution was used (approached normality of residues).

Data with single observations per animal were analyzed using a randomized block design model:

where: μ = general mean, T_i = fixed effect of treatment (breed) i, B_l = random effect of block l (combined effect of weight and paddock), and ɛ_ijk = random error with a mean 0 and variance σ 2 .

For rumen microbiota statistical analysis, Past v.4.02 (Hammer et al., 2001) software was used to create Non-Metric Multidimensional Scaling (NMDS) plots based on the Bray–Curtis dissimilarity metric. Analysis of similarities (ANOSIM) was performed (number of permutations of 10.000) to evaluate significant differences between experimental groups using the same software. Beta diversity was also evaluated using two additional distance matrices (weighted and unweighted UniFrac) and Principal Coordinate Analysis (PCoA) as the ordination method. Permanova F-statistic and p-values were calculated according to Anderson (2001). Alpha diversity data were analyzed using SPSS v.22 (Norusis, 1993) for normality data using the Kolmogorov–Smirnov test. Subsequently, the Mann–Whitney non-parametric test was performed for independent samples to correlate Holstein × Gyr and Holstein data. Relative abundances of OTUs were analyzed using the Benjamini-Hochberg method to correct the false discovery rate at phylum, family, and genus levels. A two-sided White’s non-parametric t-test was performed using STAMP v.1.2.3 (Parks et al., 2014).

To determine the OTUs most likely to explain the differences between breeds, LEfSe (Linear discriminant analysis Effect Size; Segata et al., 2011) was used to identify potential biomarkers associated with the phenotypes of interest. The statistical analyses were performed in the implementation of LEfSe on the server Galaxy (version 0.1). For this purpose, a standard table of the composition of bacterial OTUs for each breed was provided.

For all analyses, significant differences were declared when p < 0.05, and trends when 0.05 < p < 0.10.

Weather conditions varied along experimental periods, with a minimum temperature of 17.8°C and a maximum of 31.4°C. Relative humidity and Rainfall varied from 75.4% to 79.5% and 2.9 to 4.9 mm, respectively. The temperature-humidity index (THI) ranged from 68.8 to 72. Detailed information on weather data in each period was described by Quirino et al. (2022).

The pre-grazing and post-grazing heights varied during the experiment, with average heights of 71.0 ± 2.91 and 43.8 ± 1.15 cm, respectively. The accumulated herbage per paddock (88.5 ± 5.13) and herbage allowance (11.0 ± 0.63) increased during periods and the major herbage accumulation occurred in periods 3 and 4 with DM of 27.4% and 25.8%, respectively, while grazing efficiency was superior (63.3%) when neutral detergent fiber was inferior (66.0%). Thus, the intake was not limited by forage allowance (Quirino et al., 2022).

After cleaning and filtering the sequencing data, 2,062.650 high-quality sequence reads were obtained in total, with a minimum Good’s coverage of 97% and a maximum value of 99% (Table 1). NMDS visualization did not show a clear separation between the microbial communities of Holstein × Gyr and Holstein heifers (Supplementary Figure 1). However, ANOSIM revealed a significant difference (p < 0.01) between the rumen microbiota composition of each breed (Supplementary Figure 1A). The NMDS analysis carried out separately between sampling days showed similar behavior, except for the first sampling, where no significant difference was observed between the breeds (p = 0.11; Supplementary Figure 1B). A significant difference between breeds was also observed (Permanova, p < 0.05) using weighted and unweighted UniFrac as phylogenetic distance metrics represented in a Principal Coordinate Analysis (PCoA; Figure 1). Alpha diversity metrics (Chao 1, Shannon, and Simpson indexes) were not different between breeds (p = 0.77, p = 0.48, p = 0.32, respectively; Supplementary Figure 2).

The taxonomic assignment of the bacterial sequences revealed a total of 23 phyla, 41 classes, 94 orders, 156 families, and 279 genera. The phylum Firmicutes was in higher abundance for both groups, with a greater relative abundance in Holstein × Gyr heifers (53.1% ± 2.08%) compared to Holstein heifers (51.5% ± 2.13%; Figure 2). The second most abundant phylum was Bacteroidota, which was different between breeds (p = 0.03 with FDR correction). Relative abundances were 26.9% ± 1.92% and 28.7% ± 2.23% for Holstein × Gyr and Holstein heifers, respectively. Unclassified phyla corresponded to 2.24% ± 0.43% for Holstein × Gyr heifers and 2.10% ± 0.52% for Holstein heifers.

At the family level, Prevotellaceae and Lachnospiraceae were the predominant groups in samples from both breeds, but no significant differences were observed (p = 0.25 with FDR correction; Supplementary Figure 3). As for the genus level, there were no significant differences in ruminal samples from both breeds. Genus Prevotella showed the greatest relative abundance in both heifer groups, being 10.6 ± 1.46 for Holstein × Gyr and 11.9 ± 1.95 for Holstein (p = 0.33 with FDR correction; Supplementary Figure 4).

The two breeds shared 4,753 OTUs, while 1,782 OTUs were observed only in the Holstein × Gyr breed and 1,640 OTUs were unique to the Holstein breed (Figure 3). Of the 4,753 OTUs shared between the two breeds, the most abundant (>0.5% relative abundance) are represented in Figure 4. OTU 0002 was found in highest abundance (3.7% relative abundance in the Holstein breed, and 2.3% in the Holstein × Gyr breed) and was classified as a member of the genus Prevotella. Members of the genus Prevotella also represented the most abundant OTUs shared by both breeds (14 OTUs out of 27 OTUs) and differences in relative abundances were detected for this genus between the two animal groups (Figure 4).

The 10 most abundant unique OTUs identified in each breed are listed in Supplementary Table S3, which includes several OTUs of the genus Prevotella. The OTUs classified as Prevotellaceae UCG-003 (0.08% ± 0.37%), Prevotella (0.07% ± 0.16%), and Bacteroidata unclassified (0.05% ± 0.24%) were found exclusively in Holstein × Gyr heifers (Supplementary Table S3). The unique OTUs identified in Holstein heifers were classified as Prevotellaceae unclassified (0.02% ± 0.06%), Lachnospiraceae unclassified (0.01% ± 0.03%), and Prevotella (0.01 ± 0.04%; Supplementary Table S3).

The OTUs that could explain the differences associated with Holstein × Gyr steers through LDA analysis, were OTU00004 (Christensenellaceae_R_7_group), OTU00005 (Prevotella), and OTU00020 (Lachnospiraceae_XPB1014-group). The OTUs that were considered indicators of the microbiota in Holsteins steers were OTU00036 (Prevotella), OTU00016 (NK4A214_group), and OTU00014 (Methanobrevibacter; Figure 5).

Holstein heifers had a greater (p = 0.06; 9.8 ticks/animal) tick count than the Holstein × Gyr heifers (2.56 ticks/animal; Figure 6). Holstein heifers also had greater tick weight than Holstein × Gyr heifers (1.6 g/animal and 0.4 g/animal; p < 0.05). All heifers had more tick infestation in the third period (p < 0.01), resulting in greater tick weight (p < 0.01).

No mutations were identified at exon 1 of the gene IL-1β among 11 sequences analyzed in this study (five Holstein × Gyr and six Holstein heifers). However, mutations were found (Supplementary Figure 5) at Intron I in both breeds. At exon II, we confirmed the mutations described in the reference sequence from GenBank (AY851162.1) for eight of the analyzed sequences (four Holstein × Gyr and four Holstein; Table 2). For Holstein × Gyr heifers the higher presence of nucleotide variation was found in IL-1β at base position 2,126 of exon II. Holstein had the most variation at base position 2,184.

The serum urea content (p = 0.01), IGF-1 (p < 0.01), T3 (p = 0.02), T4 (p = 0.02), and Albumin (p < 0.01) were greater in Holstein × Gyr than in Holstein heifers (Table 3). Breed did not influence serum concentrations of glucose (p = 0.35), total protein (p = 0.25), or triglycerides (p = 0.16).