Regional Heterogeneity and an Unexpectedly High Abundance of Cooperia Punctata in Beef Cattle at Northern Latitude Were Revealed by ITS-2 rDNA Nemabiome Metabarcoding

Background The species composition of cattle gastrointestinal nematode (GIN) communities can vary greatly between regions. Despite this, there is remarkably little large-scale surveillance data for cattle GIN species which is due, at least in part, to a lack of scalable diagnostic tools. This lack of regional GIN species-level data represents a major knowledge gap for evidence-based parasite management and assessing the status and impact of factors such as climate change and anthelmintic drug resistance. Methods This paper presents a large-scale survey of GIN in beef herds across western Canada using ITS-2 rDNA nemabiome metabarcoding. Individual fecal samples were collected from 6–20 randomly selected heifers (n = 1665) from each of 85 herds between September 2016 and February 2017 and 10–25 rst season calves (n = 824) from each of 42 herds between November 2016 and February 2017.

Previous work has shown that the two most abundant GIN species in cattle in Northern temperate regions in North America, including western Canada, are Ostertagia ostertagi and Cooperia oncophora, with a number of other species being present at variable but typically low levels [3,12]. Ostertagia ostertagi is the most pathogenic of the cattle GIN with the potential to cause severe clinical disease, but subclinical infections are more common at northern latitudes resulting in production loss such as reduced average daily gain and carcass quality [1,2,13]. Cooperia oncophora is less pathogenic but has become more concerned in recent years due to the high prevalence of macrocyclic lactone resistance [8,14]. Cooperia punctata, the most pathogenic of the bovine Cooperia species, is typically found at high levels in more southerly regions but was recently reported in several beef cattle herds in eastern Canada, albeit as a minority of the GIN species overall [3].
In this paper, we report the application of ITS-2 nemabiome metabarcoding to survey the GIN communities in 85 beef herds across western Canada. We found an unexpectedly high relative abundance of C. punctata in Manitoba, revealing a distinct regional heterogeneity in GIN communities of cattle reared under relatively similar management and climatic conditions. The occurrence of C. punctata beyond what is generally considered its normal geographical range is of concern given its relatively high pathogenicity and propensity to develop macrocyclic lactone resistance.

Sample population
The cow-calf operations participating in this study were enrolled in the western Canadian Cow-Calf Surveillance Network (WCCCSN). The formation of this network has previously been described [15,16].
Brie y, producers for the WCCCSN were enrolled based on national agricultural census data to obtain a geographically representative sample population from western Canada, including Alberta, Saskatchewan, and Manitoba. Private veterinarians providing services to cow-calf operations in these provinces were asked to assist in recruiting producers. Inclusion criteria for participation were a minimum herd size of 100 cows, willingness to complete questionnaires, and allowing the collection of biological samples from cattle. In situations where more producers were identi ed than needed from a particular region, the required number of producers was selected on a random basis. At the time of data collection for this study, there were 111 cow-calf operations enrolled in the WCCCSN (55 in Alberta, 35 in Saskatchewan, and 21 in Manitoba). Through a survey on parasite management distributed in the summer of 2016, producers were asked to collect fecal samples from their heifers during the fall pregnancy diagnosis.
They were also asked to voluntarily collect fecal samples from their 2016 calf crop around weaning (fall 2016 and winter 2017) and submit them to the laboratory for processing [17].

Fecal sample collection
Fecal sampling from heifers took place between September 2016 and February 2017. Each herd was provided with a sampling kit. The herd veterinarian was asked to collect fecal samples from the rectum of 20 randomly selected heifers based on the availability; if 20 heifers were not available, the youngest bred cows were sampled for a total of 20 samples per herd. This sampling strategy resulted in a median age of sampled heifers of 20 months (range 12-36 months). Calves were sampled from November 2016 to February 2017. Producers were instructed to collect at least two handfuls of feces (freshly voided or rectally collected) from 20 conveniently selected individual calves. The median age of calves sampled was eight months (range 7-9 months). The number of samples collected from each herd ranged from 6 to 20 samples (median 20) for heifers and 10 to 25 (median 20) for calves. Individual fecal samples from heifers and calves were collected into labeled plastic bags, with the air expelled, stored in an insulated container at room temperature, and shipped within 24 hours of collection to the laboratory at the University of Saskatchewan, Saskatoon, Saskatchewan, Canada (heifer samples), or the laboratory at the University of Calgary, Calgary, Alberta, Canada (calf samples). Individual fecal samples were processed within three to ve days of collection. A modi ed Wisconsin sugar otation technique with minor modi cations was used to process fecal samples [18]. Five and 3 g of feces were used to obtain FEC from heifers and calves, respectively. Gastrointestinal nematode eggs were identi ed as strongyle-type spp., Nematodirus spp., or Trichuris spp.

Coproculture larval harvesting
A different pooling strategy was used for the heifer and calf samples because the low number of larvae harvested per herd for the former did not allow meaningful quantitative data to be generated at the individual herd level.

Heifer samples
A modi ed coproculture protocol was used to harvest third stage nematode larvae (L3) [19]. Brie y, a composite fecal sample was prepared for each herd by pooling 12 g of homogenized feces from the individual heifers and three coprocultures set up comprising eighty grams of the composite feces mixed with vermiculite and tap water in a 250 mL glass. Cultures were incubated at room temperature (approximately 20 to 23°C) for 21 days, after which time L3 were harvested, washed twice in tap water by centrifugation at 3725 g for 3 minutes before resuspending in 0.3 mL of tap water, and xed by the addition of 0.7 ml of 95% ethanol. The larvae derived from the three coprocultures per herd composite fecal samples were then pooled to provide one pool of harvested larvae per herd. After the enumeration of an aliquot of larvae by microscopy, samples were stored in the refrigerator (4°C) for several weeks before being sent to the University of Calgary laboratory and kept frozen at -80°C for archiving and ITS-2 rDNA nemabiome metabarcoding. Since there were insu cient larvae obtained for meaningful quantitative analysis at the individual herd level, the L3 harvested from each herd were put into two pools per province; one pool for small herds (≤ 300 cow-calf pairs) and one pool for large herds (> 300 cow-calf pairs). The small herd pools comprised 29, 16, and 10 herds, and the large herd pools comprised 15, 8, and 5 herds for Alberta, Saskatchewan, and Manitoba, respectively. This pooling strategy was chosen to give some degree of replication within the limitations of the small number of larvae harvested from many individual herds and compare the species abundance between herds of different sizes and between provinces. The total number of L3 in the small and large herd pools, respectively were: Alberta = 4,200 and 3,900 L3, Saskatchewan = 3,600 and 4,400 L3, Manitoba = 3,800 and 5,200 L3. Three separate aliquots of 1000 L3 were then taken from each pool to prepare triplicate genomic DNA samples for nemabiome metabarcoding.

Calf samples
For calves, individual animal coprocultures were prepared from 40 g of homogenized feces as per the method of Roberts and O'Sullivan [19]. Culture conditions and harvesting procedures were similar to those described for the heifer samples. Larvae were counted by microscopy, and 50% of the larvae from individual samples were pooled to create a single herd-level pool of larvae. Pooled L3 were washed, centrifuged, enumerated, and xed in 95% ethanol, similar to the procedure described for heifers. Samples were stored at -80°C until being processed for ITS-2 rDNA nemabiome metabarcoding.
There were adequate L3 counts for quantifying GIN species proportions at the herd level for 40 out of 42 herds sampled (n = 22, 10, and 8) in Alberta, Saskatchewan, and Manitoba, respectively). Genomic DNA prepared from 250 L3 from each herd-level pool was used for ITS-2 rDNA sequencing. There were su cient larvae in most herd-level pools for two or three separate aliquots of 250 L3 to prepare duplicate or triplicate genomic DNA samples. However, four herd-level pools had insu cient larvae for duplicate samples; consequently, a single aliquot of 250 L3 was processed (Additional le: Figure 1).
Molecular grade ddH 2 O was used to make 1:10 dilutions of the pooled crude lysates used as a template for rst-round PCR ampli cation of the ITS-2 rDNA target as described in Avramenko et al. [9]. Following puri cation with AMPure XP Magnetic Beads (1X) (Beckman Coulter Inc., Indianapolis, Indiana, USA), Illumina indices and P5/P7 sequencing tags (Illumina inc., San Diego, California, USA) were added using limited cycle PCR ampli cation, and the nal amplicon products puri ed using the same method as above. Approximately 50 ng of amplicon were pooled from each sample to make up the master sequencing library, quanti ed using the KAPA qPCR Library Quanti cation Kit (KAPA Biosystems inc., Wilmington, Massachusetts, USA). The nal concentration of the pooled library was 12.5 nM, with the addition of 25% PhiX Control v3 (Illumina, FC-110-3001), and it was run on an Illumina MiSeq Desktop Sequencer using a 500-cycle pair-end reagent kit (MiSeq Reagent Kits v2, MS-103-2003). Utilizing the Mothur software package, a bioinformatic pipeline was used to assign nematode species identity to each sequenced read using previously described methods [3]. Further details of the pipeline are available at https://www.nemabiome.ca/analysis.html [21]. Sequence reads were multiplied by previously validated correction factors speci c to individual GIN species [9]. The number of sequence reads mapping to each species reference sequence was divided by the total number of mapped reads per sample to determine the percentage species composition of each sample. The total sequence read number mapping to ITS-2 rDNA reference sequences for each sample ranged from 13,401 to 32,801 reads for heifer samples and 15,178 to 83,935 reads for calf samples.

Data analyses
At the individual animal level, the proportions of individual fecal samples (95% con dence interval (CI)) positive for strongyle-type, Nematodirus spp., and Trichuris spp. were determined for calves and heifers based on identifying at least one GIN egg under the microscope. The overall arithmetic mean EPG (± standard deviation (SD)) of the three morphologically different GIN egg types were also calculated for both heifers and calves. For calves, the herd-level arithmetic mean EPG (± SD) of all GIN egg types was calculated and presented with the relative herd-level GIN species proportions.
At the provincial level, the arithmetic mean EPG (± SD) of strongyle-type FEC for each province was calculated, and the statistical difference between provinces was determined using a Generalized Estimating Equations model with a negative binomial distribution and log link function, accounting for clustering at the herd-level for both calves and heifers. Nematodirus spp. and Trichuris spp. FEC were very low in both calves and heifers; therefore, their provincial-level arithmetic means were not estimated.
Alpha diversity was calculated to determine the overall species diversity of GIN populations in calves and heifers within a province. For calves, species diversity data of analytical replicates (i.e., 250 L3 aliquots; 62, 36, and 20 aliquotes for Alberta, Saskatchewan, and Manitoba, respectively) of herd-level pools were used for the mean comparison. For heifers, species diversity data of analytical replicates (i.e., 1000 L3 aliquots; 12 aliquots per province) of each provincial pool were used for the mean comparison. Analytical replicates for heifers were meaningful because they represent multiple herds within a province as larvae from many individual herds went into each provincial pool. The calculations were performed in Mothur v.1.36.1 using the built-in inverse Simpson calculation as previously reported [9]. To assess whether the inverse Simpson index differed signi cantly between each province, a one-way ANOVA, assuming nonequal variances, was performed using a Games-Howell post-hoc comparison in SPSS statistical software (IBM Corp. Released 2012. IBM SPSS Statistics for Macintosh, Version 21.0. Armonk, New York, USA).
Beta diversity estimation was performed for three major species, O. ostertagi, C. oncophora, and C. punctata, using the MetaStats plugin in 275 Mothur v. 1.36.1, using 1000 permutations and default parameters to determine whether the species composition differed in calves and heifers between two different provinces [22]. MetaStats assumed that the data was not normally distributed; therefore, modi ed non-parametric t-tests (two-tailed) for pairwise comparisons of beta diversity estimations of each GIN species in each province were used [22]. This method can overestimate the GIN species with a lower abundance; therefore, signi cance was not claimed if the species were present at less than 2% on average in both comparable groups. Signi cance was declared if P < 0.05. In terms of calves, the arithmetic mean FEC of strongyle-type spp., Nematodirus spp., and Trichuris spp. for each herd were reported with the relative GIN species proportions in each sample.

Results
Fecal egg count and ITS-2 rDNA nemabiome data for gastrointestinal nematodes in heifers ITS-2 rDNA nemabiome metabarcoding revealed that O. ostertagi was the predominant parasite species in heifers from Alberta and Saskatchewan herds (Figure 1), with an overall relative abundance of 53.9 and 52.0% in Alberta small and large herds, respectively, and 59.3 and 61% in Saskatchewan small and large herds, respectively. Cooperia oncophora was the second most abundant GIN in those two provinces, making up 39.4 and 28.6% of the parasite species proportions overall from small and large herds in Alberta, respectively, and 23.5 and 34.9% of the GIN populations from small and large herds in Saskatchewan, respectively. In contrast, Cooperia punctata was the predominant GIN species in Manitoba, with an overall relative abundance of 51.5% and 54.8 in small and large herds, respectively. Ostertagia ostertagi was the second most abundant GIN in Manitoba's small and large herds, with an overall relative abundance of 22.9 and 28.4%, followed by C. oncophora with an overall relative abundance of 10.0 and 10.2% in those herds, respectively. Other GIN species were present in smaller proportions in Alberta, Saskatchewan, and Manitoba regional pools, including Oesophagostomum radiatum (4.1-7.4%), Haemonchus placei (0.9-2.4%), Trichostrongylus longispicularis (0.1-1.6%), Haemonchus contortus (0-0.2%), and Orlo a bisonis (0-0.3%).
The comparison of inverse Simpson index for alpha diversity manifested a signi cant difference in overall GIN species diversity between provinces (F (2, 33) = 27.3, P < 0.001). The post-hoc comparisons of the inverse Simpson indexes revealed that the overall GIN species diversity in heifers was greater in Manitoba (Mean index = 3.60, SD = 0.64) herds than the Alberta (Mean index = 2.54, SD = 0.32) herds (P = 0.037). Similarly, the overall GIN species diversity of Manitoba herds was greater than Saskatchewan herds (Mean index = 2.31, SD = 0.36) herds (P = 0.009). Nevertheless, Alberta heifers' alpha GIN species diversity was not signi cantly different from Saskatchewan heifers (P = 0.765). Beta species diversity comparisons revealed that the relative abundance of C. punctata was signi cantly greater in heifers from Manitoba than Alberta (t = 35.1, P < 0.001) and Saskatchewan (t = 33.6, P < 0.001) ( Table 1) Table 2). The relative abundance of C. punctata was signi cantly lower in Alberta herds than in Saskatchewan herds (t = -3.61, P = 0.001). Gastrointestinal nematode populations in Alberta herds had a signi cantly greater abundance of C. oncophora than GIN populations found on farms from Manitoba (t = 2.50 , P = 0.022) and Saskatchewan (t = 2.05, P = 0.044). Furthermore, O. ostertagi was signi cantly lower in GIN populations found in Manitoba farms compared to GIN populations sampled in Alberta (t = -2.40, P = 0.023) and Saskatchewan farms (t = -3.17, P = 0.003). Statistical differences were not assessed for N. helvetianus, O. bisonis, H. placei, O. radiatum, and T. axei due to their low relative abundance (< 2%).

Discussion
There is a lack of large-scale studies exploring regional differences of GIN species distribution and abundance in cattle, with most existing information being largely historical or anecdotal. This paucity of studies is partly because of the specialist expertise and time required for species identi cation and the lack of scalable tools for quantifying parasite species abundance from fecal samples. The recent development of ITS-2 rDNA nemabiome sequencing now provides a powerful new approach to investigate the relative abundance of GIN species [9]. A recent ITS-2 rDNA metabarcoding study of 50 beef herds across Canada determined that O. ostertagi and C. oncophora were the two most predominant GIN species overall, consistent with the general understanding that these two species are well adapted to cold temperate climates [3,12]. The only aspect of that study that suggested regional differences was the nding that C. punctata was present at a somewhat greater abundance overall in eastern Canada than in the west [3]. However, that species was still much less abundant than O. ostertagi and C. oncophora in the eastern Canadian herds overall. Here we have extended that work to explore potential regional differences in GIN species distribution and abundance in a larger scale study of cow-calf herds across three western Canadian provinces; Alberta (45 herds), Saskatchewan (24 herds), and Manitoba (16 herds).
One of the challenges of working with cattle GIN in northern latitudes is the typically low fecal egg excretion and the consequent low recovery of eggs or larvae from fecal samples, particularly in adult cattle. Consequently, a pooling strategy was taken to provide a regional overview of the GIN species abundance in western Canadian beef operations in which samples from heifers were pooled at the provincial level (two pools per province based on herd size categorization), and samples from calves were pooled at the herd level. The ITS-2 rDNA nemabiome metabarcoding results were consistent with the previous study of Canadian beef calves nding that O. ostertagi was the predominant GIN species in both Alberta and Saskatchewan cow-calf operations, with the next most abundant species being C. oncophora [3]. However, for Manitoba, the results were markedly different from Alberta and Saskatchewan, with a very high abundance of C. punctata across multiple herds. This species was the most abundant GIN species in both the heifer and calf samples from Manitoba overall and was the single most abundant GIN species in calves from 4/8 herds sampled from Manitoba. This nding was surprising because C. punctata is considered poorly adapted to cooler climates and has only been reported to occur at high abundance at more southerly latitudes such as mid-west, southern US, and southern America, where it commonly dominates GIN communities [3,8].
The high abundance of C. punctata at such a northerly latitude is concerning since this species is signi cantly more pathogenic than C. oncophora, the species of Cooperia that generally predominates in cattle in cooler temperate regions such as western Canada. Research in experimentally infected beef calves shows that C. punctata negatively impacts appetite and nutrient utilization resulting in signi cant reductions in average daily gain and dry matter intake [23]. The lack of historical data and detailed regional studies of GIN abundance in North American beef herds over the last few decades means it is not certain whether the high abundance of C. punctata in Manitoba represents a change and, if so, how recently this might have occurred. Nevertheless, this parasite species has not been frequently reported at such northerly latitudes in the past, and to our knowledge, never as the most predominant GIN species. Therefore, there is a strong possibility that these results re ect a range expansion of this species due to several factors. Trichostrongylid GIN are potentially sensitive to climate change since the part of the life cycle outside the host is very temperature-and moisture-dependent [24]. Range expansion associated with climate change has been suggested for several parasitic nematode species in the strongylid group. For example, the ruminant GIN nematode Haemonchus contortus is currently much more common in western Canadian and UK sheep farms than in previous decades, and modeling studies suggest that signi cant range expansion is likely to occur in northern Europe in the coming decades [25,26]. A major concern is the lengthening of the grazing season predicted in western Canada associated with global warming, which could lead to range expansion of GIN more usually found in southern latitudes. Changes in cattle GIN parasite distribution could also be associated with the selection pressure for routinely used anthelmintics, particularly macrocyclic lactones. Ivermectin resistance is now widespread in Cooperia spp. in North America, and a study in the USA in 2015 suggested that the overall prevalence of Cooperia spp. was increasing due to the selection pressure for macrocyclic lactone resistance [8,14]. Whatever the explanation of the unexpectedly high relative abundance of C. punctata in Manitoba cow-calf herds, its regional distribution in western Canada is not simply related to the latitude (Figure 3). The localized distribution in Manitoba may suggest a single or small number of introduction events into that region, perhaps associated with the importation of cattle from particular sources further south. It will be interesting to map the distribution of C. punctata in Manitoba in more detail and investigate potential risk factors for its presence and the history of animal movements and importations in the region. It is also noteworthy that the sampling timing was not markedly different between provinces and so not the reason for the differences in the GIN community composition. However, all samples from heifers and calves were collected in the fall and winter. Therefore, a more detailed temporal study is warranted both across the grazing season and between years to understand further the relative abundance and the regional distribution of this parasite.
Cooperia punctata is also one GIN species that most commonly develop resistance to macrocyclic lactones, with several USA cattle reports [27,28]. Ivermectin is currently the most widely used macrocyclic lactone in western Canadian beef operations [17]. We have recently con rmed ivermectin resistance in C. oncophora, H. placei, and C. punctata and provided evidence for ivermectin resistance of hypobiotic larvae of O. ostertagi in western Canada using integrated ITS-2 rDNA nemabiome metabarcoding and FEC reduction test results [29]. This nding suggests the emergence of ivermectin resistance in multiple GIN species in western Canada, including C. punctata.

Conclusion
In summary, this study has revealed striking regional differences in GIN species abundance in western Canada. Speci cally, a previously unidenti ed high relative abundance of C. punctata in beef cattle in Manitoba. The lack of previous reports of this parasite species at such a northerly latitude raises the possibility that its range is geographically expanding in North America, perhaps under the in uence of climate change, anthelmintic use, animal movement, and/or other as yet undetermined management factors. The emergence of C. punctata as a major constituent of cattle parasite communities in northern latitudes is potentially problematic given the high pathogenicity and capacity for anthelmintic resistance in this GIN species. This study also provides an excellent illustration of the value of ITS-2 rDNA nemabiome metabarcoding as a surveillance tool for ruminant GIN parasites.   Figure 1 Fecal egg counts and gastrointestinal nematode species proportions from heifer fecal sample pools from western Canadian beef farms. Samples were pooled at the provincial level with two pools per province categorized based on herd size (small ≤ 300 cow-calf pairs, large > 300 cow-calf pairs). Three independent aliquots (R1, R2, R3) of 1000 third stage larvae were taken for ITS-2 rDNA nemabiome metabarcoding for each pool. Figure A represents the arithmetic mean of individual animal strongyle-type fecal egg counts of each provincial pool. Figure B indicates the relative gastrointestinal nematode species proportions determined by ITS-2 rDNA nemabiome metabarcoding of third-stage larvae pooled by herd size and province.

Figure 2
Fecal egg counts and gastrointestinal nematode species proportions in fecal samples of calves from 40 western Canadian cow-calf operations. Fecal samples from 10-20 individual calves were pooled at the herd level. Figure A represents the arithmetic means of 10-20 individual calf fecal egg counts of strongyle-type, Nematodirus spp. and Trichuris spp. eggs in each herd. Figure B represents the GIN species proportions in each determined by ITS-2 rDNA nemabiome metabarcoding of herd-level thirdstage larval pools. Numbers from 1-40 on the X-axis identify the individual herds. Except for herds 9, 14, 25, and 34, nemabiome metabarcoding was undertaken on duplicate or triplicate aliquots of 250 larvae per herd (Additional le 1: Figure 1).

Figure 3
Approximate locations of 40 cow-calf operations where ITS-2 rDNA nemabiome metabarcoding data in fecal samples were available for calves. Those herds with a predominance of Cooperia punctata are indicated in green. Each location identi ed by an arrow contains two farms with the same postal code.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.